text:
using jInv.Mesh;
using FactoredEikonalFastMarching;

include("../examples/getAnalyticalMediums.jl");

function runSensitivityTest()
tic()
println("########### 2D TEST SETTING ##############################");
n = [100,100]
Omega = [0.0,4.0,0.0,8.0];
Msh = getRegularMesh(Omega,n-1);

(kappaSquared,src,T_exact) = getAnalyticalConstGradInv2D(Msh.n+1,Msh.h);
T0 = getAnalytic2DeikonalSolution(n,Msh.h,src)[1];
println("First order test");
doSensTest(kappaSquared,src,T_exact,Msh,T0,false);
println("Second order test");
doSensTest(kappaSquared,src,T_exact,Msh,T0,true);
##########################################################
println("########### 3D TEST SETTING ##############################");
n = [64,64,32]
Omega = [0.0,1.0,0.0,1.0,0.0,0.5];
Msh = getRegularMesh(Omega,n-1);
(kappaSquared,src,T_exact) = getAnalyticalConstGradInv3D(Msh.n+1,Msh.h);
T0 = getAnalytic3DeikonalSolution(n,Msh.h,src)[1];
println("First order test");
doSensTest(kappaSquared,src,T_exact,Msh,T0,false);
println("Second order test");
doSensTest(kappaSquared,src,T_exact,Msh,T0,true);

##########################################################
toc()
return;
end

function doSensTest(kappaSquared,src,T_exact,Msh,T0,HO)
n = Msh.n+1;
pMem = getEikonalTempMemory(n);
pEik = getEikonalParam(Msh,kappaSquared,src,HO);
solveFastMarchingUpwindGrad(pEik,pMem);
T = T0[:].*pEik.T1[:];
println("Maximum error norm: ",maximum(abs((T_exact[:] - T))));
println("Sensitivity Check:");
delta_kappa = 0.05*randn(size(kappaSquared))*mean(kappaSquared[:]);
Tlin = zeros(prod(n));
for k = 1:10
	delta_kappa *= 0.5;
	new_kappa = kappaSquared + delta_kappa;
	pEikNew = getEikonalParam(Msh,new_kappa,src,HO);
	solveFastMarchingUpwindGrad(pEikNew,pMem);
	getSensMatVecEik(delta_kappa[:],Tlin,pEik);
	Tlin += T;
	print("||dm||: ");
	print(norm(delta_kappa[:]))
	print(". ||T(m + dm) - T(m)|| ");
	Tnew = T0[:].*pEikNew.T1[:];
	err0 = (T - Tnew);
	print(norm(err0));
	print(". ||T(m + dm) - (T(m) + J*dm)||: ");
	err = Tnew - Tlin;
	println(norm(err));

end
v = rand(prod(n));
Jv = zeros(prod(n));
getSensMatVecEik(v,Jv,pEik);
Z = rand(size(Jv));
JtZ = rand(prod(n));
getSensTMatVecEik(Z,JtZ,pEik);
if abs(dot(JtZ[:],v) - dot(Jv[:],Z[:])) > 1e-8
	println("Checking that (J^T)^T = J for eikonal solver falied: ",abs(dot(JtZ[:],v) - dot(Jv[:],Z[:])) );
end
return;
end


text:
function solve_2sum(xs, target, target_ix)
    negative_target_indices = Dict{Int64, Int64}()
    for (ix, item) in enumerate(xs[target_ix+1:end])
        if haskey(negative_target_indices, item)
            jy = negative_target_indices[item]
            return [jy, ix + target_ix, target_ix]
        end

        push!(negative_target_indices, Pair(target - item, ix + target_ix))
    end
end

function solve_3sum_problem(xs)
    for (ix, item) in enumerate(xs)
        indices = solve_2sum(xs, -item + 2020, ix)
        if indices !== nothing
            return prod(map(k -> xs[k], indices))
        end
    end
end

function main()
    lines = readlines("../resources/input_01.txt")
    xs = map(line -> parse(Int, line), lines)
    result = solve_3sum_problem(xs)
    println(result)
end

main()

ext: jl


text:
@hierarchical HierarchicalMatrix LowRankMatrix Matrix

import LinearAlgebra: matprod
function (*)(H::AbstractHierarchicalMatrix{T}, x::AbstractVector{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    mul!(zeros(TS, size(H, 1)), H, x)
end
function (*)(H::AbstractHierarchicalMatrix{T}, x::AbstractMatrix{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    mul!(zeros(TS, size(H, 1), size(x, 2)), H, x)
end

LinearAlgebra.mul!(y::AbstractVecOrMat, H::AbstractHierarchicalMatrix, x::AbstractVecOrMat) = mul!(y, H, x, 1, 1)
LinearAlgebra.rmul!(H::AbstractHierarchicalMatrix, b::Diagonal) = scale!(H, b.diag, 1)
LinearAlgebra.lmul!(b::Diagonal, H::AbstractHierarchicalMatrix) = scale!(b.diag, H, 1)

mul!(y::AbstractVecOrMat, H::AbstractHierarchicalMatrix, x::AbstractVecOrMat, istart::Int, jstart::Int) = mul!(y, H, x, istart, jstart, 1, 1)

add_col!(H::AbstractHierarchicalMatrix, u::Vector, j::Int) = add_col!(H, u, 1, j)

@generated function mul!(y::AbstractVecOrMat, H::HierarchicalMatrix, x::AbstractVecOrMat, istart::Int, jstart::Int, INCX::Int, INCY::Int)
    L = length(fieldnames(H))-1
    T = fieldname(H, 1)
    str = "
    begin
        M, N = blocksize(H)
        p = 0
        for m = 1:M
            q = 0
            for n = 1:N
                Hmn = H.assigned[m,n]
                if Hmn == 1
                    mul!(y, getindex(H.$T, m, n), x, istart + p, jstart + q, INCX, INCY)"
    for l in 2:L
        T = fieldname(H, l)
        str *= "
                elseif Hmn == $l
                    mul!(y, getindex(H.$T, m, n), x, istart + p, jstart + q, INCX, INCY)"
    end
    str *= "
                end
                q += INCX*blocksize(H, 1, n, 2)
            end
            p += INCY*blocksize(H, m, N, 1)
        end
        return y
    end"
    return Meta.parse(str)
end

@generated function scale!(H::HierarchicalMatrix, b::AbstractVector, jstart::Int)
    L = length(fieldnames(H))-1
    T = fieldname(H, 1)
    str = "
    begin
        M, N = blocksize(H)
        q = 0
        for n = 1:N
            for m = 1:M
                Hmn = H.assigned[m,n]
                if Hmn == 1
                    scale!(getindex(H.$T, m, n), b, jstart + q)"
    for l in 2:L
        T = fieldname(H, l)
        str *= "
                elseif Hmn == $l
                    scale!(getindex(H.$T, m, n), getindex(H.$T, m, n), b, jstart + q)"
    end
    str *= "
                end
            end
            q += blocksize(H, 1, n, 2)
        end
        return H
    end"
    return Meta.parse(str)
end

@generated function scale!(b::AbstractVector, H::HierarchicalMatrix, istart::Int)
    L = length(fieldnames(H))-1
    T = fieldname(H, 1)
    str = "
    begin
        M, N = blocksize(H)
        p = 0
        for m = 1:M
            for n = 1:N
                Hmn = H.assigned[m,n]
                if Hmn == 1
                    scale!(b, getindex(H.$T, m, n), istart + p)"
    for l in 2:L
        T = fieldname(H, l)
        str *= "
                elseif Hmn == $l
                    scale!(getindex(H.$T, m, n), b, getindex(H.$T, m, n), istart + p)"
    end
    str *= "
                end
            end
            p += blocksize(H, m, N, 1)
        end
        return H
    end"
    return Meta.parse(str)
end

@generated function add_col!(H::HierarchicalMatrix{S}, u::Vector{S}, istart::Int, j::Int) where S
    L = length(fieldnames(H))-1
    T = fieldname(H, 1)
    str = "
    begin
        H1 = deepcopy(H)
        M, N = blocksize(H)
        n = 1
        q = blocksize(H, 1, n, 2)
        while q < j
            n += 1
            q += blocksize(H, 1, n, 2)
        end
        q -= blocksize(H, 1, n, 2)
        p = 0
        for m = 1:M
            Hmn = H.assigned[m,n]
            if Hmn == 1
                add_col!(getindex(H.$T, m, n), u, istart + p, j - q)"
    for l in 2:L
        T = fieldname(H, l)
        str *= "
            elseif Hmn == $l
                if isimmutable(getindex(H.$T, m, n))
                    H[Block(m), Block(n)] = add_col(getindex(H.$T, m, n), u, istart + p, j - q)
                else
                    add_col!(getindex(H.$T, m, n), u, istart + p, j - q)
                end"
    end
    str *= "
            end
            p += blocksize(H, m, N, 1)
        end
        return H
    end"
    return Meta.parse(str)
end


text:
module StandardizedMatrices
export StandardizedMatrix

#-----------------------------------------------------------------------------# types
typealias AVec{T} AbstractVector{T}
typealias AMat{T} AbstractMatrix{T}
typealias VecF Vector{Float64}

#----------------------------------------------------------------# StandardizedMatrix
"""
Treat a matrix as standardized, similar to `z = StatsBase.zscore(x, 1)`, without
altering the original data.

`z = StandardizedMatrix(x::AbstractMatrix)`
`z = StandardizedMatrix(x::AbstractMatrix, μ::AbstractVector, σ::AbstractVector)`
"""
immutable StandardizedMatrix{T <: AbstractFloat, S <: AMat} <: AMat{T}
	data::S
	μ::VecF			# column means
	σinv::VecF		# inverse of column stdevs
end
# This acts as inner constructor.  All constructors should call this.  See:
# http://docs.julialang.org/en/latest/manual/constructors/#parametric-constructors
function StandardizedMatrix(x::AMat, μ::AVec, σ::AVec)
	T = eltype((x[1] - μ[1]) / σ[1])
	n, p = size(x)
	μp = length(μ)
	σp = length(σ)
	@assert p == μp "size(x, 2) == $p doesn't match length(μ) == $μp"
	@assert p == σp "size(x, 2) == $p doesn't match length(σ) == $σp"
	StandardizedMatrix{T, typeof(x)}(x, μ, one(eltype(σ)) ./ σ)
end
StandardizedMatrix(x::AMat) = StandardizedMatrix(x, vec(mean(x, 1)), vec(std(x, 1)))

#----------------------------------------------------------------------# Base methods
# http://docs.julialang.org/en/release-0.4/manual/interfaces/#man-interfaces-abstractarray
Base.size(o::StandardizedMatrix) 				= size(o.data)
Base.linearindexing(o::StandardizedMatrix)		= Base.linearindexing(o.data)
function Base.getindex(o::StandardizedMatrix, i::Int, j::Int)
	v = getindex(o.data, i, j)
	return (v - o.μ[j]) * o.σinv[j]
end
function Base.getindex(o::StandardizedMatrix, i::Int)
	v = getindex(o.data, i)
	j = ceil(Int, i / size(o, 1))
	return (v - o.μ[j]) * o.σinv[j]
end
function Base.:*{T <: Real}(A::StandardizedMatrix, B::AVec{T})
	y = zeros(typeof(A[1] * B[1]), size(A, 1))
	A_mul_B!(y, A, B)
	y
end
function Base.:*{T <: Real}(A::StandardizedMatrix, B::AMat{T})
	y = zeros(typeof(A[1] * B[1]), size(A, 1), size(B, 2))
	A_mul_B!(y, A, B)
	y
end

#------------------------------------------------------# Matrix-Vector multiplication
function Base.A_mul_B!(y::AVec, A::StandardizedMatrix, b::AVec)
	A_mul_B!(y, A.data, Diagonal(A.σinv) * b)
	center!(y)
end
function Base.At_mul_B!(y::AVec, A::StandardizedMatrix, b::AVec)
	At_mul_B!(y, A.data, b - mean(b))
	_scale!(y, A.σinv)
end

#------------------------------------------------------# Matrix-Matrix multiplication
function Base.A_mul_B!(y::AMat, A::StandardizedMatrix, b::AMat)
	A_mul_B!(y, A.data, Diagonal(A.σinv) * b)
	center!(y)
end
function Base.At_mul_B!(y::AMat, A::StandardizedMatrix, b::AMat)
	At_mul_B!(y, A.data, b .- mean(b, 1))
	_scale!(y, A.σinv)
end

#---------------------------------------------------------------------------# helpers
function center!(v::Vector)
	μ = mean(v)
	for i in eachindex(v)
		@inbounds v[i] -= μ
	end
	v
end
function center!(m::Matrix)
	μ = mean(m, 1)
	@inbounds for j in 1:size(m, 2)
		μj = μ[j]
		for i in 1:size(m, 1)
			m[i, j] -= μj
		end
	end
	m
end
function _scale!(a::Vector, b::Vector)
	for i in eachindex(a)
		@inbounds a[i] *= b[i]
	end
	a
end
function _scale!(a::Matrix, b::Vector)
	for j in 1:size(a, 2)
		for i in 1:size(a, 1)
			a[i, j] *= b[i]
		end
	end
	a
end

end #module


text:
# Exercise 7-1
function printn(s, n)
    while n > 0
        println(s)
        n -= 1
    end
end

# printn("Repeat", 5)

# Exercise 7-2
const eps = 1e-5

function mysqrt(a)
    x = a / 2
    while true
        y = (x + a/x) / 2
        if abs(y - x) < eps
            break
        end
        x = y
    end
    return x
end

function testsquareroot()
    println("a\tmysqrt\tsqrt\tdiff")
    println("-\t------\t----\t----")
    for i = 1:10
        myval, jval = mysqrt(i), sqrt(i)
        println("$i,\t$myval,\t$jval,\t$(myval - jval)")
    end
end

# testsquareroot()

# Exercise 7-3
function evalloop()
    println("Enter the next expression:")
    result = nothing
    while true
        str = readline()
        if str == "done"
            return result
        end
        result = eval(Meta.parse(str))
        println(result)
    end
    return result
end

# println(evalloop())

# Exercise 7-4

function estimatepi()
    coeff = 2*sqrt(2)/9801

    function term(k)
        return factorial(4*k) * (1103 + 26390*k) / factorial(k)^4 / 396^(4*k)
    end

    sum = curr = term(0)
    k = 1
    while coeff*curr < 1e-15
        curr = term(k)
        sum += curr
        k += 1
    end
    return 1 / (coeff*sum)
end

println(estimatepi())


text:
module FiniteDifferenceDerivatives

using ArrayViews

include("specialized.jl")

export fdd!, fdd, fddmatrix, fddat

# returns a sparse matrix D such that D*u is a der-derivative of specified order
function fddmatrix{T<:Number}(x::AbstractVector{T},der::Int,order::Int)
    npts = length(x)
    m = zeros(T,npts,npts)
    f = eye(T,npts)
    for i = 1:npts
        m[:,i]=fdd(f[:,i],x,der,order)
    end
    return sparse(m)
end

# compute the der-derivative using order-point scheme
function fdd!{T<:Number}(df::AbstractVector{T},f::AbstractVector{T},x::AbstractVector{T},der::Int,order::Int)
    npts = length(x)
    if order < der
        error("Order can not be smaller than der")
    elseif npts < order
        error("Number of mesh points shouldn't be smaller than order")
    end

    # specialized implementation for lower orders
    if der == 1 && order == 3
        fdd13!(df,f,x)
    elseif der == 2 && order == 3
        fdd23!(df,f,x)
    elseif der == 1
        fdd1n!(df,f,x,order)
    else
        # proceed with general implementation

        c = Array(T,order,der+1)

        for N = 1:npts
            # N1 is the leftmost point of the stencil
            N1    = min(max(1,N-div(order-1,2)),npts-order+1)
            df[N] = fddcore(f,x,der,x[N],N1,order,c)
        end
    end

end

function fdd{T<:Number}(f::AbstractVector{T},x::AbstractVector{T},der::Int,order::Int)
    df = Array(T,length(f))
    fdd!(df,f,x,der,order)
    return df
end

function fddcore{T<:Number}(f::AbstractVector{T},
                            x::AbstractVector{T},
                            der::Int,
                            x0::Real,
                            N1::Int,
                            order::Int,
                            c::Matrix{T})

    generatec!(c,x0,x,N1)

    dfN = zero(T)
    @simd for j = 1:order
        @inbounds dfN += c[j,der+1]*f[N1-1+j]
    end
    return dfN
end

function fddat{T<:Number}(f::AbstractVector{T},
                          x::AbstractVector{T},
                          der::Int,
                          x0::Real)

    if length(f) != length(x)
        error("x and f must have the same size")
    end

    N1 = 1
    order = length(x)
    c = Array(T,order,der+1)
    return fddcore(f,x,der,x0,N1,order,c)
end

# generate the coefficients c
function generatec!{T}(c::Matrix{T},x0::Real,x::AbstractVector{T},N1::Int)
    order = size(c,1)
    der   = size(c,2)-1

    if N1 < 1 || N1+order-1 > length(x)
        error("N1=$N1 and order=$order are out of bounds")
    end

    c1 = one(T)
    @inbounds c4 = x[N1] - x0
    c[:] = zero(T)
    c[1] = one(T)
    for i=1:order-1
        mn = min(i,der)
        c2 = one(T)
        c5 = c4
        @inbounds c4 = x[i+N1] - x0
        j = 0
        while j <= i-1
            @inbounds c3 = x[i+N1] - x[j+N1]
            c2 = c2*c3
            if j == i-1
                s = mn
                while s >= 1
                    @inbounds c[i+1,s+1] = c1*(s*c[i,s] - c5*c[i,s+1])/c2
                    s-=1
                end
                @inbounds c[i+1,1] = -c1*c5*c[i,1]/c2
            end
            s = mn
            while s >= 1
                @inbounds c[j+1,s+1] = (c4*c[j+1,s+1] - s*c[j+1,s])/c3
                s-=1
            end
            @inbounds c[j+1,1] = c4*c[j+1,1]/c3
            j+=1
        end
        c1 = c2
    end
end

end # module


text:
using HorizonSideRobots

function stair_marking(r::Robot)
    level = 0
    x = show_path_after_moving(r, West)
    y = show_path_after_moving(r, Sud)

    width = show_path_after_moving(r, Ost)
    while isborder(r, West)
        move!(r, West)
    end

    while isborder(r, Nord) == false
        n = level
        while isborder(r, Ost) == false
            move!(r, Ost)
        end
        while isborder(r, West) == false
            if n <= 0
                putmarker!(r)
            end
            move!(r, West)
            n -= 1
        end
        if level - width <= 0
            putmarker!(r)
        end
        move!(r, Nord)
        level += 1
    end

    while isborder(r, Sud) == false
        move!(r, Ost)

    end

    move_n(r, Nord, y)
    move_n(r, Ost, x)
end

function show_path_after_moving(r::Robot, side::HorizonSide)
    n = 0
    while isborder(r, side)==false
        move!(r, side)
        n = n + 1
    end
    return n
end

function move_n(r::Robot, side::HorizonSide, n::Int)
    for i in 1:n
        move!(r, side)
    end
end


text:
# https://github.com/bokeh/Bokeh.jl

supportedArgs(::BokehBackend) = [
    # :annotation,
    # :axis,
    # :background_color,
    :linecolor,
    # :color_palette,
    # :fillrange,
    # :fillcolor,
    # :fillalpha,
    # :foreground_color,
    :group,
    # :label,
    # :layout,
    # :legend,
    :seriescolor, :seriesalpha,
    :linestyle,
    :linetype,
    :linewidth,
    # :linealpha,
    :markershape,
    :markercolor,
    :markersize,
    # :markeralpha,
    # :markerstrokewidth,
    # :markerstrokecolor,
    # :markerstrokestyle,
    # :n,
    # :bins,
    # :nc,
    # :nr,
    # :pos,
    # :smooth,
    # :show,
    :size,
    :title,
    # :windowtitle,
    :x,
    # :xlabel,
    # :xlims,
    # :xticks,
    :y,
    # :ylabel,
    # :ylims,
    # :yrightlabel,
    # :yticks,
    # :xscale,
    # :yscale,
    # :xflip,
    # :yflip,
    # :z,
    # :tickfont,
    # :guidefont,
    # :legendfont,
    # :grid,
    # :surface,
    # :levels,
  ]
supportedAxes(::BokehBackend) = [:auto, :left]
supportedTypes(::BokehBackend) = [:none, :path, :scatter] #,:steppre, :steppost, :sticks, :hist2d, :hexbin, :hist, :bar, :hline, :vline, :contour]
supportedStyles(::BokehBackend) = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
supportedMarkers(::BokehBackend) = [:none, :auto, :ellipse, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross, :star5] #vcat(_allMarkers, Shape)
supportedScales(::BokehBackend) = [:identity, :ln] #, :ln, :log2, :log10, :asinh, :sqrt]
subplotSupported(::BokehBackend) = false

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

function _initialize_backend(::BokehBackend; kw...)
  @eval begin
    warn("Bokeh is no longer supported... many features will likely be broken.")
    import Bokeh
    export Bokeh
  end
end

# make255(x) = round(Int, 255 * x)

# function bokehcolor(c::Colorant)
#   @sprintf("rgba(%d, %d, %d, %1.3f)", [make255(f(c)) for f in [red,green,blue]]..., alpha(c))
# end
# bokehcolor(cs::ColorScheme) = bokehcolor(getColor(cs))

const _glyphtypes = KW(
    :ellipse    => :Circle,
    :rect       => :Square,
    :diamond    => :Diamond,
    :utriangle  => :Triangle,
    :dtriangle  => :InvertedTriangle,
    # :pentagon   =>
    # :hexagon    =>
    # :heptagon   =>
    # :octagon    =>
    :cross      => :Cross,
    :xcross     => :X,
    :star5      => :Asterisk,
  )

function bokeh_glyph_type(d::KW)
  lt = d[:linetype]
  mt = d[:markershape]
  if lt == :scatter && mt == :none
    mt = :ellipse
  end

  # if we have a marker, use that
  if lt == :scatter || mt != :none
    return _glyphtypes[mt]
  end

  # otherwise return a line
  return :Line
end

function get_stroke_vector(linestyle::Symbol)
  dash = 12
  dot = 3
  gap = 2
  linestyle == :solid && return Int[]
  linestyle == :dash && return Int[dash, gap]
  linestyle == :dot && return Int[dot, gap]
  linestyle == :dashdot && return Int[dash, gap, dot, gap]
  linestyle == :dashdotdot && return Int[dash, gap, dot, gap, dot, gap]
  error("unsupported linestyle: ", linestyle)
end

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

function _create_plot(pkg::BokehBackend, d::KW)
  # TODO: create the window/canvas/context that is the plot within the backend (call it `o`)
  # TODO: initialize the plot... title, xlabel, bgcolor, etc

  datacolumns = Bokeh.BokehDataSet[]
  tools = Bokeh.tools()
  filename = tempname() * ".html"
  title = d[:title]
  w, h = d[:size]
  xaxis_type = d[:xscale] == :log10 ? :log : :auto
  yaxis_type = d[:yscale] == :log10 ? :log : :auto
  # legend = d[:legend] ? xxxx : nothing
  legend = nothing
  extra_args = KW()  # TODO: we'll put extra settings (xlim, etc) here
  bplt = Bokeh.Plot(datacolumns, tools, filename, title, w, h, xaxis_type, yaxis_type, legend) #, extra_args)

  Plot(bplt, pkg, 0, d, KW[])
end

function _add_series(::BokehBackend, plt::Plot, d::KW)
  bdata = Dict{Symbol, Vector}(:x => collect(d[:x]), :y => collect(d[:y]))

  glyph = Bokeh.Bokehjs.Glyph(
      glyphtype = bokeh_glyph_type(d),
      linecolor = webcolor(d[:linecolor]),  # shape's stroke or line color
      linewidth = d[:linewidth],          # shape's stroke width or line width
      fillcolor = webcolor(d[:markercolor]),
      size      = ceil(Int, d[:markersize] * 2.5),  # magic number 2.5 to keep in same scale as other backends
      dash      = get_stroke_vector(d[:linestyle])
    )

  legend = nothing  # TODO
  push!(plt.o.datacolumns, Bokeh.BokehDataSet(bdata, glyph, legend))

  push!(plt.seriesargs, d)
  plt
end

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

# TODO: override this to update plot items (title, xlabel, etc) after creation
function _update_plot(plt::Plot{BokehBackend}, d::KW)
end

function _update_plot_pos_size(plt::AbstractPlot{BokehBackend}, d::KW)
end

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

# accessors for x/y data

# function getxy(plt::Plot{BokehBackend}, i::Int)
#   series = plt.o.datacolumns[i].data
#   series[:x], series[:y]
# end
#
# function setxy!(plt::Plot{BokehBackend}, xy::Tuple{X,Y}, i::Integer)
#   series = plt.o.datacolumns[i].data
#   series[:x], series[:y] = xy
#   plt
# end

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

function _add_annotations{X,Y,V}(plt::Plot{BokehBackend}, anns::AVec{@compat(Tuple{X,Y,V})})
  for ann in anns
    # TODO: add the annotation to the plot
  end
end

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

function _create_subplot(subplt::Subplot{BokehBackend}, isbefore::Bool)
  # TODO: build the underlying Subplot object.  this is where you might layout the panes within a GUI window, for example

end

function _expand_limits(lims, plt::Plot{BokehBackend}, isx::Bool)
  # TODO: call expand limits for each plot data
end

function _remove_axis(plt::Plot{BokehBackend}, isx::Bool)
  # TODO: if plot is inner subplot, might need to remove ticks or axis labels
end

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

function Base.writemime(io::IO, ::MIME"image/png", plt::AbstractPlot{BokehBackend})
  # TODO: write a png to io
  warn("mime png not implemented")
end

function Base.display(::PlotsDisplay, plt::Plot{BokehBackend})
  Bokeh.showplot(plt.o)
end

function Base.display(::PlotsDisplay, plt::Subplot{BokehBackend})
  # TODO: display/show the subplot
end


text:
using MFCC, MFCC.Mel, MFCC.Filterbanks, Plots

windows, windowstarts = mffilters(10, 22050, 22050)
filterbank = Filterbank(22050, windows, windowstarts)
melscale = mel.(range(1, length=22050, 22050))

# Plot the curve of the hertz scale by the mel scale
melhz = plot(melscale, xlabel="Frequency(Hz)", ylabel="Frequency(Mel)")
savefig(melhz, "docs/assets/mel_hz_mapping.svg")

# Plot the filterbank in the mel scale
melfb = plot(melscale, getfilters(filterbank), xlabel="Frequency(Mel)")
savefig(melfb, "docs/assets/mel_scale_filterbank.svg")

# Plot the filterbank in the hertz scale
hzfb = plot(getfilters(filterbank), xlabel="Frequency(Hz)")
savefig(hzfb, "docs/assets/hz_scale_filterbank.svg")


text:
# Simple Send and Accept with Arguments

using SimLynx
SimLynx.greet()

using Distributions: Exponential, Uniform
using Random

@process lock_process() begin
    while true
        @accept caller lock(x::String, n) begin
            println(x)
            println(n)
        end
        @accept caller unlock()
    end
end

the_lock = Nothing

@process p1(i::Integer) begin
    println("$(current_time()): Process p1($i) started")
    @send the_lock lock("message from $i", i)
    println("$(current_time()): Process p1($i) acquired lock")
    work(rand(Uniform(1.0, 10.0)))
    println("$(current_time()): Process p1($i) released lock")
    @send the_lock unlock()
    println("$(current_time()): Process p1($i) ended")
end

function run_simulation(n::Integer)
    @simulation begin
        # current_trace!(true)
        global the_lock = @schedule now lock_process()
        for i = 1:n
            @schedule at rand(Uniform(0.0, 10.0)) p1(i)
        end
        start_simulation()
    end
end

run_simulation(10)


text:
using StanVariational
using Test

if haskey(ENV, "JULIA_CMDSTAN_HOME")

  bernoulli_model = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
  ";

  bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

  stanmodel = VariationalModel("bernoulli", bernoulli_model)

  rc = stan_variational(stanmodel; data=bernoulli_data)

  if success(rc)

    @testset "Bernoulli variational example" begin
      # Read sample summary (in ChainDataFrame format)
      sdf = read_summary(stanmodel)

      @test sdf[sdf.parameters .== :theta, :mean][1] ≈ 0.32 atol=0.1
    end

  end

else
  println("\nJULIA_CMDSTAN_HOME not set. Skipping tests")
end


text:
using Mocha
srand(12345678)

############################################################
# Prepare Random Data
############################################################
N = 10000
M = 20
P = 10

X = rand(M, N)
W = rand(M, P)
B = rand(P, 1)

Y = (W'*X .+ B)
Y = Y + 0.01*randn(size(Y))

############################################################
# Define network
############################################################
backend = DefaultBackend()
init(backend)

data_layer = MemoryDataLayer(batch_size=500, data=Array[X,Y])
weight_layer = InnerProductLayer(name="ip",output_dim=P, tops=[:pred], bottoms=[:data])
loss_layer = SquareLossLayer(name="loss", bottoms=[:pred, :label])

net = Net("TEST", backend, [loss_layer, weight_layer, data_layer])
println(net)

############################################################
# Solve
############################################################
lr_policy = LRPolicy.Staged(
  (6000, LRPolicy.Fixed(0.001)),
  (4000, LRPolicy.Fixed(0.0001)),
)
method = SGD()
params = make_solver_parameters(method, regu_coef=0.0005, mom_policy=MomPolicy.Fixed(0.9), max_iter=10000, lr_policy=lr_policy)
solver = Solver(method, params)
add_coffee_break(solver, TrainingSummary(), every_n_iter=100)

solve(solver, net)

learned_b = similar(B)
copy!(learned_b, net.states[2].b)

#println("$(learned_b)")
#println("$(B)")

Mocha.dump_statistics(solver.coffee_lounge, get_layer_state(net, "loss"), true)

shutdown(backend)


text:
using HTTP, Sockets

todos = """
ToDo 1: Getting groceries
ToDo 2: Visiting my therapist
ToDo 3: Getting a haircut
"""
const HOST = ip"127.0.0.1"
const PORT = 9999
const ROUTER = HTTP.Router()                    # 1
HTTP.@register(ROUTER, "GET", "/*", req -> HTTP.Response(200, "Hello"))  # 2
HTTP.@register(ROUTER, "GET", "/list_todos", req -> HTTP.Response(200, todos))
HTTP.serve(ROUTER, HOST, PORT)                  # 3


text:
# To be able to test the language model, i.e., context and word embeddings
using JLD, Knet

const language_model = "/ai/data/nlp/conll17/competition/chmodel_converted/english_chmodel.jld"
const data_file = "/ai/data/nlp/conll17/ud-treebanks-v2.0/UD_English/en-ud-train.conllu"

include("../src/types.jl")
include("../src/lm_utils.jl")
include("../src/preprocess.jl")

function test_load_conllu(filename)
    d = load(language_model); v = create_vocab(d); wmodel = makewmodel(d);
    corpus = load_conllu(filename, v)
    ppl = fillvecs!(wmodel, corpus, v)
    unk = unkrate(corpus)
    @assert (ppl == 43.70551674466947); info("perplexity test passed!");
    @assert (unk[1]/unk[2] == 0.0938094190678691); info("unkrate test passed!");
end
!isinteractive() && test_load_conllu(data_file)

# perplexity=43.70551674466947 unkrate=0.0938094190678691


text:
module Moiety_Test

using Test, PoreMatMod

# function for a specific test case. NOT a generally useful function!
# frag and moty inputs are expected to be permuted sets of atoms with unique species
function test_is_equiv(frag::Crystal, moty::Crystal)::Bool
    for (f, f_label) in enumerate(frag.atoms.species)
        for (m, m_label) in enumerate(moty.atoms.species)
            if f_label == m_label
                for i in 1:3
                    if frag.atoms.coords.xf[i,f] ≠ moty.atoms.coords.xf[i,m]
                        return false # detected an atom w/ inconsistent coordinates
                    end
                end
            end
        end
    end
    return true
end

@testset "Moiety Tests" begin
    moiety_bcfm = moiety("S-bromochlorofluoromethane.xyz")
    moiety_2!bcfm = moiety("!-S-bromochlorofluoromethane.xyz")
    @test moiety_bcfm.atoms.species == [:C, :Cl, :F, :Br, :H]
    @test moiety_2!bcfm.atoms.species == [:C, :Cl, :F, :Br, :H!]
end

end # module


text:
# Model and fit the continuum

struct ContinuumModel
    material::Material
    e0::Float64
    takeoff::Float64
    mc::Type{<:MatrixCorrection}
    br::Type{<:NeXLBremsstrahlung}
    """
        ContinuumModel(mat::Material, e0::Float64, det::DetectorEfficiency, takeoff::Float64)

    Create a continuum model for the specified material, beam energy, detector and take-off angle.  Computes
    the *detected* quantity of continuum generated in the sample.
    """
    ContinuumModel(
        material::Material,
        e0::Float64,
        takeoff::Float64;
        matrixcorrection::Type{<:MatrixCorrection} = Riveros1993,
        bremsstrahlung::Type{<:NeXLBremsstrahlung} = Castellano2004b,
    ) = new(material, e0, takeoff, matrixcorrection, bremsstrahlung)
end

"""
    emitted(cm::ContinuumModel, ea::Float64)

Compute the intensity of the measured continuum emitted from the material and conditions specified in the continuum
model object at the specified measured energy `ea`.
"""
function emitted(cm::ContinuumModel, ea::Float64) #
    g = generated(cm, ea)
    return g > 0.0 ?
           g * correctcontinuum(
        continuumcorrection(cm.mc, cm.material, ea, cm.e0),
        cm.takeoff,
    ) : 0.0
end

"""
    generated(cm::ContinuumModel, ea::Float64)

Compute the intensity of the measured continuum generated from the material and conditions specified in the continuum
model object at the specified measured energy `ea`.
"""
generated(cm::ContinuumModel, ea::Float64) = #
    ea <= cm.e0 ? bremsstrahlung(cm.br, ea, cm.e0, cm.material) : 0.0

"""
    fitcontinuum(
      spec::Spectrum,
      resp::AbstractArray,
      rois::Vector{UnitRange};
      brem::Type{<:NeXLBremsstrahlung} = Castellano2004a,
      mc::Type{<:MatricCorrection} = Riveros1993,
    )

    Fit a continuum model to the specified range of channels (`rois`).  The `resp` argument is a matrix which describes
the detector response on a channel-by-channel basis.  It can be calculated from an `EDSDetector` and an
`DetectorEfficiency` using `resp = NeXLSpectrum.detectorresponse(det, eff)`.  The `Spectrum` object must have
the :Composition, :BeamEnergy and :TakeOffAngle properties defined.
"""
function fitcontinuum(
    spec::Spectrum,
    resp::AbstractArray{<:Real,2},
    rois::AbstractArray{<:UnitRange,1};
    brem::Type{<:NeXLBremsstrahlung} = Castellano2004a,
    mc::Type{<:MatrixCorrection} = Riveros1993,
)
    @assert haskey(spec, :Composition) "The fitcontinuum(...) function requires the spec[:Composition] property."
    @assert haskey(spec, :BeamEnergy) "The fitcontinuum(...) function requires the spec[:BeamEnergy] property."
    @assert haskey(spec, :TakeOffAngle) "The fitcontinuum(...) function requires the spec[:TakeOffAngle] property."
    cmod = ContinuumModel(
        spec[:Composition],
        spec[:BeamEnergy],
        spec[:TakeOffAngle],
        matrixcorrection = mc,
        bremsstrahlung = brem,
    )
    minEmod = max(50.0, energy(lld(spec), spec)) # Ensures
    model = resp * map(e -> e > minEmod ? emitted(cmod, e) : 0.0, energyscale(spec))
    k =
        sum(sum(model[roi]) for roi in rois) /
        sum(sum(counts(spec, roi, Float64)) for roi in rois)
    props = Dict{Symbol,Any}(
        :TakeOffAngle => spec[:TakeOffAngle],
        :BeamEnergy => spec[:BeamEnergy],
        :Composition => spec[:Composition],
        :K => k,
        :Name => "Brem[Global][$(spec[:Name])]",
    )
    return Spectrum(spec.energy, model / k, props)
end

"""
    continuumrois(elems, det::EDSDetector, minE::Float64, maxE::Float64)

Compute the ROIs for the contiguous continuum regions for the specified elements `elems` on an
`EDSDetector` for the specified range of energies.
"""
function continuumrois(
    elems,
    det::EDSDetector,
    minE::Float64,
    maxE::Float64,
)::Vector{UnitRange{Int}}
    # Calculate the un-rois
    function invert(rois, over)
        st, res = over.start, []
        for roi in rois
            low, hgh = max(st, over.start), min(roi.start, over.stop)
            if hgh > low
                push!(res, low:hgh)
            end
            st = roi.stop
        end
        if st < over.stop
            push!(res, st:over.stop)
        end
        return res
    end
    extend(roi, n) = roi.start-n:roi.stop+n
    # Average channel width between minE and maxE.
    avgwidth = let
        minC, maxC = channel(minE, det), channel(maxE, det)
        (energy(maxC, det) - energy(minC, det)) / (maxC - minC)
    end
    extra = round(Int, (2.0 * resolution(enx"Mn K-L3", det)) / avgwidth)
    rois = mapreduce(elm -> extend.(extents(elm, det, 1.0e-3), extra), append!, elems)
    # Join rois into contiguous rois
    ascontiguous = let
        srois = sort(rois)
        res = [srois[1]]
        for roi in srois[2:end]
            if length(intersect(res[end], roi)) > 0
                res[end] = min(roi.start, res[end].start):max(roi.stop, res[end].stop)
            else
                push!(res, roi)
            end
        end
        res
    end
    return invert(ascontiguous, channel(minE, det):channel(maxE, det))
end

"""
    fitcontinuum(
      spec::Spectrum,
      det::EDSDetector,
      resp::AbstractArray{<:Real,2}; #
      minE::Float64 = 1.5e3,
      maxE::Float64 = 0.95 * spec[:BeamEnergy],
      brem::Type{<:NeXLBremsstrahlung} = Castellano2004a,
      mc::Type{<:MatrixCorrection} = Riveros1993,
    )

Fit the continuum from ROIs determined from the data within the spectrum (:Composition, :BeamEnergy & :TakeOffAngle).
The ROIs are computed using `continuumrois(...)` and each roi is fit seperately.

    fittedcontinuum(
      spec::Spectrum,
      det::EDSDetector,
      resp::AbstractArray{<:Real,2}; #
      mode = :Global [ | :Local ] # Fit to all ROIs simultaneously (:Global) or to each roi independently (:Local)
      minE::Float64 = 1.5e3,
      maxE::Float64 = 0.95 * spec[:BeamEnergy],
      brem::Type{<:NeXLBremsstrahlung} = Castellano2004a,
      mc::Type{<:MatrixCorrection} = Riveros1993,
    )::Spectrum

Fit the continuum under the characteristic peaks by fitting the closest continuum ROIs.  The low energy peaks are
fit using the continuum immediately higher in energy and the high energy peaks are fit using the continuum on both
sides.
"""
function fitcontinuum(
    spec::Spectrum,
    det::EDSDetector,
    resp::AbstractArray{<:Real,2}; #
    minE::Float64 = 1.5e3,
    maxE::Float64 = 0.95 * spec[:BeamEnergy],
    brem::Type{<:NeXLBremsstrahlung} = Castellano2004a,
    mc::Type{<:MatrixCorrection} = Riveros1993,
)::Vector{Float64}
    @assert haskey(spec, :Composition) "The fitcontinuum(...) function requires the spec[:Composition] property."
    @assert haskey(spec, :BeamEnergy) "The fitcontinuum(...) function requires the spec[:BeamEnergy] property."
    @assert haskey(spec, :TakeOffAngle) "The fitcontinuum(...) function requires the spec[:TakeOffAngle] property."
    cmod = ContinuumModel(
        spec[:Composition],
        spec[:BeamEnergy],
        spec[:TakeOffAngle],
        matrixcorrection = mc,
        bremsstrahlung = brem,
    )
    # meas is the raw continuum shape.  It needs to be scaled to the unknown.
    minEmod = max(50.0, energy(lld(spec), spec))
    model = resp * map(e -> e > minEmod ? emitted(cmod, e) : 0.0, energyscale(spec))
    prevroi, brem = missing, counts(spec, Float64)
    for roi in continuumrois(elms(spec, true), det, minE, maxE)
        currrois = ismissing(prevroi) ? (roi,) : (prevroi, roi)
        k =
            sum(sum(counts(spec, roi, Float64)) for roi in currrois) /
            sum(sum(model[roi]) for roi in currrois)
        peak = (ismissing(prevroi) ? lld(det) : prevroi.stop):roi.start
        brem[peak] = k * model[peak]
        prevroi = roi
    end
    return brem
end

function fittedcontinuum(
    spec::Spectrum,
    det::EDSDetector,
    resp::AbstractArray{<:Real,2}; #
    mode::Symbol = :Global,
    minE::Float64 = 1.5e3,
    maxE::Float64 = 0.95 * spec[:BeamEnergy],
    brem::Type{<:NeXLBremsstrahlung} = Castellano2004a,
    mc::Type{<:MatrixCorrection} = Riveros1993,
)::Spectrum
    function localfittedcontinuum(spec, det, resp, minE, maxE, brem, mc)::Spectrum
        res = Spectrum(
            spec.energy,
            fitcontinuum(spec, det, resp, minE = minE, maxE = maxE, brem = brem, mc = mc),
            copy(spec.properties),
        )
        res[:Name] = "Brem[Local][$(spec[:Name])]"
        return res
    end
    globalfittedcontinuum(spec, det, resp, minE, maxE, brem, mc)::Spectrum = fitcontinuum(
        spec,
        resp,
        continuumrois(elms(spec, true), det, minE, maxE),
        brem = brem,
        mc = mc,
    )
    @assert (mode == :Global) || (mode == :Local) "mode must equal :Global | :Local in fitted continuum"
    return mode == :Global ? globalfittedcontinuum(spec, det, resp, minE, maxE, brem, mc) : #
           localfittedcontinuum(spec, det, resp, minE, maxE, brem, mc)
end

"""
    subtractcontinuum(
      spec::Spectrum,
      det::EDSDetector,
      resp::AbstractArray{<:Real,2}; #
      minE::Float64 = 1.5e3,
      maxE::Float64 = 0.95 * spec[:BeamEnergy],
      brem::Type{<:NeXLBremsstrahlung} = Castellano2004a,
      mc::Type{<:MatrixCorrection} = Riveros1993,
    )::Spectrum

Computes the characteristic-only spectrum by subtracting the .
"""
function subtractcontinuum(
    spec::Spectrum,
    det::EDSDetector,
    resp::AbstractArray{<:Real,2}; #
    minE::Float64 = 1.5e3,
    maxE::Float64 = 0.95 * spec[:BeamEnergy],
    brem::Type{<:NeXLBremsstrahlung} = Castellano2004a,
    mc::Type{<:MatrixCorrection} = Riveros1993,
)::Spectrum
    res = Spectrum(
        spec.energy,
        counts(spec, Float64) -
        fitcontinuum(spec, det, resp, minE = minE, maxE = maxE, brem = brem, mc = mc),
        copy(spec.properties),
    )
    res[:Name] = "CharOnly[$(res[:Name])]"
    return res
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
################################################################################################
# OpenCL support for GPU-accelerated OpenCV
################################################################################################
#  See http://docs.opencv.org/ref/master/dc/d83/group__core__opencl.html

# Set up OpenCL
setUseOpenCL(select=true) = @cxx cv::ocl::setUseOpenCL(select)

# ocl::getOpenCLPlatforms
# Returns the list of OpenCL platforms
# int ocl::getOpenCLPlatforms(PlatformsInfo& platforms)
# cxx"""
# int getOCLPlatforms(void)
# {
#     cv::ocl::PlatformInfo& platforms;
#     return(getOpenCLPlatforms(platforms));
# }
# """
getOpenCLPlatforms() = @cxx getOpenCLPlatforms()

# ocl::getOpenCLDevices
# Returns the list of devices
# int ocl::getOpenCLDevices(DevicesInfo& devices, int deviceType=CVCL_DEVICE_TYPE_GPU,
#    const PlatformInfo* platform=NULL )
# Parameters:
# devices – Output variable
# deviceType – Bitmask of TYPE_GPU, TYPE_CPU or TYPE_DEFAULT
# deviceType
# TYPE_DEFAULT
# TYPE_CPU
# TYPE_GPU
# TYPE_ACCELERATOR
# TYPE_DGPU
# TYPE_IGPU
# platform – Specifies preferrable platform
# cxx"""
# int getOCLDevices(int deviceType)
# {
#     cv::cuda::DeviceInfo& devices;
#     const cv::ocl::PlatformInfo* platform=NULL;

#     return(getOpenCLDevices(devices, deviceType, platform));
# }
# """
getOpenCLDevices(C_NULL, deviceType=TYPE_GPU) =  @cxx getOpenCLDevices(C_NULL, deviceType)

# ocl::setDevice
# Initialize OpenCL computation context
# void ocl::setDevice(const DeviceInfo* info)
# cxx"""
# void setDevice(void)
# {
#     const cv::cuda::DeviceInfo* info;
#     cv::ocl::setDevice(info);
# }
# """
setDevice() = @cxx setDevice()

# ocl::initializeContext
# Alternative way to initialize OpenCL computation context
# This function can be used for context initialization with D3D/OpenGL interoperability.
# void ocl::initializeContext(void* pClPlatform, void* pClContext, void* pClDevice)
# Parameters:
# pClPlatform – selected platform_id (via pointer, parameter type is cl_platform_id*)
# pClContext – selected cl_context (via pointer, parameter type is cl_context*)
# pClDevice – selected cl_device_id (via pointer, parameter type is cl_device_id*)

initializeContext() = @cxx initializeContext()

# ocl::setBinaryPath
# void ocl::setBinaryPath(const char* path)
# Parameters:
# path – the path of OpenCL kernel binaries
setBinaryPath(path::String) = @cxx cv::ocl::setBinaryPath(pointer(path))

################################################################################################
# Data structures Universal Mat => "UMat"

# UMat::UMat()
UMat() = @cxx cv::UMat()

# UMat::UMat(int rows, int cols, int type)
UMat(rows::Int, cols::Int, matType::CV_MatType) = @cxx cv::UMat(rows, cols, matType)

# UMat::UMat(Size size, int type)
#cxx""" cv::UMat UMat(int rows, int cols, int matType){ cv::UMat img(cv::Size(rows,cols), matType); return(img); }"""
UMat(size, matType::CV_MatType) = @cxx cv::UMat(size, matType)

# UMat::UMat(int rows, int cols, int type, const Scalar& s)
#cxx""" cv::UMat UMat(int rows, int cols, int matType){ cv::UMat img(cv::Size(rows,cols), matType); return(img); }"""
UMat(rows::Int, cols::Int, matType::CV_MatType, s) = @cxx cv::UMat(rows, cols, matType, s)

# UMat::UMat(Size size, int type, const Scalar& s)
UMat(size, matType::CV_MatType, s) = @cxx cv::UMat(size, matType, s)

# UMat::UMat(const UMat& m)
UMat(m) = @cxx cv::UMat(m)

# UMat::UMat(int ndims, const int* sizes, int type)
# const psizes(sizes) = pointer([sizes::Int])
UMat(ndims::Int, psizes::Ptr{Int}, matType::CV_MatType) = @cxx cv::UMat(ndims, psizes, matType)

# UMat::UMat(int ndims, const int* sizes, int type, const Scalar& s)
UMat(ndims::Int, psizes::Ptr{Int}, matType::CV_MatType, s) = @cxx cv::UMat(ndims, psizes, matType, s)

# UMat::UMat(const UMat& m, const Rect& roi)
UMat(img, roi) = @cxx cv::UMat(img, roi)

# UMat::UMat(const UMat& m, const Range* ranges)
# const ranges = pointer(range)
UMat(img, ranges) = @cxx cv::UMat(img, ranges)

# UMat class methods
# addition
# cxx""" cv::UMat add(cv::UMat img1, cv::UMat img2) { return(img1 + img2); } """
# imadd(img1, img2) = @cxx add(img1, img2)

# substract
# cxx""" cv::UMat substract(cv::UMat img1, cv::UMat img2) { return(img1 - img2); } """
# imsubstract(img1, img2) = @cxx substract(img1, img2)

# UMat::row
cxx""" cv::UMat row(cv::UMat img, int y) { return(img.row(y)); } """
row(img, x::Int) = @cxx row(img, x)

# UMat::col
cxx""" cv::UMat col(cv::UMat img, int x) { return(img.col(x)); } """
col(img, y::Int) = @cxx col(img, y)

# UMat::rowRange
cxx""" cv::UMat rowRange(cv::UMat img, int startrow, int endrow) { return(img.rowRange(startrow, endrow)); } """
rowRange(img, startrow::Int, endrow::Int) = @cxx rowRange(img, startrow, endrow)

cxx""" cv::UMat rowRange(cv::UMat img, const cv::Range& r) { return(img.rowRange(r)); } """
rowRange(img, range) = @cxx rowRange(img, range)

# UMat::colRange
# const range = cvRange(start::Int, tend::Int)
cxx""" cv::UMat colRange(cv::UMat img, int startcol, int endcol) { return(img.colRange(startcol, endcol)); } """
colRange(img, startcol::Int, endcol::Int) = @cxx colRange(img, startcol, endcol)

cxx""" cv::UMat colRange(cv::UMat img, const cv::Range& r) { return(img.colRange(r)); } """
colRange(img, range) = @cxx colRange(img, range)

# UMat::diag
# d=0 is the main diagonal
# d>0 is a diagonal from the lower half
# d=1 below the main one
# d<0 is a diagonal from the upper half

cxx""" cv::UMat diag(cv::UMat img, int d=0) { return(img.diag(d)); } """
diag(img, d::Int) = @cxx diag(img, d)

cxx""" cv::UMat diag(cv::UMat img, const cv::UMat& m) { return(img.diag(m)); } """
diag(img, m) = @cxx diag(img, m)

# UMat::clone()
cxx""" cv::UMat clone(cv::UMat img) { return(img.clone()); } """
clone(img) = @cxx clone(img)

# UMat::copyTo
cxx""" void copy(cv::UMat &out, cv::UMat &img) { img.copyTo(out); } """
cxx""" void copyToGPU(cv::UMat &out, cv::Mat &img) { img.copyTo(out); } """
cxx""" void copyToCPU(cv::Mat &out, cv::UMat &img) { img.copyTo(out); } """
cxx""" void copyTomask(cv::UMat &img, cv::UMat &mask, cv::UMat &out) { img.copyTo(out, mask); } """
copy(out, img) = @cxx copy(out, img)
copyToGPU(out, img) = @cxx copyToGPU(out, img)
copyToCPU(out, img) = @cxx copyToCPU(out, img)
copyTomask(img, mask, out) = @cxx copyTomask(img, mask, out)

cxx""" cv::UMat imageROI(cv::UMat img, cv::Rect roi) { return img(roi); } """
ROImage(img, ROI) = @cxx imageROI(img, ROI)

# UMat::convertTo

# UMat::assignTo
cxx""" void  assignTo(cv::UMat img, cv::UMat& m, int type) { img.assignTo(m, type); } """
assignTo(img, m, rtype=-1) = @cxx assignTo(img, m, rtype)

# UMat::reshape
# rows = 0 (no change)
cxx""" cv::UMat reshape(cv::UMat img, int cn, int rows) { return(img.reshape(cn, rows)); } """
reshape(img, ch::Int, rows=0) = @cxx reshape(img, ch, rows)

# UMat::t() (transpose)
# cxx""" cv::UMat transpose(cv::UMat img, double lambda) { return(img.t()*lambda); } """
# transpose(img, lambda) = @cxx transpose(img, lambda)

# UMat::inv (invert)
cxx""" cv::UMat inv(cv::UMat img, int method) { return(img.inv(method)); } """
inv(img, method=DECOMP_LU) = @cxx inv(img, method)

# UMat::mul (mutiply)   # element-wise only
cxx""" cv::UMat mul(cv::UMat img, double scale) { return(img.mul(scale)); } """
mul(img, scale=1) = @cxx mul(img, scale)

# UMat::cross (cross-product of 2 Vec<float,3>)
# cxx""" cv::UMat cross(cv::UMat img, cv::UMat m) { return(img.cross(m)); } """
# cross(img, m) = @cxx cross(img, m)

# UMat::dot (Computes a dot-product of two equally sized matrices)
cxx"""  double dot(cv::UMat img, cv::UMat m) { return(img.dot(m)); } """
dot(img, m) = @cxx dot(img, m)

# UMat::zeros()
# ndims –> Array dimensionality.
# rows  –> Number of rows.
# cols  –> Number of columns.
# size  –> Alternative to the matrix size specification Size(cols, rows) .
# sz    –> Array of integers specifying the array shape.
# type  –> Created matrix type.

cxx""" cv::UMat zerosUM(int rows, int cols, int matType) {cv::UMat A = cv::UMat::zeros(rows, cols, matType); return (A); }"""
cxx""" cv::UMat zerosUM(cv::Size size, int matType) {cv::UMat A = cv::UMat::zeros(size, matType); return (A); }"""
cxx""" cv::UMat zerosUM(int ndims, int* sz, int matType) {cv::UMat A = cv::UMat::zeros(ndims, sz, matType); return (A); }"""
zerosUM(rows::Int, cols::Int, matType::CV_MatType) = @cxx zerosUM(rows, cols, matType)
zerosUM(size, matType::CV_MatType) = @cxx zerosUM(size, matType)
zerosUM(ndims::Int, sz::Ptr{Int32}, matType::CV_MatType) = @cxx zerosUM(ndims, sz, matType)

# UMat::ones()
onesUM(rows::Int, cols::Int, matType::CV_MatType) = @cxx cv::UMat::ones(rows, cols, matType)
onesUM(size, matType::CV_MatType) = @cxx cv::UMat::ones(size, matType)
onesUM(ndims::Int, sz::Ptr{Int32}, matType::CV_MatType) = @cxx cv::UMat::ones(ndims, sz, matType)

# UMat::eye()
eyeUM(rows::Int, cols::Int, matType::CV_MatType) = @cxx cv::UMat::eye(rows, cols, matType)
eyeUM(size, matType::CV_MatType) = @cxx cv::UMat::eye(size, matType)

cxx""" cv::Size size(cv::UMat img) {  return(img.size()); } """
sizeU(img) = @cxx size(img)                                   # returns Size(cols, rows), if matrix > 2d, size = (-1,-1)

cxx""" int channels(cv::UMat img) { return(img.channels()); } """
channelsU(img) = @cxx channels(img)                           # number of matrix channels

cxx""" int cvtype(cv::UMat img) { return(img.type()); } """
cvtypevalU(img) = @cxx cvtype(img)


text:
module SeisBeamforming

using SeisIO, SeisNoise, JLD2, FileIO, PlotlyJS, Geodesy, Statistics, LinearAlgebra, Dates

export seisnoisebeamform

#include("./src/Types/BeamData.jl")
#include("./src/Types/show.jl")

include("./Types/BeamData.jl")
include("./Types/show.jl")

function seisnoisebeamform(InputDict::Dict)

    finame 				= InputDict["finame"]
	fopath        		= InputDict["fopath"]
	IsPlotFigure = InputDict["IsPlotFigure"]

    stationlist 		= load(finame, "info/stationlist")
	DLtimestamplist 	= load(finame, "info/DLtimestamplist")
	DL_time_unit		= load(finame, "info/DL_time_unit")
	starttime 			= load(finame, "info/starttime")
	endtime 			= load(finame, "info/endtime")

	channel			= InputDict["channel"]
	bp_len 			= InputDict["bp_len"] #hourly beamforming
	bp_step 		= InputDict["bp_step"] # 30 minutes overlap
	freq_min 		= InputDict["freq_min"]
	freq_max 		= InputDict["freq_max"]
	sampling_freq 	= InputDict["sampling_freq"]
	freq_range 		= InputDict["freq_range"]
	fexamin_num 	= InputDict["fexamin_num"] # num of examined freqency
	vel_range 		= InputDict["vel_range"]
	Δvel   			= InputDict["Δvel"]
	Δθ    			= InputDict["Δθ"]

    # store station location
    lat = Float64[]; lon = Float64[]; x = Float64[]; y = Float64[];

    for i = 1:length(stationlist)
        Stemp = load(finame, joinpath(DLtimestamplist[1], stationlist[i]))
        if any(occursin.(channel, stationlist[i][end-2:end]))
            push!(lat, Stemp.loc.lat)
            push!(lon, Stemp.loc.lon)
            #- ignore elevation -#
        end
    end

	if isempty(lat); error("no available station to compute beampower in this array. abort."); end

	# project lat and lon to x-y coordinate
    origin_lla = LLA(mean(lat), mean(lon))
    trans = ENUfromLLA(origin_lla, wgs84)

    for i = 1:length(stationlist)
        point_lla = LLA(lat[i], lon[i])
        point_enu = trans(point_lla)
        push!(x, point_enu.e)
        push!(y, point_enu.n)
    end

    # check with plotting
	if IsPlotFigure
	    trace1 = PlotlyJS.scatter(;x=x./1e3, y=y./1e3,
	                          mode="markers", name="Stations",
							  textposition="top center",
	                          text=stationlist,
	                          marker_size=12)
	    trace2 = PlotlyJS.scatter(;x=[0.0], y=[0.0],
							mode="markers",
	 						marker=attr(symbol=:star),
							name="Origin",
							marker_size=12)
	    layout = Layout(;width=500, height=500,
					 showlegend=true,
	                 #xaxis_range=[-15, 15],
	                 #yaxis_range=[-15, 15],
	                 xaxis=attr(title="x [km]"),
					 yaxis=attr(title="y [km]"))
	    p = PlotlyJS.plot([trace1, trace2], layout)
		display(p)
	end

	# save info into jld2
	jldopen(fopath, "w") do file
		file["info/DLtimestamplist"]  = DLtimestamplist;
		file["info/stationlist"]  = stationlist;
		file["info/starttime"]   = starttime;
		file["info/endtime"]     = endtime;
	end

	# Loop with DL time sample; this becomes unit duration of bpstack
	for i = 1:length(DLtimestamplist)
		println(DLtimestamplist[i])

		# make beampower timestamp: bptimestamp
		t1 = parse.(Int64, split(DLtimestamplist[i], ".")[1:2])
		hm = split(DLtimestamplist[1], ".")[3]
		m1, d1 = SeisIO.j2md(t1[1], t1[2])
		st = DateTime(join([t1[1], "-", m1, "-",  d1, hm]))

		bptimestamp = []
		ttemp = deepcopy(st)

		while ttemp <= st + Second(DL_time_unit) - Second(bp_len)
			push!(bptimestamp, ttemp)
			ttemp += Second(bp_step)
		end

		# compute fft for all stations
		FFT1 = Array{FFTData, 1}(undef, length(stationlist))

		for j = 1:length(stationlist)
			Stemp = load(finame, joinpath(DLtimestamplist[i], stationlist[j]))
			Stemp.t[1,2]  = round(Stemp.t[1,2], sigdigits=13) #to avoid t_len skip
			try
				FFT1[j] = compute_fft(SeisData(Stemp), freq_min, freq_max,
			 	sampling_freq, Int(bp_step), Int(bp_len))
			catch # unknow error with FFT1

				FFT1[j] = FFTData(t = [0.0])
				println("skip station $(stationlist[j])")
			end
		end

		# assume all station have same sampling frequency
		L = sampling_freq * bp_len
		fvec0 = sampling_freq * collect(0:L/2)/L

		# compute examined frequencies
		fexam = exp10.(range(log10(freq_range[1]),stop=log10(freq_range[2]),length=fexamin_num))
		exfreq = []
		#@doc "exfreq contains (id, frequency) in the FFT.fft array to compute Y array." exfreq
		for j = 1:length(fexam)
			id1 = findfirst(x -> x>=fexam[j], fvec0)
			push!(exfreq, (Int(id1), fvec0[id1]))
		end

		# compute examined phase velocity and θ
		exvel = collect(vel_range[1]:Δvel:vel_range[2])
		exθ = deg2rad.(collect(0:Δθ:360))

		# assign size of BeamData.bp
		size_v = length(exvel)
		size_θ = length(exθ)
		size_t = length(bptimestamp)

		Beam = Array{BeamData, 1}(undef, length(exfreq))

		for l = 1:length(exfreq) #loop every examined frequency

			Btemp = BeamData()
			Btemp.name = join(channel, ",")
			Btemp.bp_len = bp_len
			Btemp.bp_step = bp_step
			Btemp.f = [exfreq[l][2]]
			Btemp.c = exvel
			Btemp.θ = exθ
			Btemp.t = d2u.(bptimestamp)
			Btemp.bp = Array{Complex{Float32},3}(undef, size_θ, size_v, size_t)
			Btemp.tstack = d2u(bptimestamp[round(Int64,length(bptimestamp)/2)+1]) #mean time of stack duration
			Btemp.tstack_len = parse(Int64, DL_time_unit)
			Btemp.bpstack = zeros(ComplexF32, size_θ, size_v)

			stack_count = 0
			for k = 1:length(bptimestamp) #loop every bp_len

				# compute size of e and Y at this loop: i.e. gather available stations at this timewindow
				Y_station_and_timestampid = []
				#@doc "Y_station_and_timestampid contains the station id and timestamp id corresponding to bptimestamp[k]
				# to compute Y array." Y_station_and_timestampid

				for stid = 1:length(stationlist)
					for timestampid = 1:length(FFT1[stid].t)
						if u2d(FFT1[stid].t[timestampid]) == bptimestamp[k]
							#this station is used as one station in seismic array
							push!(Y_station_and_timestampid, (stid, timestampid))
						end
					end
				end

				if length(Y_station_and_timestampid) >= 2 # if stations are more than 2 in array

					# loop below is to compute one polar coord figure
					stack_count += 1

					for eθid in 1:length(exθ)
						eθ = exθ[eθid]

						for evid = 1:length(exvel) #loop examined phase velocity
							ev = exvel[evid]

							# examined frequency
							ω_id = exfreq[l][1]
							ω	 = 2 * pi * exfreq[l][2] # ω = angular frequency
							kv   = (ω / ev) .* [cos(eθ), sin(eθ)]

							#compute the array plane wave response [exp(-ikv.xN)]^*T
							xv = Float32[]
							ee = ComplexF32[]

							for stid = 1:length(Y_station_and_timestampid)
								stationid = Y_station_and_timestampid[stid][1]
								xv= [x[stationid]; y[stationid]]
								push!(ee, exp(-im*dot(kv, xv)))
							end

							# compute spatial correlation matrix YY
							Y = ComplexF32[]

							for stid = 1:length(Y_station_and_timestampid)
								stationid = Y_station_and_timestampid[stid][1]
								timestampid = Y_station_and_timestampid[stid][2]
								f1 = FFT1[stationid].fft[ω_id, timestampid]
								push!(Y, f1)
							end

							# compute e'YY'e
							YY = Y*transpose(conj.(Y))

							Btemp.bp[eθid, evid, k] = transpose(conj.(ee)) * YY * ee

							#add into stack
							Btemp.bpstack[eθid, evid] += transpose(conj.(ee)) * YY * ee
						end
					end

				else
					println("$(bptimestamp[k]): num of station in array is less than two; thus skip this beam forming.")
					Btemp.bp[:, :, k] = zeros(ComplexF32, size_θ, size_v)
				end
			end

			# normalize stacked beam power by number of stacking
			Btemp.bpstack = Btemp.bpstack ./ stack_count

			# store into beam data struct
			Beam[l] = Btemp
		end

		# save info into jld2
		jldopen(fopath, "r+") do file
			varname = joinpath("BeamData", DLtimestamplist[i])
			file[varname]  = Beam;
		end
	end

	print("Beam forming is successfully done at: ")
	println(now())

end

end # module

max_stars_repo_name: kura-okubo/SeisBeamforming

max_issues_repo_name: kura-okubo/SeisBeamforming

max_forks_repo_name: kura-okubo/SeisBeamforming


text:
using Distances
using LinearAlgebra
using Distributions
using Mmap

width = 5
height = 10

function write_mmap(width, height)
    io = open("/tmp/mmap.bin", "w+")
    compressed_pixels = Mmap.mmap(io, Matrix{Float32}, (width, height))
    compressed_pixels = rand(Float32, width, height)
    println("write_mmap: ", compressed_pixels)
    Mmap.sync!(compressed_pixels)
    close(io)
end

function write_mmap2(width, height)
    io = open("/tmp/mmap.bin", "w+")
    compressed_pixels = rand(Float32, width, height)
    println("write_mmap: ", compressed_pixels)
    write(io, compressed_pixels)
    close(io)
end

function read_mmap(width, height)
    io = open("/tmp/mmap.bin")
    compressed_pixels = Mmap.mmap(io, Matrix{Float32}, (width, height))
    close(io)

    println("read_mmap: ", compressed_pixels)
end

write_mmap2(width, height)
read_mmap(width, height)

function testmmap()
    data = Float32.(rand(10000, 15))
    Eucldist = pairwise(Euclidean(), data, dims = 1)
    D = maximum(Eucldist .^ 2)
    sigma2hat = mean(((Eucldist .^ 2)./D)[tril!(trues(size((Eucldist .^ 2) ./ D)), -1)])
    L = exp.(-(Eucldist .^ 2 / D) / (2 * sigma2hat))
    s = open("/tmp/mmap.bin", "w+")
    write(s, size(L, 1))
    write(s, size(L, 2))
    write(s, L)
    close(s)

    # deref and gc collect
    Eucldist = data = L = zeros(Float32, 2, 2)
    GC.gc()

    s = open("/tmp/mmap.bin", "r+") # allow read and write
    m = read(s, Int)
    n = read(s, Int)
    L = Mmap.mmap(s, Matrix{Float32}, (m, n))  # now L references the file contents
    K = eigen(L)
    println(K)
end

# testmmap()
# @time testmmap()  # 109.657995 seconds (17.48 k allocations: 4.673 GiB, 0.73% gc time)


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, "..", ".."))
using ZigZagBoomerang
Pkg.activate(@__DIR__)
cd(@__DIR__)
#using Revise
using DataStructures
using LinearAlgebra
using Random
using SparseArrays
using Test
using FileIO
using Statistics
using Makie, AbstractPlotting
using ProgressMeter

function gridlaplacian(T, m, n)
    S = sparse(T(0.0)I, n*m, n*m)
    linear = LinearIndices((1:m, 1:n))
    for i in 1:m
        for j in 1:n
            for (i2, j2) in ((i + 1, j), (i, j + 1))
                if i2 <= m && j2 <= n
                    S[linear[i, j], linear[i2, j2]] -= 1.
                    S[linear[i2, j2], linear[i, j]] -= 1.

                    S[linear[i, j], linear[i, j]] += 1.
                    S[linear[i2, j2], linear[i2, j2]] += 1.
                end
            end
        end
    end
    S
end

Random.seed!(1)

# Define precision operator of a Gaussian random field (sparse matrix operating on `vec`s of `n*n` matrices)
#n = 100
n = 300
const σ2 = 0.5
Γ0 = 10*gridlaplacian(Float64, n, n)^2
Γ = 0.1I + Γ0
mat(x) = reshape(x, (n, n)) # vector to matrix
function mat0(y)
    mat(y  .- 0.1)
end
# Γ is very sparse
@show nnz(Γ)/length(Γ) # 0.000496

# Corresponding Gaussian potential
# ϕ(x, Γ, y) = 0.5*x'*Γ*x  + dot(x - y, x - y)/(2*σ2) # not used by the program

# Define ∇ϕ(x, i, Γ) giving the partial derivative of ϕ(x) with respect to x[i]
∇ϕ(x, i, Γ, y) = ZigZagBoomerang.idot(Γ, i, x)  + (x[i]-y[i])/σ2 # more efficient that dot(Γ[:, i], x)

# Random initial values
t0 = 0.0
sz = 1.0
h(x, y) = x^2+(5y/4-sqrt(abs(x)))^2
heart = [ max.(1 - h(x, y), 0) for x in range(-1.5-sz,1.5+sz, length=n), y   in range(-1.1-sz,1.9+sz, length=n)]
image(heart)
μ0 = 5.0*vec(heart)
y = μ = μ0 + randn(n*n)
μpost = yhat = (I/σ2 + Γ)\y
Γpost = (Γ + I/σ2)

x0 = μpost
θ0 = rand([-1.0,1.0], n*n)

# Rejection bounds
c = 0.05*[norm(Γpost[:, i], 2) for i in 1:n*n]

# Define ZigZag
Z = ZigZag(Γpost, μpost)
# or try the FactBoomerang
#Z = FactBoomerang(Γ, x0*0, 0.1)

κ = 2.0*ones(length(x0))
# Run sparse ZigZag for T time units and collect trajectory
T = 200.0
#@time trace0, (tT, xT, θT), (acc, num) = spdmp(∇ϕ, t0, x0, θ0, T, c, Z, Γ, μ; structured = true, adapt = false)

prg = Progress(20, 1)
su = false
adapt = false
trace, (t, x, θ), (acc, num), c = @time sspdmp(∇ϕ, t0, x0, θ0, T, c, Z, κ, Γ, μ; structured = true,
                                                strong_upperbounds = su ,
                                                adapt = adapt,
                                                progress = ()->next!(prg))

function fastmean(trace::ZigZagBoomerang.Trace, T)
    x = copy(trace.x0)
    θ = copy(trace.θ0)
    y = 0*x
    t = fill(trace.t0, length(x0))
    k = 1
    scale = 1/(2T)
    while k < length(trace.events)
        t2, i, xi, θi = trace.events[k]
        k += 1
        y[i] += (x[i]+xi)*(t2-t[i])*scale
        t[i] = t2
        x[i] = xi
        θ[i] = θi
    end
    y
end

@time est = fastmean(trace, T)
#display(scene)

#scene1, _ = heatmap(abs.([mat(mean(last.(traj0[end÷2:end]))) mat(est) ]))

scene2b, _ = heatmap(abs.(mat(μ0 - est)))

#error("stop")
if false

# Prepare surface plot
    M = Node(mat0(traj[end].second))
    scene = Scene()
    surface!(scene, M; shading=false, show_axis=false, colormap = :oleron, colorrange = (-3,3))
    Makie.scale!(scene, 1.0, 1.0, 2.0)
    zlims!(scene, -5.0, 5.0)
    # Movie frames
    mkpath(joinpath(@__DIR__, "output"))
    for i in eachindex(traj)
        M[] = mat0(traj[i].second)
        FileIO.save(joinpath(@__DIR__, "output", "surfs$i.png"), scene)
    end

    # Make video
    dir = joinpath(@__DIR__, "output")
    run(`ffmpeg -y -r 40 -f image2 -i $dir/surfs%d.png -vcodec libx264 -crf 25  -pix_fmt yuv420p $(typeof(Z).name)field.mp4`)
end

# Save hearts
save("hearttruec.png", image(mat(μ0), interpolate=false))
save("heartc.png", image(mat(y), interpolate=false))
save("hearthatc.png", image(mat(yhat), interpolate=false))
save("heartpostmeanc.png", image(mat(est), interpolate=false))
save("hearterrorc.png", scene2b)

mean((mat(μ0 - est)).^2)

@show mean(abs.(mat(μ0 - est)))
@show extrema(μ0)

@show  extrema(y)

@show extrema(yhat)

@show extrema((mat(est)))

@show extrema(abs.(mat(μ0 - est)))

using ZigZagBoomerang: FactTrace, Trace, sep

function move_forward1(t, x, θ, t2, ::ZigZag)
    t2, x + θ*(t2-t), θ
end

struct Marginal{T}
    FT::T
    i::Int
end
Base.IteratorSize(::Marginal) = Iterators.SizeUnknown()

function Base.iterate(D::Marginal{<:FactTrace})
    FT = D.FT
    t, x, θ = FT.t0, FT.x0[D.i],  FT.θ0[D.i]
    t => x, (t, x, θ, 1)
end

function Base.iterate(D::Marginal{<:FactTrace{T}}, (t, x, θ, k)) where {T<:Union{FactBoomerang,ZigZag}}
    FT = D.FT
    k > length(FT.events) && return nothing
    local t2, i, xi, θi
    while true
        k == length(FT.events) && return t => x, (t, x, θ, k + 1)
        t2, i, xi, θi = FT.events[k]
        k += 1
        i == D.i && break
    end
    t, x, θ = move_forward1(t, x, θ, t2, FT.F)
    t = t2
    x = xi
    θ = θi
    return t => x, (t, x, θ, k + 1)
end

scene3 = lines(sep(Marginal(trace, 44543))..., color=:red)
lines!(sep(Marginal(trace, 44544))..., color=:blue)
lines!(sep(Marginal(trace, 44524))..., color=:black)
lines!(sep(Marginal(trace, 44644))..., color=:orange)

save("hearttracec.png", scene3)


text:
const MISSING_INITIAL_CONDITIONS_TIME_COUNT = 999.0

"""
Status Initis is always calculated based on the Power Output of the device
This is to make it easier to calculate when the previous model doesn't
contain binaries. For instance, looking back on an ED model to find the
IC of the UC model
"""
function status_init(canonical::Canonical,
                     devices::IS.FlattenIteratorWrapper{PSD}) where {PSD<:PSY.ThermalGen}

    key = ICKey(DeviceStatus, PSD)
    parameters = model_has_parameters(canonical)
    length_devices = length(devices)
    ini_conds = get_initial_conditions(canonical, key)
    # Improve here
    ref_key = parameters ? Symbol("P_$(PSD)") : :activepower

    if isempty(ini_conds)
        @info("Setting $(key.quantity) initial conditions for the status of all devices $(PSD) based on system data")
        ini_conds = canonical.initial_conditions[key] = Vector{InitialCondition}(undef, length_devices)
        for (ix, g) in enumerate(devices)
            status_value = (PSY.get_activepower(g) > 0) ? 1.0 : 0.0
            ini_conds[ix] = InitialCondition(canonical,
                                            g,
                                            ref_key,
                                            status_value)
        end
    else
        ic_devices = (ic.device for ic in ini_conds)
        for g in devices
            g in ic_devices && continue
            @info("Setting $(key.quantity) initial conditions for the status device $(g.name) based on system data")
                status_value = (PSY.get_activepower(g) > 0) ? 1.0 : 0.0
                push!(ini_conds, InitialCondition(canonical,
                                                  g,
                                                  ref_key,
                                                  status_value))
        end
    end

    @assert length(ini_conds) == length(devices)

    return

end

function output_init(canonical::Canonical,
                    devices::IS.FlattenIteratorWrapper{PSD}) where {PSD<:PSY.ThermalGen}

    key = ICKey(DevicePower, PSD)
    parameters = model_has_parameters(canonical)
    length_devices = length(devices)
    ini_conds = get_initial_conditions(canonical, key)
    # Improve this
    ref_key = parameters ? Symbol("P_$(PSD)") : :activepower

    if isempty(ini_conds)
        @info("Setting $(key.quantity) initial_condition of all devices $(PSD) based on system data")
        ini_conds = canonical.initial_conditions[key] = Vector{InitialCondition}(undef, length_devices)
        for (ix, g) in enumerate(devices)
                ini_conds[ix] = InitialCondition(canonical,
                                                g,
                                                ref_key,
                                                PSY.get_activepower(g),
                                                TimeStatusChange)
        end

    else
        ic_devices = (ic.device for ic in ini_conds)
        for g in devices
            g in ic_devices && continue
            @info("Setting $(key.quantity) initial_condition of device $(g.name) based on system data")
                push!(ini_conds, InitialCondition(canonical,
                                                g,
                                                ref_key,
                                                PSY.get_activepower(g),
                                                TimeStatusChange))
        end

    end

    @assert length(ini_conds) == length(devices)

    return

end

function duration_init(canonical::Canonical,
                        devices::IS.FlattenIteratorWrapper{PSD}) where {PSD<:PSY.ThermalGen}

    keys = [ICKey(TimeDurationON, PSD), ICKey(TimeDurationOFF, PSD)]
    parameters = model_has_parameters(canonical)
    length_devices = length(devices)
    ref_key = parameters ? Symbol("P_$(PSD)") : :activepower

    for (ik, key) in enumerate(keys)
        ini_conds = get_initial_conditions(canonical, key)

        if isempty(ini_conds)
            @info("Setting $(key.quantity) initial_condition of all devices $(PSD) based on system data")
            ini_conds = canonical.initial_conditions[key] = Vector{InitialCondition}(undef, length_devices)

            for (ix, g) in enumerate(devices)
                time_on = (PSY.get_activepower(g) > 0) ? MISSING_INITIAL_CONDITIONS_TIME_COUNT : 0.0
                time_off = PSY.get_activepower(g) <= 0 ? MISSING_INITIAL_CONDITIONS_TIME_COUNT : 0.0
                times = [time_on, time_off]
                ini_conds[ix] = InitialCondition(canonical,
                                                    g,
                                                    ref_key,
                                                    times[ik],
                                                    TimeStatusChange)
            end

        else

            ic_devices = (ic.device for ic in ini_conds if !isnothing(ic.cache))

            for g in devices
                g in ic_devices && continue
                time_on = (PSY.get_activepower(g) > 0) ? MISSING_INITIAL_CONDITIONS_TIME_COUNT : 0.0
                time_off = PSY.get_activepower(g) <= 0 ? MISSING_INITIAL_CONDITIONS_TIME_COUNT : 0.0
                times = [time_on, time_off]
                @info("Setting $(key.quantity) initial_condition of device $(g.name) based on system data")
                push!(ini_conds, InitialCondition(canonical,
                                                  g,
                                                  ref_key,
                                                  times[ik],
                                                  TimeStatusChange))
            end

        end

        @assert length(ini_conds) == length(devices)

    end

    return

end

######################### Initialize Functions for Storage #################################

function storage_energy_init(canonical::Canonical,
                             devices::IS.FlattenIteratorWrapper{PSD}) where {PSD<:PSY.Storage}

    key = ICKey(DeviceEnergy, PSD)
    parameters = model_has_parameters(canonical)
    length_devices = length(devices)
    ini_conds = get_initial_conditions(canonical, key)
    ref_key = parameters ? Symbol("E_$(PSD)") : :energy

    if isempty(ini_conds)
        @info("Setting $(key.quantity) initial_condition of all devices $(PSD) based on system data")
        ini_conds = canonical.initial_conditions[key] = Vector{InitialCondition}(undef, length_devices)

        for (ix, g) in enumerate(devices)
            ini_conds[ix] = InitialCondition(canonical,
                                                g,
                                                ref_key,
                                                PSY.get_energy(g))
        end
    else
        ic_devices = (ic.device for ic in ini_conds)
        for g in devices
            g in ic_devices && continue
            @info("Setting $(key.quantity) initial_condition of device $(g.name) based on system data")
            push!(ini_conds, InitialCondition(canonical,
                                            g,
                                            ref_key,
                                            PSY.get_energy(g)))
        end
    end

        @assert length(ini_conds) == length(devices)

        return

end


text:
print("Testing ACQR_simula: ")
nt = 1000
seed = 99
m = ACQR_simula(nt,seed)
y = m["y"]

if sum(y.^2) > 1e-6
	println("OK")
end


text:
"""
```
dlnMobility(
    dlnParams::DislocationParameters{T1,T2,T3,T4} where {T1 <: mobBCC,T2,T3,T4},
    matParams::MaterialParameters,
    network::DislocationNetwork,
    idx = nothing,
)
```
Compute the nodal force and velocities for a `BCC` material.

Original by Bruce Bromage at the Department of Materials of the University of Oxford, @brucebromage on github.

This is outdated, new capabilities include rotating the frame of reference and better handling of cross-slip.
"""
function dlnMobility(
    dlnParams::DislocationParameters{T1, T2, T3, T4} where {T1 <: mobBCC, T2, T3, T4},
    matParams::MaterialParameters,
    network::DislocationNetwork,
    idx = nothing,
)
    # Peierls-Nabarro stress for the bcc material.
    σPN = matParams.σPN
    # Drag coefficients.
    edgeDrag = dlnParams.dragCoeffs[:edge]
    screwDrag = dlnParams.dragCoeffs[:screw]
    climbDrag = dlnParams.dragCoeffs[:climb]
    lineDrag = dlnParams.dragCoeffs[:line]

    links = network.links
    bVec = network.bVec
    coord = network.coord
    maxConnect = network.maxConnect
    segForce = network.segForce
    connectivity = network.connectivity
    elemT = eltype(network.bVec)
    I3 = SMatrix{3, 3, elemT}(I)

    # Do it for all nodes if no list is provided.
    if isnothing(idx)
        numNode = network.numNode[1]
        idx = 1:numNode
    else
        numNode = length(idx)
    end

    nodeForce = zeros(elemT, 3, numNode)
    nodeVel = zeros(elemT, 3, numNode)

    # Loop through nodes.
    for (i, node1) in enumerate(idx)
        totalDrag = zeros(SMatrix{3, 3, elemT})
        iNodeForce = zeros(SVector{3, elemT})
        iNodeVel = zeros(SVector{3, elemT})
        numConnect = connectivity[1, node1]

        # Loop through number of connections.
        for connect in 1:numConnect
            link = connectivity[2 * connect, node1]
            colLink = connectivity[2 * connect + 1, node1]
            colOppLink = 3 - colLink
            node2 = links[colOppLink, link]

            # Line direction, continue to next iteration if norm is 0 and normalise line direction.
            t = SVector{3, elemT}(
                coord[1, node2] - coord[1, node1],
                coord[2, node2] - coord[2, node1],
                coord[3, node2] - coord[3, node1],
            )
            nT = norm(t)
            # If the segment length is zero, skip to the next iteration.
            nT < eps(elemT) && continue
            tN = 1 / nT
            t *= tN # Normalise t.

            # Segment force relevant to node1.
            iNodeConForce = SVector{3, elemT}(
                segForce[1, colLink, link],
                segForce[2, colLink, link],
                segForce[3, colLink, link],
            )
            # If the sress is lower than the Peierls-Nabarro stress, the node behaves as if the force acting on it is zero.
            norm(iNodeConForce) * tN < σPN ? iNodeConForce = zeros(SVector{3, elemT}) :
            nothing
            # Add force from this connection to the total force on the node.
            iNodeForce += iNodeConForce

            # Burgers Vector
            b = SVector{3, elemT}(bVec[1, link], bVec[2, link], bVec[3, link])
            nB = norm(b)
            # Burgers vector families of screw dislocations in BCC, in case there have been collisions and the Burgers vector norm is greater than 1.
            absB = abs.(b)
            absB /= minimum(absB)
            sqrtBType = sqrt(eps(elemT))
            bTypeArr1 = @. abs(1 - absB) < sqrtBType
            bTypeArr2 = @. abs(2 - absB) < sqrtBType
            bTypeArr3 = @. abs(3 - absB) < sqrtBType
            bType1 = sum(bTypeArr1)
            bType2 = sum(bTypeArr2)
            bType3 = sum(bTypeArr3)
            screw = 1 - nB < sqrtBType
            # Normalise Burgers vector.
            b /= nB

            nBnT_2 = nB * nT / 2
            # If it's not a screw-ish dislocation, do edge calculation and continue to next iteration.
            if !(
                bType1 == 3 && screw ||
                bType1 == 2 && bType2 == 1 && screw ||
                bType1 == 1 && bType2 == 1 && bType3 == 1 && screw
            )
                # ξ += |b|| * ||t||/2 * (ξ_climb * I + (ξ_line - ξ_climb) * t ⊗ t)
                totalDrag += nBnT_2 * (climbDrag * I3 + (lineDrag - climbDrag) * t ⊗ t)
                continue
            end

            # cos(θ) = (a ⋅ b)/(||a|| ||b||)
            # Both vectors have been normalised already.
            tDb = t ⋅ b
            cosSqθ = tDb^2
            sinSqθ = 1 - cosSqθ

            # ξ += ||b|| * ||t||/2 (ξ_screw * I + (ξ_line - ξ_screw) * t ⊗ t)
            totalDrag += nBnT_2 * (screwDrag * I3 + (lineDrag - screwDrag) * t ⊗ t)

            # If dislocation segment is pure screw skip to next iteration.
            sinSqθ < eps(typeof(tDb)) && continue

            # P, Q vectors as new local axes.
            p = b × t / sqrt(sinSqθ)
            q = p × t

            # Eqn (112) from Arsenlis et al 2007 MSMSE 15 553
            screwDragSq = screwDrag^2
            # screw_ξ_glide = 1/sqrt(sin^2(θ) / ξ_edge^2 + cos^2(θ) / ξ_screw^2)
            ScrewGlideDrag = 1 / sqrt(sinSqθ / edgeDrag^2 + cosSqθ / screwDragSq)
            # screw_ξ_climb = sqrt(sin^2(θ) * ξ_climb^2 + cos^2(θ) * ξ_screw^2)
            ScrewClimbDrag = sqrt(sinSqθ * climbDrag^2 + cosSqθ * screwDragSq)

            # ξ += ||b|| * ||t||/2 *
            #      ((screw_ξ_glide - ξ_screw) * q ⊗ q + (screw_ξ_climb - ξ_screw))
            totalDrag +=
                nBnT_2 * (
                    (ScrewGlideDrag - screwDrag) * q ⊗ q +
                    (ScrewClimbDrag - screwDrag) * p ⊗ p
                )
        end
        # Solve for velocity.
        # ξ v = f
        iNodeVel = totalDrag \ iNodeForce
        #=
        try
            iNodeVel = totalDrag \ iNodeForce
        catch SingularSystem
            origTotalDrag = totalDrag
            while true
                try
                    totalDrag += I3 * maximum(abs.(origTotalDrag)) * sqrt(eps(dType))
                    iNodeVel = totalDrag \ iNodeForce
        break
        catch SingularSystem
        end
            end
        end =#
        nodeForce[:, i] = iNodeForce
        nodeVel[:, i] = iNodeVel
    end

    return nodeForce, nodeVel
end
"""
```
dlnMobility!(
    dlnParams::DislocationParameters{T1,T2,T3,T4} where {T1 <: mobBCC,T2,T3,T4},
    matParams::MaterialParameters,
    network::DislocationNetwork,
    idx = nothing,
)
```
In-place computation of the nodal force and velocities for a `BCC` material.

Original by Bruce Bromage at the Department of Materials of the University of Oxford, @brucebromage on github.

This is outdated, new capabilities include rotating the frame of reference and better handling of cross-slip.
"""
function dlnMobility!(
    dlnParams::DislocationParameters{T1, T2, T3, T4} where {T1 <: mobBCC, T2, T3, T4},
    matParams::MaterialParameters,
    network::DislocationNetwork,
    idx = nothing,
)
    # Peierls-Nabarro stress for the bcc material.
    σPN = matParams.σPN
    # Drag coefficients.
    edgeDrag = dlnParams.dragCoeffs[:edge]
    screwDrag = dlnParams.dragCoeffs[:screw]
    climbDrag = dlnParams.dragCoeffs[:climb]
    lineDrag = dlnParams.dragCoeffs[:line]

    links = network.links
    bVec = network.bVec
    coord = network.coord
    nodeVel = network.nodeVel
    nodeForce = network.nodeForce
    maxConnect = network.maxConnect
    segForce = network.segForce
    connectivity = network.connectivity
    elemT = eltype(network.bVec)
    I3 = SMatrix{3, 3, elemT}(I)

    # Do it for all nodes if no list is provided.
    if isnothing(idx)
        numNode = network.numNode[1]
        idx = 1:numNode
    else
        numNode = length(idx)
    end

    # Loop through nodes.
    for node1 in idx
        totalDrag = zeros(SMatrix{3, 3, elemT})
        iNodeForce = zeros(SVector{3, elemT})
        iNodeVel = zeros(SVector{3, elemT})
        numConnect = connectivity[1, node1]

        # Loop through number of connections.
        for connect in 1:numConnect
            link = connectivity[2 * connect, node1]
            colLink = connectivity[2 * connect + 1, node1]
            colOppLink = 3 - colLink
            node2 = links[colOppLink, link]

            # Line direction, continue to next iteration if norm is 0 and normalise line direction.
            t = SVector{3, elemT}(
                coord[1, node2] - coord[1, node1],
                coord[2, node2] - coord[2, node1],
                coord[3, node2] - coord[3, node1],
            )
            nT = norm(t)
            # If the segment length is zero, skip to the next iteration.
            nT < eps(elemT) && continue
            tN = 1 / nT
            t *= tN # Normalise t.

            # Segment force relevant to node1.
            iNodeConForce = SVector{3, elemT}(
                segForce[1, colLink, link],
                segForce[2, colLink, link],
                segForce[3, colLink, link],
            )
            # If the sress is lower than the Peierls-Nabarro stress, the node behaves as if the force acting on it is zero.
            norm(iNodeConForce) * tN < σPN ? iNodeConForce = zeros(SVector{3, elemT}) :
            nothing
            # Add force from this connection to the total force on the node.
            iNodeForce += iNodeConForce

            # Burgers Vector
            b = SVector{3, elemT}(bVec[1, link], bVec[2, link], bVec[3, link])
            nB = norm(b)
            # Burgers vector families of screw dislocations in BCC, in case there have been collisions and the Burgers vector norm is greater than 1.
            absB = abs.(b)
            absB /= minimum(absB)
            sqrtBType = sqrt(eps(elemT))
            bTypeArr1 = @. abs(1 - absB) < sqrtBType
            bTypeArr2 = @. abs(2 - absB) < sqrtBType
            bTypeArr3 = @. abs(3 - absB) < sqrtBType
            bType1 = sum(bTypeArr1)
            bType2 = sum(bTypeArr2)
            bType3 = sum(bTypeArr3)
            screw = 1 - nB < sqrtBType
            # Normalise Burgers vector.
            b /= nB

            nBnT_2 = nB * nT / 2
            # If it's not a screw-ish dislocation, do edge calculation and continue to next iteration.
            if !(
                bType1 == 3 && screw ||
                bType1 == 2 && bType2 == 1 && screw ||
                bType1 == 1 && bType2 == 1 && bType3 == 1 && screw
            )
                # ξ += |b|| * ||t||/2 * (ξ_climb * I + (ξ_line - ξ_climb) * t ⊗ t)
                totalDrag += nBnT_2 * (climbDrag * I3 + (lineDrag - climbDrag) * t ⊗ t)
                continue
            end

            # cos(θ) = (a ⋅ b)/(||a|| ||b||)
            # Both vectors have been normalised already.
            tDb = t ⋅ b
            cosSqθ = tDb^2
            sinSqθ = 1 - cosSqθ

            # ξ += ||b|| * ||t||/2 (ξ_screw * I + (ξ_line - ξ_screw) * t ⊗ t)
            totalDrag += nBnT_2 * (screwDrag * I3 + (lineDrag - screwDrag) * t ⊗ t)

            # If dislocation segment is pure screw skip to next iteration.
            sinSqθ < eps(typeof(tDb)) && continue

            # P, Q vectors as new local axes.
            p = b × t / sqrt(sinSqθ)
            q = p × t

            # Eqn (112) from Arsenlis et al 2007 MSMSE 15 553
            screwDragSq = screwDrag^2
            # screw_ξ_glide = 1/sqrt(sin^2(θ) / ξ_edge^2 + cos^2(θ) / ξ_screw^2)
            ScrewGlideDrag = 1 / sqrt(sinSqθ / edgeDrag^2 + cosSqθ / screwDragSq)
            # screw_ξ_climb = sqrt(sin^2(θ) * ξ_climb^2 + cos^2(θ) * ξ_screw^2)
            ScrewClimbDrag = sqrt(sinSqθ * climbDrag^2 + cosSqθ * screwDragSq)

            # ξ += ||b|| * ||t||/2 *
            #      ((screw_ξ_glide - ξ_screw) * q ⊗ q + (screw_ξ_climb - ξ_screw))
            totalDrag +=
                nBnT_2 * (
                    (ScrewGlideDrag - screwDrag) * q ⊗ q +
                    (ScrewClimbDrag - screwDrag) * p ⊗ p
                )
        end
        # Solve for velocity.
        # ξ v = f

        iNodeVel = totalDrag \ iNodeForce
        #=
        try
            iNodeVel = totalDrag \ iNodeForce
        catch SingularSystem
        origTotalDrag = totalDrag
            while true
                try
                    totalDrag += I3 * maximum(abs.(origTotalDrag)) * sqrt(eps(dType))
                    iNodeVel = totalDrag \ iNodeForce
                    break
                catch SingularSystem
                end
            end
        end =#

        nodeForce[:, node1] = iNodeForce
        nodeVel[:, node1] = iNodeVel
    end

    return network
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
using NeuralQuantum, Test
using NeuralQuantum: lookup, init_lut!, flipat_fast!, update_lut!, logψ_Δ, lut, Δ_logψ_and_∇logψ!
num_types = [Float32, Float64]

re_machines = Dict()
im_machines = Dict()

ma = (T, N) -> NDM(T, N, 2, 3)
re_machines["NDM"] = ma

ma = (T, N) -> RBMSplit(T, N, 2)
im_machines["RBMSplit"] = ma

all_machines = merge(re_machines, im_machines)

N = 4
T = last(num_types)

@testset "Test LUT: Update_LUT $name" for name=keys(all_machines)
    for T=num_types
        net = all_machines[name](T,N)
        isnothing(lookup(net)) && continue

        cnet = cached(net)

        v = state_lut(T, SpinBasis(1//2)^N, net)
        zero!(v)
        init_lut!(v, cnet, true)
        # compute exact
        results = [];
        for i=1:spacedimension(v)
            set_index!(v, i)
            init_lut!(v, cnet, true)

            for j=unique(rand(1:nsites(v), nsites(v)))
                old, new = flipat_fast!(v, j)
            end
            update_lut!(v, cnet)
            u_lut = deepcopy(lut(v))
            init_lut!(v, cnet, true)
            i_lut = lut(v)
            good = Bool[]
            for f=propertynames(i_lut)
                push!(good, getproperty(u_lut, f) ≈ getproperty(i_lut, f))
            end
            push!(results, all(good))
        end
        @test all(results)
    end
end

@testset "Test LUT: Logψ_Δ $name" for name=keys(all_machines)
    for T=num_types
        net = all_machines[name](T,N)
        isnothing(lookup(net)) && continue

        cnet = cached(net)

        vr = state(T, SpinBasis(1//2)^N, net)
        zero!(vr)

        v = state_lut(T, SpinBasis(1//2)^N, net)
        zero!(v)
        init_lut!(v, cnet)
        # compute exact
        Δ_vals = []; exact_vals = [];
        eval_vals = []; exact = [];
        for i=1:spacedimension(v)
            set_index!(vr, i)
            set_index!(v, i)
            init_lut!(v, cnet, true)

            lnψ_old = cnet(vr)
            for j=unique(rand(1:nsites(v), nsites(v)))
                old, new = flipat_fast!(v, j)
                setat!(vr, j, new)
            end
            lnψ_new  = cnet(vr)
            lnψ_Δ_ex = lnψ_new-lnψ_old
            lnψ_Δ    = logψ_Δ(cnet, v)
            push!(Δ_vals, lnψ_Δ)
            push!(exact_vals, lnψ_new-lnψ_old)
            push!(eval_vals, logψ(cnet, v))
            push!(exact, lnψ_new)
        end
        @test Δ_vals ≈ exact_vals
        @test eval_vals ≈ exact
    end
end

@testset "Test LUT: Cached Gradient $name" for name=keys(merge(re_machines, im_machines))
    for T=num_types
        net = all_machines[name](T,N)
        isnothing(lookup(net)) && continue

        cnet = cached(net)
        cder = grad_cache(net)

        vr = state(T, SpinBasis(1//2)^N, net)
        zero!(vr)

        v = state_lut(T, SpinBasis(1//2)^N, net)
        zero!(v)
        init_lut!(v, cnet)

        # compute exact
        grads = [];
        for i=1:spacedimension(v)
            set_index!(vr, i)
            set_index!(v, i)
            init_lut!(v, cnet)

            for j=1:nsites(v)
                old, new = flipat_fast!(v, j)
                setat!(vr, j, new)
            end

            der_c = ∇logψ(cnet, vr)
            der_l = ∇logψ(cnet, v)
            for f=propertynames(der_l)
                ∇  = getproperty(der_c,  f)
                c∇ = getproperty(der_l, f)
                push!(grads, ∇ ≈ c∇)
            end
        end
        @test all(grads)
    end
end

@testset "Test LUT: Cached Gradient $name" for name=keys(merge(re_machines, im_machines))
    for T=num_types
        net = all_machines[name](T,N)
        isnothing(lookup(net)) && continue

        cnet = cached(net)
        cder = grad_cache(net)

        vr = state(T, SpinBasis(1//2)^N, net)
        zero!(vr)

        v = state_lut(T, SpinBasis(1//2)^N, net)
        zero!(v)
        init_lut!(v, cnet)

        # compute exact
        grads_ex   = Bool[]
        grads_diff = Bool[]
        diffs      = Bool[]
        for i=1:spacedimension(v)
            set_index!(vr, i)
            set_index!(v, i)
            init_lut!(v, cnet)

            for j=1:nsites(v)
                old, new = flipat_fast!(v, j)
                setat!(vr, j, new)
            end

            der_c = ∇logψ(cnet, vr)
            der_l = ∇logψ(cnet, v)
            diffval_l = logψ_Δ(cnet, v)
            der_a = grad_cache(cnet)
            diffval, der_a = Δ_logψ_and_∇logψ!(der_a, cnet, v)
            push!(diffs, diffval_l ≈ diffval )
            for f=propertynames(der_l)
                ∇  = getproperty(der_c,  f)
                c∇ = getproperty(der_l, f)
                a∇ = getproperty(der_a, f)
                push!(grads_ex, ∇ ≈ c∇)
                push!(grads_diff, a∇ ≈  c∇)
            end
        end
        @test all(grads_ex)
        @test all(grads_diff)
        @test all(diffs)
    end
end


text:
eM_box = eMesh_box()

@testset "extensions" begin
    @test_throws ErrorException volume(as_tri_eMesh(eM_box))
    @test_throws ErrorException area(as_tet_eMesh(eM_box))
end


text:
module NestedArrays

import Base
import Base: undef, UndefInitializer, @propagate_inbounds
import ElasticArrays
const ea = ElasticArrays

struct NestedArray{T} <: DenseVector{Vector{T}}
    data::ea.ElasticArray{T, 2, 1}
end

function NestedArray{T}(::UndefInitializer, dim::Integer) where {T}
    data = ea.ElasticArray{T}(undef, dim, 0)
    NestedArray{T}(data)
end

function NestedArray{T}(::UndefInitializer, innerdim::Integer, outerdim::Integer) where {T}
    data = ea.ElasticArray{T}(undef, innerdim, outerdim)
    NestedArray{T}(data)
end

function NestedArray(A::AbstractMatrix{T}) where {T}
    NestedArray{T}(ea.ElasticArray(A))
end

Base.size(A::NestedArray) = size(A.data)[end:end]

Base.@propagate_inbounds Base.getindex(A::NestedArray, i::Int) = A.data[:,i]

#TODO: Implement getindex for UnitRange and Colon to return a NestedArray
Base.@propagate_inbounds function Base.getindex(A::NestedArray{T}, I::OrdinalRange) where {T}
    NestedArray{T}(A.data[:,I])
end

Base.@propagate_inbounds function Base.getindex(A::NestedArray{T}, I::Colon) where {T}
    NestedArray{T}(A.data[:,I])
end

Base.@propagate_inbounds function Base.setindex!(
    A::NestedArray{T},
    x::Vector{T},
    i::Integer,
) where {T}
    A.data[:,i] = x
end

Base.@propagate_inbounds function Base.setindex!(
    A::NestedArray{T},
    x::Matrix{T},
    I
) where {T}
    A[I] = NestedArray(x)
end

#push! and append!
function Base.push!(A::NestedArray{T}, v::Vector{T}) where {T}
    append!(A.data, v)
    A
end

end # module


text:
using Test
using BibTeX

@testset "Parsing invalid specs" begin

    @test_throws ErrorException parse_bibtex("@book")

    @test_throws ErrorException parse_bibtex("@book@")

end

@testset "Minimal functionality test" begin

    preamble, result = parse_bibtex("""
            @preamble{some instructions}
            @comment blah blah
            @string{short = long}
            @a{b,
              c = { {c} c},
              d = "d d",
              e = f # short
            }
            """);

    @test preamble == "some instructions"

    @test result["b"]["type"] == "a"

    @test result["b"]["c"] == "{ c } c"

    @test result["b"]["d"] == "d d"

    @test result["b"]["e"] == "f short"

end


text:
"""
	normalizeAndSetNorm!

# arguments:

* simulator: comparable for
  * `Quantum2dFractalSimulator`
  * `Quantum2dFractalinititerSimulator`
* tensor::Array
* normname
"""
function normalizeAndSetNorm!{T}(simulator::Quantum2dFractalSimulator{T}, tensor::Array, normname::AbstractString)
	normedTensor, normalizationFactor = normalizeTensor(tensor)
	setNorm!(simulator, normalizationFactor, normname)
	return normedTensor, normalizationFactor
end

function normalizeAndSetNorm!{T}(simulator::Quantum2dFractalInititerSimulator{T}, tensor::Array, normname::AbstractString)
	normedTensor, normalizationFactor = normalizeTensor(tensor)
	setNorm!(simulator, normalizationFactor, normname)
	return normedTensor, normalizationFactor
end


text:
const loopbackip = IPv4(127,0,0,1)

function addloopback!(block_dist)
    localip = getipaddr()
    for iplist in block_dist
        # add loopback interface for localip
        (localip in iplist) && !(loopbackip in iplist) && push!(iplist, loopbackip)
    end
    block_dist
end

function Block(hfile::HDFSFile)
    isdir(hfile) && (return Block(hfile, true, 0))

    worker_ids = workers()
    localip = getipaddr()
    worker_ips = map(x->getaddrinfo(string(isa(x, LocalProcess)?localip:x.config.host.value)), map(x->Base.worker_from_id(x), worker_ids))

    blks = hdfs_blocks(hfile)
    #logmsg(blks)
    block_dist = [map(IPv4,x[2]) for x in blks]
    addloopback!(block_dist)
    block_wrkr_ids = map(ips->worker_ids[findin(worker_ips, ips)], block_dist)
    filestat = stat(hfile)
    block_sz = filestat.block_sz
    file_sz = filestat.size

    data = [(hfile, ((x-1)*block_sz+1):(min(file_sz,x*block_sz))) for x in 1:length(block_dist)]
    Block(hfile, data, block_wrkr_ids, as_it_is, as_it_is)
end

max_stars_repo_name: tanmaykm/HadoopBlocks

max_issues_repo_name: tanmaykm/HadoopBlocks

max_forks_repo_name: tanmaykm/HadoopBlocks


text:
import Sobol.next!

"""
    A NumberGenerator must always return a vector of uniforms when called with the function next!(::NumberGenerator).
    It can contain whatever it wants as long as it does that method.
"""
abstract type NumberGenerator end

"""
    A NumberGenerator wrapper for SobolSeq.
"""
struct SobolGen <: NumberGenerator
    seq_::SobolSeq
end
function next!(number_generator::SobolGen)
    return next!(number_generator.seq_)
end

"""
    A NumberGenerator wrapper for MersenneTwister.
"""
struct Mersenne <: NumberGenerator
    twister_::MersenneTwister
    number_of_itos_::UInt
end
function next!(number_generator::Mersenne)
    return rand(number_generator.twister_, number_generator.number_of_itos_)
end

"""
    ImportanceSampler(underlying_generator_::NumberGenerator, conversion_function_::Function)
    The importance sampler takes in a NumberGenerator that is used to get a random vector of uniforms,  [0,1]^d.
    Then it inputs this vector into a function (the conversion_function): [0,1]^d -> [0,1]^d.
    Note that if the identity function is input then the numbers generated are the same as the input
    number generator. If on the other hand a convex function (in some dimension) is input it is less likely to
    get a high number for that dimension.
    If you are doing Importance Sampling remember to multiply the results by the inverse of the pdf later. This
    package doesn't do that for you.
"""
struct ImportanceSampler <: NumberGenerator
    underlying_generator_::NumberGenerator
    conversion_function_::Function
end
function next!(number_generator::ImportanceSampler)
    array_of_uniforms = next!(number_generator.underlying_generator_)
    return number_generator.conversion_function_(array_of_uniforms)
end


text:
selector(items) = @rpc select(items)

input() = @rpc input()

info(s) = @msg info(string(s))

sendnotify(s::AbstractString) = @msg notify(s)

macro !(ex)
  quote
    result = $(esc(ex))
    display(span(c(render(Inline(), $(Expr(:quote, ex))),
                   " = ",
                   render(Inline(), result))))
    result
  end
end

clearconsole() = @rpc clearconsole()

function plotsize()
  info = @rpc plotsize()
  size = if info isa Array
    info
  elseif haskey(info, "size")
    info["size"]
  else
    [401, 501]
  end
  size .- 1
end

ploturl(url::String) = @msg ploturl(url)

const SELECTORS = Dict(
  "number" => ["syntax--constant", "syntax--numeric", "syntax--julia"],
  "symbol" => ["syntax--constant", "syntax--other", "syntax--symbol", "syntax--julia"],
  "string" => ["syntax--string", "syntax--quoted", "syntax--double", "syntax--julia"],
  "macro" => ["syntax--support", "syntax--function", "syntax--macro", "syntax--julia"],
  "keyword" => ["syntax--keyword", "syntax--control", "syntax--julia"],
  "funccall" => ["syntax--support", "syntax--function", "syntax--julia"],
  "funcdef" => ["syntax--entity", "syntax--name", "syntax--function"],
  "operator" => ["syntax--operator", "syntax--julia"],
  "comment" => ["syntax--comment", "syntax--julia"],
  "variable" => ["syntax--julia"],
  "type" => ["syntax--support", "syntax--type", "syntax--julia"]
)

function syntaxcolors(selectors = SELECTORS)
  colorsstr = @rpc syntaxcolors(selectors)
  colors = Dict{String, UInt32}()
  for (k, v) in colorsstr
    colors[k] = parse(UInt32, v, base=16)
  end
  colors
end