text:
struct OffsetPrecalc{T,N,C,B,R,X,M,P<:AbstractStridedPointer{T,N,C,B,R,X,M},I} <: AbstractStridedPointer{T,N,C,B,R,X,M}
    ptr::P
    precalc::I
end
@inline Base.pointer(ptr::OffsetPrecalc) = pointer(ptr.ptr)
@inline Base.similar(ptr::OffsetPrecalc, p::Ptr) = OffsetPrecalc(similar(ptr.ptr, p), ptr.precalc)
@inline pointerforcomparison(p::OffsetPrecalc) = pointerforcomparison(p.ptr)
# @inline pointerforcomparison(p::OffsetPrecalc, i::Tuple) = pointerforcomparison(p.ptr, i)
@inline offsetprecalc(x, ::Any) = x
@inline pointerforcomparison(p::AbstractStridedPointer) = pointer(p)
@inline pointerforcomparison(p::AbstractStridedPointer, i) = gep(p, i)
@inline ArrayInterface.offsets(p::OffsetPrecalc) = offsets(p.ptr)

@inline Base.strides(p::OffsetPrecalc) = strides(p.ptr)

@inline function similar_no_offset(sptr::OffsetPrecalc{T}, ptr::Ptr{T}) where {T}
    OffsetPrecalc(similar_no_offset(sptr.ptr, ptr), sptr.precalc)
end

"""

Basically:

if I ∈ [3,5,7,9]
    c[(I - 1) >> 1]
else
    b * I
end

because

c = b .* [3, 5, 7, 9]

"""
@generated function lazymul(::StaticInt{I}, b, c::Tuple{Vararg{Any,N}}) where {I,N}
    Is = (I - 1) >> 1
    if (isodd(I) && 1 ≤ Is ≤ N) && (c.parameters[Is] !== nothing)
        Expr(:block, Expr(:meta,:inline), Expr(:call, :getindex, :c, Is))
    elseif (I ∈ (6, 10) && (I >> 2 ≤ N)) && (c.parameters[I >> 2] !== nothing)
        Expr(:block, Expr(:meta,:inline), Expr(:call, :lazymul, Expr(:call, Expr(:curly, :StaticInt, 2)), Expr(:call, :getindex, :c, I >> 2)))
    else
        Expr(:block, Expr(:meta,:inline), Expr(:call, :lazymul, Expr(:call, Expr(:curly, :StaticInt, I)), :b))
    end
end
@inline lazymul(a, b, c) = lazymul(a, b)
@inline lazymul(a::StaticInt, b, ::Nothing) = lazymul(a, b)

@generated function lazymul_no_promote(::Type{T}, ::StaticInt{I}, b, c::Tuple{Vararg{Any,N}}) where {T, I, N}
    Is = (I - 1) >> 1
    if (isodd(I) && 1 ≤ Is ≤ N) && (c.parameters[Is] !== nothing)
        Expr(:block, Expr(:meta,:inline), Expr(:call, :getindex, :c, Is))
    elseif (I ∈ (6, 10) && (I >> 2 ≤ N)) && (c.parameters[I >> 2] !== nothing)
        Expr(:block, Expr(:meta,:inline), Expr(:call, :lazymul_no_promote, Expr(:call, Expr(:curly, :StaticInt, 2)), Expr(:call, :getindex, :c, I >> 2)))
    else
        Expr(:block, Expr(:meta,:inline), Expr(:call, :lazymul_no_promote, Expr(:call, Expr(:curly, :StaticInt, I)), :b))
    end
end
@inline lazymul_no_promote(::Type{T}, a, b, c) where {T} = lazymul_no_promote(T, a, b)
@inline lazymul_no_promote(::Type{T}, a::StaticInt, b, ::Nothing) where {T} = lazymul_no_promote(T, a, b)

# @inline tdotc(a::Tuple{StaticInt{I1},Vararg}, b::Tuple{I2,Vararg}, c::Tuple{I3,Vararg}) where {I1,I2,I3} = lazymul_maybe_cached(StaticInt{I1}(),b[1],c[1])
# @inline tdotc(a::Tuple{StaticInt{I1},Vararg}, b::Tuple{I2,Vararg}, c::Tuple{Nothing,Vararg}) where {I1,I2} = lazymul(StaticInt{I1}(),b[1])
# @inline tdotc(a::Tuple{I1,Vararg}, b::Tuple{I2,Vararg}, c::Tuple{I3,Vararg}) where {I1,I2,I3} = lazymul(a[1],b[1])

# @inline tdotc(a::Tuple{StaticInt{I1},I4,Vararg}, b::Tuple{I2,I5,Vararg}, c::Tuple{I3,I6,Vararg}) where {I1,I2,I3,I4,I5,I6} = lazyadd(lazymul_maybe_cached(StaticInt{I1}(),b[1],c[1]), tdotc(Base.tail(a), Base.tail(b), Base.tail(c)))
# @inline tdotc(a::Tuple{StaticInt{I1},I4,Vararg}, b::Tuple{I2,I5,Vararg}, c::Tuple{Nothing,I6,Vararg}) where {I1,I2,I4,I5,I6} = lazyadd(lazymul(StaticInt{I1}(),b[1]), tdotc(Base.tail(a), Base.tail(b), Base.tail(c)))
# @inline tdotc(a::Tuple{I1,I4,Vararg}, b::Tuple{I2,I5,Vararg}, c::Tuple{I3,I6,Vararg}) where {I1,I2,I3,I4,I5,I6} = lazyadd(lazymul(a[1],b[1]), tdotc(Base.tail(a), Base.tail(b), Base.tail(c)))

unwrap(::Type{StaticInt{N}}) where {N} = N
unwrap(_) = nothing
# descript is a tuple of (unrollfactor) for each ind; if it shouldn't preallocate, unrollfactor may be set to 1
function precalc_quote_from_descript(descript, contig, X)
    precalc = Expr(:tuple)
    anyprecalcs = anydynamicprecals = false
    pstrideextracts = Expr(:block)
    for (i,uf) ∈ enumerate(descript)
        if i == contig || uf < 3
            push!(precalc.args, nothing)
        else
            t = Expr(:tuple)
            Xᵢ = unwrap(X[i])
            anyprecalcs = true
            if Xᵢ === nothing
                anydynamicprecals = true
                pstride_i = Symbol(:pstride_, i)
                push!(pstrideextracts.args, Expr(:(=), pstride_i, Expr(:ref, :pstride, i)))
                foreach(u -> push!(t.args, Expr(:call, :vmul_fast, u, pstride_i)), 3:2:uf)
            else
                foreach(u -> push!(t.args, u * Xᵢ), 3:2:uf)
            end
            push!(precalc.args, t)
        end
    end
    q = Expr(:block, Expr(:meta,:inline))
    if anydynamicprecals
        push!(q.args, :(pstride = strides(p)))
        push!(q.args, pstrideextracts)
    end
    if anyprecalcs
        push!(q.args, Expr(:call, :OffsetPrecalc, :p, precalc))
    else
        push!(q.args, :p)
    end
    q
end
@generated function offsetprecalc(p::AbstractStridedPointer{T,N,C,B,R,X,O}, ::Val{descript}) where {T,N,C,B,R,X,O,descript}
    precalc_quote_from_descript(descript, C, X.parameters)
end

# offsetprecalc(p::StridedBitPointer) = p

# @inline stride1offset(ptr::OffsetPrecalc, i) = stride1offset(ptr.ptr, i)

# @inline stridedoffset(ptr::OffsetPrecalc{T,<:AbstractColumnMajorStridedPointer}, ::Tuple{}) where {T} = Zero()
# @inline stridedoffset(ptr::OffsetPrecalc{T,<:AbstractColumnMajorStridedPointer{T}}, i::Tuple{I}) where {I,T} = Zero()
# # @inline stridedoffset(ptr::OffsetPrecalc{T,<:AbstractColumnMajorStridedPointer{T}}, i::Integer) where {T} = vmulnp(static_sizeof(T), i)
# @inline function stridedoffset(ptr::OffsetPrecalc{T,<:AbstractColumnMajorStridedPointer{T}}, i::Tuple{I1,I2,Vararg}) where {I1,I2,T}
#     @inbounds tdotc(Base.tail(i), ptr.ptr.strides, ptr.precalc)
# end

# @inline offsetprecalc(p::AbstractRowMajorStridedPointer{T,0}, ::Val{<:Any}) where {T} = p
# @generated function offsetprecalc(p::AbstractRowMajorStridedPointer, ::Val{descript}) where {descript}
#     precalc_quote_from_descript(Base.tail(reverse(descript)))
# end

# @inline function stridedoffset(ptr::OffsetPrecalc{T,<:AbstractRowMajorStridedPointer{T,1}}, i::Tuple{I1,I2}) where {I1,I2,T}
#     @inbounds vmul_maybe_cached(i[1], ptr.ptr.strides[1], ptr.precalc[1])
# end
# @inline function stridedoffset(ptr::OffsetPrecalc{T,<:AbstractRowMajorStridedPointer{T}}, i::Tuple{I1,I2,I3,Vararg}) where {I1,I2,I3,T}
#     ri = reverse(i)
#     @inbounds tdotc(Base.tail(ri), ptr.ptr.strides, ptr.precalc)
# end

# @generated function offsetprecalc(p::AbstractSparseStridedPointer, ::Val{descript}) where {descript}
#     precalc_quote_from_descript(descript)
# end
# @inline stridedoffset(ptr::OffsetPrecalc{T,<:AbstractSparseStridedPointer{T}}, i::Tuple{I}) where {I,T} = @inbounds vmul_maybe_cached(i[1], ptr.ptr.strides[1], ptr.precalc[1])
# # @inline stridedoffset(ptr::OffsetPrecalc{T,<:AbstractSparseStridedPointer{T}}, i::Integer) where {T} = @inbounds vmul_maybe_cached(i, ptr.ptr.strides[1], ptr.precalc[1])
# @inline function stridedoffset(ptr::OffsetPrecalc{T,<:AbstractSparseStridedPointer{T}}, i::Tuple{I1,I2,Vararg}) where {I1,I2,T}
#     @inbounds tdotc(i, ptr.ptr.strides, ptr.precalc)
# end

@inline tdot(ptr::OffsetPrecalc{T}, a, b, c) where {T} = tdot(T, a, b, ptr.precalc, c)


text:
using Test, Pkg
using DistStat
using LinearAlgebra

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    using CUDA
    ArrayType = CuArray
else
    ArrayType = Array
end

type = [Float64, Float32, Int64, Int32]

for T in type
    data = ArrayType{T}(reshape(collect(1:25), 5, 5))
    a = distribute(data)
    b = MPIArray{eltype(a), 2, ArrayType}(undef, 1, 5)
    c = ArrayType{eltype(a)}(undef, 5)

    cols = b.partitioning[DistStat.Rank()+1][2]

    @test data[:, cols] == a.localarray

    val = 500

    fill_diag!(a, val)
    data[diagind(data)] .= val

    @test data[:, cols] == a.localarray

    DistStat.diag!(b, a)
    DistStat.diag!(c, a)

    @test all(b.localarray .== val)
    @test all(c .== val)
end

ext: jl


text:
@testset "Samplers" begin

@testset "CameraSample" begin
    s = GRAYTR.CameraSample(4.0, 0.5, 0.75, 1.0)
    t = normalize(s, 2.0, 0.5)
    @test t.imgX ≈ 2.00
    @test t.imgY ≈ 1.0
    @test t.lensU ≈ 0.75
    @test t.lensV ≈ 1.0
end

@testset "StratifiedSampler" begin
    @testset "Creating" begin
        S = GRAYTR.StratifiedSampler(64, 64, 2)
        @test S.xstart == 1
        @test S.xend == 64
        @test S.ystart == 1
        @test S.yend == 64
        @test S.xs == 2
        @test S.ys == 2
        @test S.jitter == true
        @test S.xnorm == 64
        @test S.ynorm == 64
    end

    @testset "Subwindow computation" begin
        S = GRAYTR.StratifiedSampler(32, 64, 2)
        @test GRAYTR.compute_subwindow(S, 1, 4) == (1, 16, 1, 32)
        @test GRAYTR.compute_subwindow(S, 2, 4) == (17, 32, 1, 32)
        @test GRAYTR.compute_subwindow(S, 3, 4) == (1, 16, 33, 64)
        @test GRAYTR.compute_subwindow(S, 4, 4) == (17, 32, 33, 64)

        S = GRAYTR.StratifiedSampler(2, 2, 2)
        @test GRAYTR.compute_subwindow(S, 1, 4) == (1, 1, 1, 1)
        @test GRAYTR.compute_subwindow(S, 2, 4) == (2, 2, 1, 1)
        @test GRAYTR.compute_subwindow(S, 3, 4) == (1, 1, 2, 2)
        @test GRAYTR.compute_subwindow(S, 4, 4) == (2, 2, 2, 2)

        #S = GRAYTR.StratifiedSampler(2, 4, 2)
        #@test GRAYTR.compute_subwindow(S, 1, 4) == (1, 1, 1, 2)
        #@test GRAYTR.compute_subwindow(S, 2, 4) == (2, 2, 1, 2)
        #@test GRAYTR.compute_subwindow(S, 3, 4) == (1, 1, 3, 4)
        #@test GRAYTR.compute_subwindow(S, 4, 4) == (2, 2, 3, 4)
        #
        #S = GRAYTR.StratifiedSampler(4, 2, 2)
        #@test GRAYTR.compute_subwindow(S, 1, 4) == (1, 2, 1, 1)
        #@test GRAYTR.compute_subwindow(S, 2, 4) == (3, 4, 1, 1)
        #@test GRAYTR.compute_subwindow(S, 3, 4) == (1, 2, 2, 2)
        #@test GRAYTR.compute_subwindow(S, 4, 4) == (3, 4, 2, 2)

    end

    @testset "Subsampler" begin
        S = GRAYTR.StratifiedSampler(32, 64, 2)
        Sub = GRAYTR.get_subsampler(S, 2, 4)
        @test Sub.xstart == 17
        @test Sub.xend == 32
        @test Sub.ystart == 1
        @test Sub.yend == 32
        @test Sub.xs == 2
        @test Sub.ys == 2
        @test Sub.jitter == true
        @test Sub.xnorm == 32
        @test Sub.ynorm == 64

    end

    @testset "Finished state" begin
        S = GRAYTR.StratifiedSampler(4, 4, 2)
        @test !GRAYTR.finished(S, 2)
        @test !GRAYTR.finished(S, 15)
        @test GRAYTR.finished(S, 16)
        @test GRAYTR.finished(S, 17)
    end

    @testset "Getting samples" begin
        S = GRAYTR.StratifiedSampler(4, 4, 2)
        state = 0
        samples, state = GRAYTR.get_samples(S, state)

        @test length(samples) == 4
        @test eltype(samples) == GRAYTR.CameraSample
        @test state == 1
    end

end

end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
@eval begin
    export $fname

    function $fname(nlp :: AbstractNLPModel,
                    x0 :: Array{Float64,1};
                    kwargs... )
        n = nlp.meta.nvar

        c = $c

        ARCK = TrustRegion(10.0 )

        x, f, gNorm, gnorm2, iter, calls, optimal, status = TRARC(nlp,
                                                                  x0,
                                                                  ARCK,
                                                                  c;
                                                                  kwargs...)

        return (x, f, gNorm, iter, optimal, status == :UserLimit, status)
    end
    #
    function $fname(nlp :: AbstractNLPModel; kwargs ...)
        x0 = copy(nlp.meta.x0)
        return  $fname(nlp, x0; kwargs... )
    end
end


text:
import RCall.sexp
import Base.get

# struct JuliaObjectID
#     i :: Int32
# end
#
# function nextid(id :: JuliaObjectID)
#     JuliaObjectID(id.i + 1)
# end
#
# JuliaObjectID() = JuliaObjectID(0)

struct JuliaObject
    id :: Ptr{RCall.ExtPtrSxp}
    typ :: String
    JuliaObject(id :: Ptr{RCall.ExtPtrSxp}, typ = "Regular") = new(id, typ)
end

function getPlainID(x :: JuliaObject)
    x.id
end

function getType(x :: JuliaObject)
    x.typ
end

# mutable struct JuliaObjectContainer
#     object_dict :: Dict{JuliaObjectID, Any}
#     ind :: JuliaObjectID
# end
#
# JuliaObjectContainer() = JuliaObjectContainer(Dict(), JuliaObjectID())
#
# function add!(container :: JuliaObjectContainer, x, typ = "Regular")
#     container.ind = nextid(container.ind)
#     container.object_dict[container.ind] = x
#     JuliaObject(container.ind, typ)
# end

# function remove!(container :: JuliaObjectContainer, id)
#     delete!(container.object_dict, id)
# end

# function get(container :: JuliaObjectContainer, id :: JuliaObjectID)
#     container.object_dict[id]
# end
#
# function get(container :: JuliaObjectContainer, id)
#     container.object_dict[JuliaObjectID(id)]
# end

## As long as the interface stays the same, the following code should be fine.
## The global JuliaObjectContainer julia_object_stack

# const julia_object_stack = JuliaObjectContainer()

function new_obj(obj, typ = "Regular")
    # add!(julia_object_stack, obj, typ)

    # wrap in a `Ref`
    refj = Ref(obj)
    jptr = pointer_from_objref(refj)
    s = RCall.makeExternalPtr(jptr)
    RCall.jtypExtPtrs[s] = refj
    RCall.registerCFinalizerEx(s)

    JuliaObject(s, typ)
end

## debugging for freeing
# function rm_obj(id)
#     RCall.decref_extptr(RCall.sexp(x))
# end

JuliaObject(x :: JuliaObject, typ = "Regular") = x
JuliaObject(x :: RObject, typ = "Regular") = new_obj(rcopy(x), typ)
JuliaObject(x :: RCall.Sxp, typ = "Regular") = new_obj(rcopy(x), typ)
JuliaObject(x, typ = "Regular") = new_obj(x, typ)

## Conversion related to JuliaObject

const makeJuliaObjectInR = reval("JuliaCall:::juliaobjectnew")

function sexp(x :: JuliaObject)
    RCall.rcall_p(makeJuliaObjectInR, getPlainID(x), getType(x))
end

import RCall.rcopy

function rcopy(::Type{JuliaObject}, x::Ptr{EnvSxp})
    try
        # get(julia_object_stack, rcopy(Int32, RObject(x)[:getID]()))
        # get(julia_object_stack, rcopy(Int32, RObject(x)[:id]))
        RCall.jtypExtPtrs[RCall.sexp(x[:id])].x
    catch e
        nothing
    end
end

import RCall: RClass, rcopytype

rcopytype(::Type{RClass{:JuliaObject}}, x::Ptr{EnvSxp}) = JuliaObject

## Fallback conversions

@static if julia07
    import RCall.sexpclass
    sexpclass(x) = RClass{:JuliaObject}
    sexp(::Type{RClass{:JuliaObject}}, x) = sexp(JuliaObject(x))
else
    sexp(x) = sexp(JuliaObject(x))
end

## Regarding to issue #12, #13 and #16,
## we should use JuliaObject for general AbstractArray
@static if julia07
    @suppress_err begin
        JuliaCall.sexpclass(x :: AbstractArray{T}) where {T} = RClass{:JuliaObject}
    end

    ## AbstractArray{Any} should be converted to R List
    sexpclass(x :: AbstractArray{Any}) = RClass{:list}
else
    @suppress_err begin
        JuliaCall.sexp(x :: AbstractArray{T}) where {T} = sexp(JuliaObject(x))
    end

    ## AbstractArray{Any} should be converted to R List
    sexp(x :: AbstractArray{Any}) = sexp(VecSxp, x)
end

## Preserve BigFloat precision,
## as the design decision in issue #16
# sexp(x::AbstractArray{BigFloat}) = sexp(JuliaObject(x))
# sexp(x::BigFloat) = sexp(JuliaObject(x))

function setfield1!(object, name, value1)
    value = fieldtype(typeof(object), name)(value1)
    setfield!(object, name, value)
end

max_stars_repo_name: tinagui/JuliaCall

max_issues_repo_name: tinagui/JuliaCall

max_forks_repo_name: tinagui/JuliaCall


text:
using GLMakie

x_max = 4
δ = 0.02
x = inv(x_max):δ:x_max
y = inv.(x)
axis = (;
  xlabel="Number of qubits",
  xlabelsize=50,
  ylabel="Circuit depth/time",
  ylabelsize=50,
  xticks=[0],
  yticks=[0],
)

f = lines(x, y;
  axis=axis,
  linewidth=10,
)

xlims!(f.axis, 0, x_max)
ylims!(f.axis, 0, x_max)

textsize = 50

text!(f.axis, "Classical\ncomputer"; position=Point2f(0.3, 0.3), textsize=textsize)
text!(f.axis, "Quantum\ncomputer"; position=Point2f(0.5x_max, 0.5x_max), textsize=textsize)
text!(f.axis, "Full state"; position=Point2f(0.3, 0.5x_max), rotation=π/2, textsize=textsize)
text!(f.axis, "Tensor networks"; position=Point2f(0.5x_max, 0.01), textsize=textsize)

arrows!(f.axis, [0.5x_max], [0.3], [0.0], [0.5]; linewidth=15, arrowsize=50)
text!(f.axis, "Noise"; position=Point2f(0.53x_max, 0.5), textsize=textsize, align=(:left, :bottom))

display(f)

#save("test.png", f)


text:
# Generic newton method

function newton(F_DF, x₀; tol::Real=1e-12, maxiter::Int=15, verbose::Bool=true)
    verbose && return _newton_verbose(F_DF, x₀, tol, maxiter)
    return _newton_silent(F_DF, x₀, tol, maxiter)
end

function _newton_silent(F_DF, x₀::Number, tol, maxiter)
    F, DF = F_DF(x₀)
    nF = abs(F)
    nF ≤ tol && return x₀, true
    x = x₀
    i = 1
    while i ≤ maxiter && isfinite(nF)
        x -= DF \ F
        F, DF = F_DF(x)
        nF = abs(F)
        nF ≤ tol && return x, true
        i += 1
    end
    return x, false
end

function _newton_silent(F_DF, x₀, tol, maxiter)
    F, DF = F_DF(x₀)
    nF = norm(F)
    nF ≤ tol && return x₀, true
    x = copy(x₀)
    i = 1
    while i ≤ maxiter && isfinite(nF)
        x .-= ldiv!(lu!(DF), F)
        F, DF = F_DF(x)
        nF = norm(F)
        nF ≤ tol && return x, true
        i += 1
    end
    return x, false
end

function _newton_verbose(F_DF, x₀::Number, tol, maxiter)
    _display_newton_infos(tol, maxiter)
    F, DF = F_DF(x₀)
    nF = abs(F)
    AF = DF \ F
    nAF = abs(AF)
    _display_newton_iteration(0, nF, nAF)
    nF ≤ tol && return x₀, true
    x = x₀
    i = 1
    while i ≤ maxiter && isfinite(nF)
        x -= AF
        F, DF = F_DF(x)
        nF = abs(F)
        AF = DF \ F
        nAF = abs(AF)
        _display_newton_iteration(i, nF, nAF)
        nF ≤ tol && return x, true
        i += 1
    end
    return x, false
end

function _newton_verbose(F_DF, x₀, tol, maxiter)
    _display_newton_infos(tol, maxiter)
    F, DF = F_DF(x₀)
    nF = norm(F)
    ldiv!(lu!(DF), F)
    nAF = norm(F)
    _display_newton_iteration(0, nF, nAF)
    nF ≤ tol && return x₀, true
    x = copy(x₀)
    i = 1
    while i ≤ maxiter && isfinite(nF)
        x .-= F
        F, DF = F_DF(x)
        nF = norm(F)
        ldiv!(lu!(DF), F)
        nAF = norm(F)
        _display_newton_iteration(i, nF, nAF)
        nF ≤ tol && return x, true
        i += 1
    end
    return x, false
end

# in-place version

function newton!(F_DF!, x₀, F, DF; tol::Real=1e-12, maxiter::Int=15, verbose::Bool=true)
    verbose && return _newton_verbose!(F_DF!, x₀, F, DF, tol, maxiter)
    return _newton_silent!(F_DF!, x₀, F, DF, tol, maxiter)
end

function _newton_silent!(F_DF!, x₀, F, DF, tol, maxiter)
    F_DF!(F, DF, x₀)
    nF = norm(F)
    nF ≤ tol && return x₀, true
    i = 1
    while i ≤ maxiter && isfinite(nF)
        x₀ .-= ldiv!(lu!(DF), F)
        F_DF!(F, DF, x₀)
        nF = norm(F)
        nF ≤ tol && return x₀, true
        i += 1
    end
    return x₀, false
end

function _newton_verbose!(F_DF!, x₀, F, DF, tol, maxiter)
    _display_newton_infos(tol, maxiter)
    F_DF!(F, DF, x₀)
    nF = norm(F)
    ldiv!(lu!(DF), F)
    nAF = norm(F)
    _display_newton_iteration(0, nF, nAF)
    nF ≤ tol && return x₀, true
    i = 1
    while i ≤ maxiter && isfinite(nF)
        x₀ .-= F
        F_DF!(F, DF, x₀)
        nF = norm(F)
        ldiv!(lu!(DF), F)
        nAF = norm(F)
        _display_newton_iteration(i, nF, nAF)
        nF ≤ tol && return x₀, true
        i += 1
    end
    return x₀, false
end

#

function _display_newton_infos(tol, maxiter)
    println(string("\n Newton's method: tol = ", tol, ", maxiter = ", maxiter))
    println("      iteration        |F(x)|             |AF(x)|")
    println("-------------------------------------------------")
end

_display_newton_iteration(i, nF, nAF) = @printf("%11d %19.4e %19.4e\n", i, nF, nAF)


text:
# ---
# title: 203. Remove Linked List Elements
# id: problem203
# author: zhwang
# date: 2022-02-16
# difficulty: Easy
# categories: Linked List
# link: <https://leetcode.com/problems/remove-linked-list-elements/description/>
# hidden: true
# ---
#
# Remove all elements from a linked list of integers that have value **_val_**.
#
# **Example:**
#
#
#
#     Input:  1->2->6->3->4->5->6, _**val**_ = 6
#     Output: 1->2->3->4->5
#
#
#
## @lc code=start
using LeetCode

function remove_elements!(head::Union{ListNode, Nothing}, value::Int)::Union{ListNode, Nothing}
    cur = fake_head = ListNode()
    fake_head.next = head
    while !isnothing(cur.next)
        if cur.next.val == value
            cur.next = cur.next.next
        else
            cur = cur.next
        end
    end
    fake_head.next
end
## @lc code=end
## @lc test=start
@testset "203.remove-linked-list-elements.jl" begin
    @test remove_elements!(ListNode{Int}([1,2,6,3,4,5,6]), 6) == ListNode{Int}([1,2,3,4,5])
    @test remove_elements!(nothing, 1) === nothing
    @test remove_elements!(ListNode{Int}([7,7,7,7]), 7) === nothing
end
## @lc test=end

## @lc info=start
# link: [solution to 203](https://leetcode-cn.com/problems/remove-linked-list-elements/solution/203-yi-chu-lian-biao-yuan-su-python-juli-awfr/)
## @lc info=end


text:
include("header.jl")
include("simulator.jl")
include("visual.jl")

# load mech
gas = CreateSolution(mech);
const ns = gas.n_species;
const nr = gas.n_reactions;
const nu = ns + 2;
const ind_diag = diagind(ones(nu - 1, nu - 1));
const ones_nu = ones(nu - 1);

include("sensitivity.jl")

phi = 1.0;          # equivalence ratio
P = 40.0 * one_atm; # pressure, atm
T0 = 1200.0;        # initial temperature, K

@elapsed ts, pred = get_Tcurve(phi, P, T0, zeros(nr);
                                dT=dT, doplot=true, dTabort=dTabort);
# ts, pred = downsampling(ts, pred; dT=2.0, verbose=false);

npr = 16
p = zeros(npr*3);

idt = ts[end]
Tign = pred[end, end]
ng = length(ts)
np = length(p)
Fy = BandedMatrix(Zeros(ng * nu, ng * nu), (nu, nu));
Fp  = zeros(ng * nu, np);
i = 1;
i_F = 1 + (i - 1) * nu:i * nu - 1;
@view(Fy[i_F, i_F])[ind_diag] .= ones_nu;
Fy[i * nu, i * nu] = -1.0;
Fy[i * nu, (i + 1) * nu] = 1.0;
du = similar(pred[:, 1])

dts = @views(ts[2:end] .- ts[1:end - 1]) ./ idt;

function dudtp!(du, u, p, t)
    Y = @view(u[1:ns])
    T = u[end]
    mean_MW = 1.0 / dot(Y, 1 ./ gas.MW)
    ρ_mass = P / R / T * mean_MW
    X = Y2X(gas, Y, mean_MW)
    C = Y2C(gas, Y, ρ_mass)
    cp_mole, cp_mass = get_cp(gas, T, X, mean_MW)
    h_mole = get_H(gas, T, Y, X)
    S0 = get_S(gas, T, P, X)
    _p = reshape(p, npr, 3)
    kp = @. @views(exp(_p[:, 1] + _p[:, 2] * log(T) - _p[:, 3] * 4184.0 / R / T))
    qdot::typeof(p) = wdot_func(gas.reaction, T, C, S0, h_mole; get_qdot=true)
    @. qdot[1:npr] *= kp
    wdot = gas.reaction.vk * qdot
    Ydot = wdot / ρ_mass .* gas.MW
    Tdot = -dot(h_mole, wdot) / ρ_mass / cp_mass
    du .= vcat(Ydot, Tdot)
    return du
end

@threads for i = 1:ng
    u = @view(pred[:, i])
    du = similar(u)
    i_F = 1 + (i - 1) * nu:i * nu - 1
    @view(Fp[i_F, :]) .= jacobian((du, x) -> dudtp!(du, u, x, 0.0),
                                du, p)::Array{Float64,2} .* (-idt)
end

@threads for i = 2:ng
    u = @view(pred[:, i])
    du = similar(u)
    i_F = 1 + (i - 1) * nu:i * nu - 1
    @view(Fy[i_F, i_F]) .= jacobian((du, x) -> dudt!(du, x, zeros(nr), 0.0),
                                du, u)::Array{Float64,2} .* (-idt)
end

@elapsed for i = 2:ng
    u = @view(pred[:, i])
    i_F = 1 + (i - 1) * nu:i * nu - 1
    @view(Fy[i_F, i * nu]) .= - dudt!(du, u, zeros(nr), 0.0)
    @view(Fy[i_F, i_F])[ind_diag] .+= ones_nu ./ (dts[i - 1])
    @view(Fy[i_F, i_F .- nu])[ind_diag] .+= ones_nu ./ (-dts[i - 1])
    if i < ng
        Fy[i * nu, i * nu] = -1.0
        Fy[i * nu, (i + 1) * nu] = 1.0
    else
        Fy[i * nu, i * nu - 1] = 1.0
    end
end

@time dydp = - Fy \ sparse(Fp)
grad = @view(dydp[end, :]) ./ idt

max_stars_repo_name: DENG-MIT/ArrheniusActiveSubspace

max_issues_repo_name: DENG-MIT/ArrheniusActiveSubspace

max_forks_repo_name: DENG-MIT/ArrheniusActiveSubspace


text:
"""
A ProjMPOSum computes and stores the projection of an
implied sum of MPOs into a basis defined by an MPS,
leaving a certain number of site indices of each MPO
unprojected. Which sites are unprojected can be shifted
by calling the `position!` method. The MPOs used as
input to a ProjMPOSum are *not* added together beforehand;
instead when the `product` method of a ProjMPOSum is invoked,
each projected MPO in the set of MPOs is multiplied by
the input tensor one-by-one in an efficient way.

Drawing of the network represented by a ProjMPOSum
`P([H1,H2,...])`, showing the case of `nsite(P)==2`
and `position!(P,psi,4)` for an MPS `psi` (note the
sum Σⱼ on the left):

```
     o--o--o-      -o--o--o--o--o--o <psi|
     |  |  |  |  |  |  |  |  |  |  |
 Σⱼ  o--o--o--o--o--o--o--o--o--o--o Hⱼ
     |  |  |  |  |  |  |  |  |  |  |
     o--o--o-      -o--o--o--o--o--o |psi>
```
"""
mutable struct ProjMPOSum
  pm::Vector{ProjMPO}
end

ProjMPOSum(mpos::Vector{MPO}) = ProjMPOSum([ProjMPO(M) for M in mpos])

ProjMPOSum(Ms::MPO...) = ProjMPOSum([Ms...])

nsite(P::ProjMPOSum) = nsite(P.pm[1])

Base.length(P::ProjMPOSum) = length(P.pm[1])

"""
    product(P::ProjMPOSum,v::ITensor)

    (P::ProjMPOSum)(v::ITensor)

Efficiently multiply the ProjMPOSum `P`
by an ITensor `v` in the sense that the
ProjMPOSum is a generalized square matrix
or linear operator and `v` is a generalized
vector in the space where it acts. The
returned ITensor will have the same indices
as `v`. The operator overload `P(v)` is
shorthand for `product(P,v)`.
"""
function product(P::ProjMPOSum, v::ITensor)::ITensor
  Pv = product(P.pm[1], v)
  for n in 2:length(P.pm)
    Pv += product(P.pm[n], v)
  end
  return Pv
end

function Base.eltype(P::ProjMPOSum)
  elT = eltype(P.pm[1])
  for n in 2:length(P.pm)
    elT = promote_type(elT, eltype(P.pm[n]))
  end
  return elT
end

(P::ProjMPOSum)(v::ITensor) = product(P, v)

"""
    size(P::ProjMPOSum)

The size of a ProjMPOSum are its dimensions
`(d,d)` when viewed as a matrix or linear operator
acting on a space of dimension `d`.

For example, if a ProjMPOSum maps from a space with
indices `(a,s1,s2,b)` to the space `(a',s1',s2',b')`
then the size is `(d,d)` where
`d = dim(a)*dim(s1)*dim(s1)*dim(b)`
"""
Base.size(P::ProjMPOSum) = size(P.pm[1])

"""
    position!(P::ProjMPOSum, psi::MPS, pos::Int)

Given an MPS `psi`, shift the projection of the
MPO represented by the ProjMPOSum `P` such that
the set of unprojected sites begins with site `pos`.
This operation efficiently reuses previous projections
of the MPOs on sites that have already been projected.
The MPS `psi` must have compatible bond indices with
the previous projected MPO tensors for this
operation to succeed.
"""
function position!(P::ProjMPOSum, psi::MPS, pos::Int)
  for M in P.pm
    position!(M, psi, pos)
  end
end

"""
    noiseterm(P::ProjMPOSum,
              phi::ITensor,
              ortho::String)

Return a "noise term" or density matrix perturbation
ITensor as proposed in Phys. Rev. B 72, 180403 for aiding
convergence of DMRG calculations. The ITensor `phi`
is the contracted product of MPS tensors acted on by the
ProjMPOSum `P`, and `ortho` is a String which can take
the values `"left"` or `"right"` depending on the
sweeping direction of the DMRG calculation.
"""
function noiseterm(P::ProjMPOSum, phi::ITensor, dir::String)
  nt = noiseterm(P.pm[1], phi, dir)
  for n in 2:length(P.pm)
    nt += noiseterm(P.pm[n], phi, dir)
  end
  return nt
end


text:
# visuals
function visualize!(vis, p, q; Δt = 0.1)
	setvisible!(vis["/Background"], true)
	setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
	setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
	# setvisible!(vis["/Axes"], false)
	# setvisible!(vis["/Grid"], false)

	# GeometryBasics.Sphere(Point3f0(0), r)

    setobject!(vis[:satellite],
    	# Rect(Vec(-0.25, -0.25, -0.25),Vec(0.5, 0.5, 0.5)),
		GeometryBasics.Sphere(Point3f0(0), 0.25),
    	MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

    arrow_x = ArrowVisualizer(vis[:satellite][:arrow_x])
    mat = MeshPhongMaterial(color=RGBA(1.0, 0.0, 0.0, 1.0))
    setobject!(arrow_x, mat)
    settransform!(arrow_x,
    	Point(0.0, 0.0, 0.0),
    	Vec(0.75, 0.0, 0.0),
    	shaft_radius=0.05,
    	max_head_radius=0.1)

    arrow_y = ArrowVisualizer(vis[:satellite][:arrow_y])
    mat = MeshPhongMaterial(color=RGBA(0.0, 1.0, 0.0, 1.0))
    setobject!(arrow_y, mat)
    settransform!(arrow_y,
    	Point(0.0, 0.0, 0.0),
    	Vec(0.0, 0.75, 0.0),
    	shaft_radius=0.05,
    	max_head_radius=0.1)

    arrow_z = ArrowVisualizer(vis[:satellite][:arrow_z])
    mat = MeshPhongMaterial(color=RGBA(0.0, 0.0, 1.0, 1.0))
    setobject!(arrow_z, mat)
    settransform!(arrow_z,
    	Point(0.0, 0.0, 0.0),
    	Vec(0.0, 0.0, 0.75),
    	shaft_radius=0.05,
    	max_head_radius=0.1)

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	for t = 1:length(q)
		MeshCat.atframe(anim, t) do
			settransform!(vis["satellite"],
			# compose(Translation((q[t][1:3] + [-0.25; -0.25; -0.25])...),
			# 	LinearMap(quaternion_rotation_matrix(q[t][4:7]))))
			compose(Translation((q[t][1:3])...),
				LinearMap(quaternion_rotation_matrix(q[t][4:7]))))
		end
	end

    MeshCat.setanimation!(vis, anim)
end


text:
mutable struct AMGPreconditioner{Tv, Ti} <: AbstractPreconditioner{Tv,Ti}
    A::ExtendableSparseMatrix{Tv,Ti}
    fact::AlgebraicMultigrid.Preconditioner
    max_levels::Int
    max_coarse::Int
    function AMGPreconditioner{Tv,Ti}(;max_levels=10, max_coarse=10) where {Tv,Ti}
        precon=new()
        precon.max_levels=max_levels
        precon.max_coarse=max_coarse
        precon
    end
end

"""
```
AMGPreconditioner(;max_levels=10, max_coarse=10, valuetype=Float64,indextype=Int64)
AMGPreconditioner(matrix;max_levels=10, max_coarse=10)
```

Create the  [`AMGPreconditioner`](@ref) wrapping the Ruge-Stüben AMG preconditioner from [AlgebraicMultigrid.jl](https://github.com/JuliaLinearAlgebra/AlgebraicMultigrid.jl)
"""
AMGPreconditioner(;valuetype::Type=Float64, indextype::Type=Int64,kwargs...)=AMGPreconditioner{valuetype,indextype}(;kwargs...)

@eval begin
    @makefrommatrix AMGPreconditioner
end

function update!(precon::AMGPreconditioner{Tv,Ti}) where {Tv,Ti}
    @inbounds flush!(precon.A)
    precon.fact=AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(precon.A.cscmatrix))
end


text:
# indicator of the L2 norm sphere with given radius

"""
  IndSphereL2(r::Real=1.0)

Returns the function `g = ind{x : ||x|| = r}`, for a real parameter `r > 0`.
"""

immutable IndSphereL2{R <: Real} <: IndicatorNonconvex
  r::R
  function IndSphereL2(r::R)
    if r <= 0
      error("parameter r must be positive")
    else
      new(r)
    end
  end
end

IndSphereL2{R <: Real}(r::R=1.0) = IndSphereL2{R}(r)

function (f::IndSphereL2){T <: RealOrComplex}(x::AbstractArray{T})
  if abs(vecnorm(x) - f.r)/f.r > 1e-14
    return +Inf
  end
  return 0.0
end

function prox!{T <: RealOrComplex}(f::IndSphereL2, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  normx = vecnorm(x)
  if normx > 0 # zero-zero?
    scal = f.r/normx
    for k in eachindex(x)
      y[k] = scal*x[k]
    end
  else
    normy = 0.0
    for k in eachindex(x)
      y[k] = randn()
      normy += y[k]*y[k]
    end
    normy = sqrt(normy)
    y[:] *= f.r/normy
  end
  return 0.0
end

fun_name(f::IndSphereL2) = "indicator of an L2 norm sphere"
fun_dom(f::IndSphereL2) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndSphereL2) = "x ↦ 0 if ||x|| = r, +∞ otherwise"
fun_params(f::IndSphereL2) = "r = $(f.r)"

function prox_naive{T <: RealOrComplex}(f::IndSphereL2, x::AbstractArray{T}, gamma::Real=1.0)
  normx = vecnorm(x)
  if normx > 0
    y = x*f.r/normx
  else
    y = randn(size(x))
    y *= f.r/vecnorm(y)
  end
  return y, 0.0
end


text:
# TODO: move these to net.jl

### General utilities:

import Base: isequal

function isequal(a::Net, b::Net)
    for n in fieldnames(a)
        if isdefined(a,n) && isdefined(b,n)
            isequal(a.(n), b.(n)) || return false
        elseif isdefined(a,n) || isdefined(b,n)
            return false
        end
    end
    return true
end

# Do not copy everything to gpu
function gpucopy_internal(x::Net, s::ObjectIdDict)
    haskey(s,x) && return s[x]
    n = Net()
    for f in (:netinputs, :sp)
        n.(f) = x.(f)
    end
    for f in (:inputs, :outputs, :tosave, :toback, :toincr, :sparse)
        n.(f) = deepcopy_internal(x.(f),s)
    end
    for f in (:op, :out, :dif, :out0, :dif0, :tmp, :stack)
        n.(f) = gpucopy_internal(x.(f),s)
    end
    n.params = filter(x->isa(x,Par), n.op)
    return n
end

cpucsc(xtype, dims)=spzeros(xtype, dims...)
cpuarr(xtype, dims)=zeros(xtype, dims)
gpucsc(xtype, dims)=CudaSparseMatrixCSC(spzeros(xtype, dims...))
gpucsr(xtype, dims)=CudaSparseMatrixCSR(spzeros(xtype, dims...))
gpuarr(xtype, dims)=fillsync!(CudaArray(xtype, dims), 0)

stype(a)=nothing
stype(a::SparseMatrixCSC)=:csc
stype(a::CudaSparseMatrixCSC)=:csc
stype(a::CudaSparseMatrixCSR)=:csr

# function newarray(ongpu, stype, xtype, dims)
#     (ongpu ?
#      (stype==:csc ? gpucsc(xtype, dims) :
#       stype==:csr ? gpucsr(xtype, dims) :
#       stype==nothing ? gpuarr(xtype, dims) : error()) :
#      (stype==:csc ? cpucsc(xtype, dims) :
#       stype==nothing ? cpuarr(xtype, dims) : error()))
# end

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

issimilar2(i,o)=(eltype(i) == eltype(o) && size(i) == size(o))
issimilar3(i,o)=(eltype(i) == eltype(o) && size(i) == size(o) && issparse(i) == issparse(o))

nothings(n::Integer)=fillsync!(cell(n),nothing)

# ### DEBUGGING

ptr16(x)=hex(x==nothing ? 0 : hash(pointer(x)) % 0xffff, 4)
ptr8(x)=hex(x==nothing ? 0 : hash(pointer(x)) % 0xff, 2)
idx1(x)=(x==nothing ? -1 : atype(x)==CudaArray ? to_host(x)[1] : atype(x)==Array ? x[1] : error("$(typeof(x))"))

vecnorm0(x,y...)=map(vecnorm0,(x,y...))
vecnorm0(x::Vector)=map(vecnorm0,x)
vecnorm0(x::Tuple)=map(vecnorm0,x)
vecnorm0(x::Par)= ((isdefined(x,:out)? vecnorm0(x.out) : 0),
                   (isdefined(x,:dif)? vecnorm0(x.dif) : 0))
vecnorm0(::Void)=0
vecnorm0(x)=(@sprintf("%.8f",vecnorm(x))) #floor(1e6*vecnorm(x))/1e6

### DEAD CODE

# function Base.summary(r::Net)
#     """
#     op: $(map(typeof,r.op))
#     inputs: $(r.inputs)
#     outputs: $(r.outputs)
#     params: $(vecnorm0(r.params))
#     tosave: $(find(r.tosave))
#     toback: $(find(r.toback))
#     toincr: $(find(r.toincr))
#     sparse: $(find(r.sparse))
#     out: $(vecnorm0(r.out))
#     dif: $(vecnorm0(r.dif))
#     out0: $(vecnorm0(r.out0))
#     dif0: $(vecnorm0(r.dif0))
#     tmp: $(vecnorm0(r.tmp))
#     stack: $(vecnorm0(r.stack[1:r.sp]))
#     sp: $(r.sp)
#     """
# end

# function dbg(r,f,n)
#     r.dbg || return
#     a = r.(f)[n]
#     print("\n==> $f[$n]($(ptr16(a))) stack:")
#     println(map(ptr16, r.stack[1:r.sp]))
#     a != nothing && display(convert(Array,a))
#     if n <= nops(r) && !isempty(params(r.op[n]))
#         p = params(r.op[n])[1]
#         println("\nw[$n]($(ptr16(p.arr)))")
#         display(convert(Array,p.arr))
#         if isdefined(p,:diff)
#             println("\ndw[$n]($(ptr16(p.diff)))")
#             display(convert(Array,p.diff))
#         end
#     end
# end

# forw{T<:Number}(r::Net,  x::Vector{T}; a...)=error("forw expects a minibatch") # forw(r, reshape(x,  length(x),  1); a...)
# back{T<:Number}(r::Net, dy::Vector{T}; a...)=error("back expects a minibatch") # back(r, reshape(dy, length(dy), 1); a...)

# function back(r::Net, dy; o...)                 # dy[1:T] loss gradients for output y
#     (N,T) = (length(r.net), length(dy))
#     @assert length(r.x) == T
#     inithidden(r, :dh, N, 1)                    # r.dh[1:N] are loss gradients for hidden r.h[1:N,t]
#     initdiff(r)                                 # initialize weight gradients dw to zero
#     for t = T:-1:1
#         r.dh[N] = dy[t]                         # y[t]=r.h[N,t] so dy[t]=r.dh[N]
#         for n = N:-1:1
#             r.dh[n] == nothing && continue      # 'nothing' represents 0 loss gradient
#             ni = ninputs(r.net[n])
#             x = rnninput(r,n,t)
#             ni==1 && (x=x[1])                   # handle single input vs tuple input layers
#             dx = back(r.net[n], r.dh[n]; x=x, y=r.h[n,t], incr=true) # incr=true increments dw
#             ni==1 && (dx=(dx,))                 # handle single input vs tuple input layers
#             for j=1:ni                          # set r.dh using dx for each input of n
#                 i=r.inputs[n][j]                # i is the j'th input of n
#                 if (r.dh[i] == nothing)         # dh[i] should be updated using dx[j]
#                     r.dh[i] = dx[j]
#                 elseif length(r.freads[i]) == 1
#                     r.dh[i] = dx[j]
#                 else
#                     axpy!(1, dx[j], r.dh[i])    # increment dh[i] if i has more than one output
#                 end
#             end
#             if length(r.freads[n]) > 1
#                 fillsync!(r.dh[n], 0)               # if n has multiple outputs, reset dh[n]
#             end
#         end
#     end
# end

# function initparams(r::Net, inputs...; seq=false)
#     for n = 1:nops(r)
#         p = params(r.op[n])
#         if !isempty(p) && findfirst(isempty,p)>0
#             forw(r.op[n], r.out0[r.inputs[n]]...; y=r.out0[n])
#         end
#     end
#     for w in params(r)
#         similar!(w, :diff, w.arr)
#         seq && similar!(w, :inc, w.arr)
#     end
# end

# # r.buf[n] is the array pre-allocated for register r.reg[n]
# # r.buf is initialized in three steps:
# # - initbuf determines the indices for array reuse.
# #   i=0 means allocate, i>0 means reuse r.buf[i]
# # - allocbuf allocates 0 size KUdense arrays
# #   this is a separate step because we need input x for type info
# # - the actual op resize the bufs as needed.

# # DONE: this should also look at back needs when allocating

# function initbuf(r::Net)
#     N = length(r.op)
#     index = fillsync!(Array(Any,N), 0)
#     for n=1:N
#         if overwrites(r.op[n])
#             index[n] = i = r.inputs[n][1] # tentatively prepare to overwrite first input
#             k = n
#             while true # see if anybody else uses i before its next update
#                 k = mod1(k+1, N)
#                 in(i, r.inputs[k]) && (index[n] = 0; break)
#                 k == i && break
#             end # fencepost check: N==1:OK, i==n:OK
#         end
#     end
#     # info("$(sum(index.==0)) bufs will be allocated.")
#     return index
# end

# function allocbuf(r::Net, x)
#     isa(r.buf[1], Int) || return
#     N = length(r.op)
#     for n=1:N
#         if r.buf[n] == 0
#             r.buf[n] = KUdense(atype(x[1]), eltype(x[1]), (0,))
#         end
#     end
#     for n=1:N
#         isa(r.buf[n], Int) || continue
#         k = n
#         while isa(r.buf[k], Int)
#             k = r.buf[k]
#         end
#         r.buf[n] = r.buf[k]
#     end
# end

# # r.output[n] points to the last output of r.op[n]
# # its value could be an array (provided by r.buf[n]) or 'nothing'
# # 'nothing' is used to represent the zero matrix at time=0
# # and others that result from them at time>0

# function initoutput(r::Net)
#     r.output = fillsync!(Array(Any,length(r.op)), nothing)
# end

# # r.stack[n<=N][t] is a copy of the output of r.op[n] at r.time=t
# # r.stack[n>N][t] is a copy of the n-N'th MLP input at r.time=t
# # - where N is the number of op
# # - remember MLP inputs are indicated by 0,-1,-2 etc. in r.inputs
# # We only keep the copies necessary for the back calculation
# # - initstack initializes necessary r.stack[n] with empty arrays
# #   the other elements in r.stack are left unassigned
# # - allocstack makes sure at least r.time arrays are allocated
# #   for each assigned r.stack[n]
# #   allocstack also copies inputs x as necessary

# function initstack(r::Net)
#     N = length(r.op)
#     hist = Array(Any,N+ninputs(r))
#     for n=1:N
#         back_reads_y(r.op[n]) && !isassigned(hist,n) && (hist[n]=Any[])
#         if back_reads_x(r.op[n])
#             for i in r.inputs[n]
#                 i <= 0 && (i = N+1-i)
#                 !isassigned(hist,i) && (hist[i]=Any[])
#             end
#         end
#     end
#     return hist
# end

# function allocstack(r::Net, x)
#     N = length(r.op)
#     for n=1:N
#         isassigned(r.stack,n) || continue
#         length(r.stack[n]) >= r.time && continue
#         length(r.stack[n]) == r.time - 1 || error("Stack corruption")
#         push!(r.stack[n], KUdense(atype(x[1]), eltype(x[1]), (0,)))
#     end
#     for i=1:length(x)
#         n=N+i
#         isassigned(r.stack,n) || continue
#         if length(r.stack[n]) < r.time
#             length(r.stack[n]) == r.time - 1 || error("Stack corruption")
#             push!(r.stack[n], copy(x[i]))
#         else
#             copysync!(r.stack[n][r.time], x[i])
#         end
#     end
# end

# function forw1(r::Net, x; train=true, o...)
#     train ?
#     forwtrain(r, x; o...) :
#     forwpredict(r, x; o...)
# end

# function forwtrain(r::Net, x; o...)
#     r.t += 1
#     copyto(r.x, x, r.t)
#     for n = 1:N
#         r.y[n] = forw(r.net[n], forwinput(r,n,x)...; o...)
#         if !predict

#             if ((n == N) ||                     # preserve the last layer
#                 (length(r.freads[n]) > 1) ||    # preserve if multiple forw reads
#                 (length(r.breads[n]) > 0))      # preserve if any back reads
#                 copyto(r.h, h, n, r.t)
#             else
#                 r.h[n,r.t] = h
#             end
#         end
#         return r.y[N]
#     end
# end

# function copyto(a, x, i...)
#     # resize a if necessary here...
#     ai = getindex(a, i...)
#     if ai == nothing
#         setindex!(a, x, i...)
#     else
#         resize!(ai, size(x))
#         copysync!(ai, x)
#     end
#     return getindex(a, i...)
# end

# # DONE: the semantics of forw should be simpler: each unit reads the
# # values of its inputs as they are at that point in time, whether they
# # are ahead in the network or behind.

# function inithidden(r::Net, x, N, T)
#     if !isdefined(r,x)
#         r.(x) = fillsync!(Array(Any, N, T), nothing)
#     elseif size(r.(x), 1) != N
#         error("Size mismatch")
#     elseif size(r.(x), 2) != T
#         isa(r.(x), SubArray) && (r.(x) = r.(x).parent)
#         if size(r.(x), 2) > T
#             r.(x) = sub(r.(x), 1:N, 1:T)
#         elseif size(r.(x), 2) < T
#             h = Array(Any, N, T)
#             copysync!(h, 1, r.(x), 1, length(r.(x)))
#             h[length(r.(x))+1:end] = nothing
#             r.(x) = h
#         end
#     end
#     @assert size(r.(x)) == (N,T)
#     return r.(x)
# end

# OK: pass x,y
## OK: To pass x, y take a look at the current interface
## OK: Return dx and param x for back may be tuples
# OK: figure out dw increment
## OK: dw gets zeroed once in the beginning
## OK: incr=(t<T): what if at t==T dh was nothing so it didn't get zeroed out?
# OK: figure out dx increment
## OK: To do dx increment, rnn back should handle it
## OK: what if we copy same dx to two inputs and one wants to increment
## OK: dy for layers with forwreads > 1 needs to be zeroed every time step
## NO: it needs to be zeroed by its last output, which is the largest index smaller than it, or if not largest index in forwreads: remember that and zero with that, check reset array and if target nonempty increment or reset.
## OK: actually it can be just zeroed as soon as it is used.

# OK: Initialize w the first time?
# OK: back needs x, y, and incr options.
# OK: get rid of forw.y and back.dx parameters

# OK: dy should be like x, only t component, and refer to the output of net[N,t]
# OK: if dy[t] is nothing do not propagate any error.
# OK: we dont want to alloc dy per time step, we want one dy per layer, let layers manage them.

# function rnnget(rnn, x, y, t, l)
#     isa(rnn[l],Tuple) ?
#     rnnget(rnn, x, y, t, l, rnn[l][1], rnn[l][2:end]) :
#     rnnget(rnn, x, y, t, l, rnn[l], tuple((l-ninputs(rnn[l])):(l-1)...))
# end

# function rnnget(rnn, x, y, t, l, layer, index)
#     n = ninputs(layer)          # DONE: ninputs not defined (done for ninputs=1 default)
#     length(index) == n || error("Wrong number of inputs")
#     input = map(index) do i
#         i == 0 ? x[t] :
#         i < l  ? y[i,t] :
#         y[i,t-1]  # DONE: problem with t=1
#     end
#     n == 1 && (input = input[1])
#     return (layer, input)
# end

# function rnninit(rnn, x)        # DONE: avoid allocation every call: have an rnn type?
#     # the first row is reserved for the input
#     # the first column is reserved for empty arrays?
#     # or should we try to accept nothing input?
#     # but what about batch size?  should be hdim x nbatch
#     # hinton says train h0 like the weights: copy column nbatch times
#     y = fillsync!(Array(Any, 1+length(rnn), 1+length(x)), nothing)
# end

# nothings(m,n)=fillsync!(Array(Any,m,n),nothing)

# function forwreads(r::Net)
#     N = length(r.inputs)
#     freads = [ Int[] for n=1:N ]
#     for n=1:N
#         for x in r.inputs[n]
#             (x>0) && push!(freads[x], n)
#         end
#     end
#     freads
# end

# function backreads(r::Net)
#     N = length(r.inputs)
#     freads = forwreads(r)
#     breads = [ Int[] for n=1:N ]
#     for n=1:N
#         for o in freads[n]
#             back_reads_x(r.net[o]) && push!(breads[n], o)
#         end
#         back_reads_y(r.net[n]) && push!(breads[n], n)
#     end
#     breads
# end

# function rnnsaves(r::Net)
#     map(forwreads(r), backreads(r)) do f,b
#         length(f) > 1 || length(b) > 0
#     end
# end

# type RNN1; net; inputs; freads; breads; t; x; y; Y; dy;
#     function RNN1(a...)
#         r = new()
#         n = length(a)
#         r.net = map(ai->(isa(ai,Tuple) ? ai[1] : ai), a)
#         r.inputs = [ Int[] for i=1:n ] # inputs[j]={i: h[j,t] needs h[i,t']}
#         # DONE: we only need the counts not the elements for freads and breads
#         r.freads = [ Int[] for i=1:n ] # freads[i]={j: h[j,t] needs h[i,t']}
#         r.breads = [ Int[] for i=1:n ] # breads[i]={j:dy[j,t] needs h[i,t']}
#         for j=1:n
#             inputs = isa(a[j],Tuple) ? a[j][2:end] : ((j-ninputs(a[j])):(j-1))
#             ninputs(r.net[j]) == length(inputs) || error("Wrong number of inputs")
#             overwrites(r.net[j]) && in(0, inputs) && error("$j:$(typeof(r.net[j])) overwrites Net input")
#             push!(r.inputs[j], inputs...)
#             for i in r.inputs[j]
#                 (0 <= i <= n) || error("Bad input index")
#                 (i > 0) && push!(r.freads[i], j)
#             end
#         end
#         # DONE: there is no need for this separate for loop
#         for i=1:n
#             for j in r.freads[i]
#                 back_reads_x(r.net[j]) && push!(r.breads[i], j)
#             end
#             back_reads_y(r.net[i]) && push!(r.breads[i], i)
#         end
#         r.t = 0
#         return r
#     end
# end

# # NO: merge x, h?  x has too many special exceptions.
# # DONE: apply the same dependency analysis to the first layer and the last layer too.

# function forwinput(r::Net, x, n)
#     indices = r.inputs[n]
#     arrays = Array(Any, length(indices))
#     for j = 1:length(indices)
#         i = indices[j]
#         arrays[j] = ((i == 0) ?
#                      (overwrites(r.net[n],j) ? copyto(r.y, x, n) : x) :
#                      (overwrites(r.net[n],j) && (length(r.freads[i])>1) ? copyto(r.y, r.y[i], n) : r.y[i]))
#     end
#     return arrays
# end

# copy or point to x from output[n+1]?
# how to handle t=0
# NO: have r.inputs contain the actual matrices? sometimes inputs are "nothing"
# ASSUME: do we still need protection from overwrites or should the compiler handle that?
# OK: all op need to support y keyword option.
# DONE: make sure the y option is not treated optional!

        # if !predict && isassigned(r.stack,n)
        #     if r.output[n] != nothing
        #         copysync!(r.stack[n][r.time], r.output[n])
        #     else
        #         warn("Recording 'nothing' in stack at $n:$(typeof(r.op[n]))")
        #         r.stack[n][r.time] = nothing
        #     end
        # end
#     # incr w
#     # incr dy
#     # multiple dy
#     # storage for dy: initalization, increment, sharing?
#     # cannot use the same registers as the output, may need simultaneously, may have different sharing
#     # pop for previous reg values
#     # r.dy[n] for diff registers

# function forwinput(r::Net, x, n)
#     map(r.inputs[n]) do i               # return the input matrix or tuple of input matrices to produce r.output[n]
#         i <= 0     ? x[1-i] :           # i <= 0 is used to indicate input x[1-i]
#         i < n      ? r.output[i] :      # 0 < i < n are outputs from the current time slice
#         r.time > 1 ? r.output[i] :	# i >= n are from the previous time slice if r.time > 1
#         nothing                         # nothing represents zero matrix from r.time=0
#     end
# end

# DONE: pushinput takes a lot of time, do we need it for ffnn?
# Yes, we don't know if the next move is forw or back.

# DONE: l.y and l.dx for layers should point to internal storage, not externally provided ones.
# right now an op may not use the provided y as a result of resize.
# ops should not resize externally provided output arrays? (at least without a warning)

# DONE: r.dy: last op has an extra output to the outside!
# DONE: r.dy: what to do with the network inputs: alloc for now but do not use makes code simpler.
# DONE: r.dy: implement dx option for back for all ops

# DONE: add reset (reset taken, called it init)

# # registers can be (1) nothing, (2) pointer to net input, (3) pointer to net.buf

# regbuf(r,n)=(r.reg[r.out[n]]!=nothing ? r.reg[r.out[n]] : r.buf[r.out[n]])

# # r.reg[i] is the i'th register of the MLP
# # Each register points to an array allocated elsewhere
# # or is set to 'nothing' representing the zero matrix
# # These registers are used for network inputs, op outputs, op gradients.
# # initreg(r) just creates an empty array of registers r.reg
# # other init functions below add registers as needed to r.reg

# function initreg(r::Net)
#     r.reg = Any[]
# end

# function initbuf(r::Net)
#     r.buf = copy(r.reg) # array of nothings, initialized at initforw
# end

# function initregi(r::Net, i::Int, y, dims=size(y))
#     dy = r.reg[i]
#     dy == nothing && (dy = r.buf[i])
#     if dy == nothing
#         dy = r.buf[i] = similar(y, dims)
#     end
#     if size(dy) != dims
#         warn("Resizing $(size(dy))->$dims")
#         dy = r.buf[i] = resize!(dy,dims)
#     end
#     r.reg[i] == dy || r.reg[i] == nothing || error("reg $i mismatch")
#     atype(dy) == atype(y) || error("atype mismatch")
#     eltype(dy) == eltype(y) || error("eltype mismatch")
#     r.reg[i] = dy
# end

        # x = xx[1]
        # i = r.out[n]
        # y = getybuf(r,n)
        # if y == nothing
        #     y = r.buf[i] = similar(x, s)
        # end
        # if size(y) != s
        #     warn("Resizing $(size(y))->$s")
        #     y = r.buf[i] = resize!(y,s)
        # end
        # r.reg[i] == y || r.reg[i] == nothing || error("reg mismatch")
        # atype(y) == atype(x) || error("atype mismatch")
        # eltype(y) == eltype(x) || error("eltype mismatch")
            # DONE: (back) buf gets the result but reg is still nothing? no need for buf/reg distinction in back
            # DONE: (back) zero out the dw
# op(r::Net,n::Int)=r.op[n]

# gety(r::Net,n::Int)=r.reg[r.out[n]]
# getdy(r::Net,n::Int)=r.reg[r.dif[n]]
# getx(r::Net,n::Int)=r.reg[r.out[r.inputs[n]]]
# getdx(r::Net,n::Int)=r.reg[r.dif[r.inputs[n]]]

# getbuf(r::Net,i::Int)=(r.reg[i]!=nothing?r.reg[i]:r.buf[i])
# getybuf(r::Net,n::Int)=getbuf(r,r.out[n])
# getdybuf(r::Net,n::Int)=getbuf(r,r.dif[n])
# getxbuf(r::Net,n::Int)=map(i->getybuf(r,i),r.inputs[n])
# getdxbuf(r::Net,n::Int)=map(i->getdybuf(r,i),r.inputs[n])

# get1x(r::Net,n::Int)=get1(getx(r,n))
# get1dxbuf(r::Net,n::Int)=get1(getdxbuf(r,n))

# sety(r::Net,n::Int,x)=(r.reg[r.out[n]]=x)
# setdy(r::Net,n::Int,x)=(r.reg[r.dif[n]]=x)
# setinput(r::Net,n::Int,x)=(r.reg[r.out[n+nops(r)]]=x)
# getinput(r::Net,n::Int)=r.reg[r.out[n+nops(r)]]

# pushy(r::Net,n::Int)=(r.save[n] && pushreg(r,r.out[n]))
# popy(r::Net,n::Int)=(r.save[n] ? popreg(r,r.out[n]) : r.reg[r.out[n]])
# pushinput(r::Net,n::Int)=(r.save[n+nops(r)] && pushreg(r,r.out[n+nops(r)]))
# popinput(r::Net,n::Int)=(r.save[n+nops(r)] && popreg(r,r.out[n+nops(r)]))

# function pushreg(r::Net,i::Int)
#     length(r.stack) <  r.sp && error("Stack error")
#     length(r.stack) == r.sp && push!(r.stack, :newcell)
#     r.sp += 1
#     if r.reg[i] == nothing
#         r.stack[r.sp] == nothing || r.stack[r.sp] == :newcell || warn("pushing nothing")
#         r.stack[r.sp] = nothing
#     elseif r.stack[r.sp] == nothing
#         warn("copying over nothing")
#         r.stack[r.sp] = copy(r.reg[i])
#     elseif r.stack[r.sp] == :newcell
#         r.stack[r.sp] = copy(r.reg[i])
#     elseif size(r.reg[i]) != size(r.stack[r.sp])
#         warn("resizing during push")
#         copysync!(r.stack[r.sp], r.reg[i])
#     else
#         copysync!(r.stack[r.sp], r.reg[i])
#     end
# end

# function popreg(r::Net,i::Int)
#     r.sp > 0 || error("Stack error")
#     if r.reg[i] == nothing
#         r.stack[r.sp] == nothing || warn("popping array over nothing")
#     elseif r.stack[r.sp] == nothing
#         # warn("popping nothing over array")
#     elseif size(r.reg[i]) != size(r.stack[r.sp])
#         warn("resizing during pop")
#     end
#     r.reg[i] = r.stack[r.sp]
#     r.sp -= 1
#     return r.reg[i]
# end

# function initsequence(r::Net, x::Vector; trn=false, a...)
#     r.sp == 0 || error("Stack corruption")
#     inputs = isa(x[1],Tuple) ? x[1] : (x[1],)
#     initbatch(r, inputs...; trn=trn, seq=true, a...)
#     fillsync!(r.out, nothing)                               # to represent zero matrices at t=0
#     if trn
#         fillsync!(r.dif, nothing)                           # why? (todo)
#         for n=1:length(r.dif0)
#             r.toincr[n] && fillsync!(r.dif0[n], 0)           # zeroed by back at every item in sequence
#         end
#         for w in params(r)
#             fillsync!(w.diff, 0)                            # zeroed only once at the beginning of the sequence
#         end
#     end
# end

# function initbatch(r::Net, inputs...; trn=false, seq=false)
#     length(inputs) == ninputs(r) || error("Wrong number of inputs")
#     initout0(r, inputs...)
#     !seq && fillsync!(r.out, nothing)
#     if trn
#         initparams(r, inputs...; seq=seq)
#         initdif0(r)
#         if !seq
#             fillsync!(r.dif, nothing)                           # why? (todo)
#             for n=1:length(r.dif0)
#                 r.toincr[n] && fillsync!(r.dif0[n], 0)           # zeroed by back at every item in sequence
#             end
#             for w in params(r)
#                 fillsync!(w.diff, 0)                            # zeroed only once at the beginning of the sequence
#             end
#         end
#     end
# end

# DONE: incremental updates
 # DONE: forw: do we really need inputs to be tuple here? YES, have the option of toincr-input nets like multi-input ops.
# DONE: forw: (minor) this ends up using a lot of nothing pushes first minibatch, that's ok, 'nothing' is a legit value.

# Terminology:
# An instance is an item or a sequence.
# An item is a contiguous feature array: x[d...]
# A sequence is a sequence of items in time: x[t][d...]
# An item batch is a contiguous array with an extra dimension: x[d...,i]
# inputs[i] should be an item batch.

# function initbatch(r::Net, inputs...)

# end

# function init(r::Net, x::Vector; a...)
#     if isempty(x) || x[1] == nothing
#         error("Got nothing as input")
#     elseif isa(x[1], Tuple)
#         init(r, x[1]...; seq=true, a...)
#     elseif isbits(eltype(x[1]))
#         init(r, x[1]; seq=true, a...)
#     end
# end

# DONE:
# + try the adding problem with irnn
# + gradient analysis with adding, should fail
# + add incremental update of dy and confirm gradient testing is now ok
# + let forw and back write the input to their own registers
# x ops with forw=nothing should accept back=nothing input
# + minibatching
# x rethink the 'nothing' optimization, removed from dif but not out?
# + return value for back, returndx option: no need, return if dx specified.
# + add gradient clipping
# + implement train/predict and try ffnn mnist experiments: how do we treat sequences and minibatches?
# + we need to solve predict() and accuracy() problems first: use MLP for now?

### initforw:

# We have a batch or a sequence of batches as input.
# The batch could be for an FNN or an RNN.
# Before each forw-batch we make sure out0 has the right size.
# Before each back-batch we make sure dif0 has the right size.
# dw needs to be zeroed only at the beginning of a sequence.
# dw does not need incr if not part of a sequence.
# out and dif registers need to be nothinged only at the beginning of a sequence.
# we need to preserve their content from one rnn-batch to the next.
# dy for multi ops are zeroed during back so init should not worry.
# so we need to know if we are in a sequence or not: add a keyword arg to forw and back.

# The input x::Vector can be x[i][t][d...,b] or x[t][d...,b]
# where i:instance, t:time, d:dims, b:batch
# We only want to init if x[t][d...,b]
# In that case x[1] will be a numeric array
# or a tuple of numeric arrays (to support multiple inputs)

# init before a sequence zeroes out everything:
# out[:] is set to nothing to represent zero matrices at t=0
# dif[:] is set to nothing why?  (todo)
# dif0[n] is set to zero for r.multi[n]
# but first we need to allocate using the batch init
# the input is x[t][d...,i] for single input
# or x[t][input][d...,i] for multiple inputs

# at the beginning of a sequence:
# we want out[:]=nothing, dif[:]=nothing, dw=0, dif0[i] for multi.
# before every itembatch:
# we just want to fix the sizes.

# """
# Life cycle of registers:

# 		out	out0	dif	dif0	tmp	dw	stack	calls
# Net		-	-	-	-	-	-	-	initout,initdif
# initout		n	[]	-	-	-	-	-
# initdif		-	-	n	[]	[]+m	-	-
#  forwseq	-	-	-	-	-	-	-	initseq,forw+s*
#   initseq	n	-	n+t	0+t+m	-	0+t	-	initbat
#  forw		wr	w	-	-	-	-	w	initbat,seq passed to initbat
#   initbat	n-s	-	n+t-s	0+t-s+m	-	0+t-s	-	initpar2(trn),initdif0(trn),initout0,seq passed to initpar2
#    initout0	-	a	-	-	-	-	-
#    initpar2	-	r	-	-	-	a	-
#    initdif0	-	r	-	a	a+m	-	-
#  backseq	-	-	-	-	-	-	-	back+s*
#   back		rw(pop)	-	nwr0	w	w	-	r	seq passed as incr to op.back
#  loss		r	-	-	-	-	-	-

# (+m:multi, +s:seq, +t:trn, r:read, w:write, a:alloc, 0:zero-out, n:fill-nothing)
# """

# # Not init w here any more...
#     # for w in params(r)
#     #     if !isdefined(w,:diff)
#     #         if issparse(x)      # todo: is this the right thing for all ops?
#     #             w.diff = CudaSparseMatrixCSR(spzeros(eltype(x), size(w)...))
#     #         else
#     #             w.diff = similar(w.arr)
#     #         end
#     #     end
#     # end
#     # if seq
#     #     for w in params(r)
#     #         similar!(w, :inc, w.diff)
#     #         fillsync!(w.diff, 0)                            # zeroed only once at the beginning of the sequence
#     #     end
#     # end

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

# # Do not copy inputs to gpu:

# function gpucopy_internal(x::Net, stackdict::ObjectIdDict)
#     if haskey(stackdict, x)
#         return stackdict[x]
#     end
#     y = ccall(:jl_new_struct_uninit, Any, (Any,), Net)
#     stackdict[x] = y
#     for i in fieldnames(x)
#         if isdefined(x,i)
#             y.(i) = (i == :inputs ?
#                      cpucopy_internal(x.(i), stackdict) :
#                      gpucopy_internal(x.(i), stackdict))
#         end
#     end
#     return y
# end

    # if keepstate
    #     copysync!(f.out, f.out0)   ## why did we copy here?
    # else
    #     fillsync!(f.out, nothing)
    # end
    # fillsync!(f.dif, nothing)
    # for n=1:length(f.op)    ## todo: fix reset dif for back
    #     f.toback[n] && f.toincr[n] && f.dif0[n]!=nothing && fillsync!(f.dif0[n], 0)
    # end
# ### Stack functions: push, pop

# function push(f::Net,r::Reg)
#     length(f.stack) <  f.sp && error("Stack error")
#     length(f.stack) == f.sp && push!(f.stack, :newcell)
#     f.sp += 1
#     if r.out == nothing                             # todo: (minor) remove these checks once code is tested
#         f.stack[f.sp] == nothing || f.stack[f.sp] == :newcell || Base.warn_once("pushing nothing over array")
#         f.stack[f.sp] = nothing
#     elseif f.stack[f.sp] == nothing
#         Base.warn_once("pushing array over nothing")  # This actually happens with rnnlm keepstate
#         f.stack[f.sp] = copy(r.out)
#     elseif f.stack[f.sp] == :newcell
#         f.stack[f.sp] = copy(r.out)
#     elseif size(r.out) != size(f.stack[f.sp])
#         Base.warn_once("pushing array of different size")
#         resize!(f.stack[f.sp], size(r.out))
#         copysync!(f.stack[f.sp], r.out)
#     else
#         copysync!(f.stack[f.sp], r.out)
#     end
# end

# function pop(f::Net,r::Reg)
#     f.sp > 0 || error("Stack error")
#     if r.out == nothing
#         f.stack[f.sp] == nothing || warn("popping array over nothing")
#     elseif f.stack[f.sp] == nothing
#         # warn("popping nothing over array")
#     elseif size(r.out) != size(f.stack[f.sp])
#         warn("popping different sized array")
#     end
#     r.out = f.stack[f.sp]
#     f.sp -= 1
# end

        # get(r,:incr) && # TODO: incr not ready here, can we do something more efficient than zeroing everything?
        # r.dif0!=nothing &&
        # isdefined(r,:dif0) && get(r,:incr) && fillsync!(r.dif0, 0)


text:
function molecule(::Molecule"formaldehyde")
  return """
O	0.0000000	0.0000000	0.6744930
C	0.0000000	0.0000000	-0.5297240
H	0.0000000	0.9347280	-1.1087990
H	0.0000000	-0.9347280	-1.1087990
  """
end


text:
module MeshBasisFunctions

using Meshes
using StaticArrays
using SparseArrays

export
    MeshBasisFunction,
    PiecewiseLinearBasis,
    BarycentricCoordTransform,
    tobarycoord,
    tonatcoord,
    interpolation_wts,
    interpolation_operator,
    interpolate

abstract type MeshBasisFunction{D} end

struct PiecewiseLinearBasis{D} <: MeshBasisFunction{D} end

struct BarycentricCoordTransform{D, F, L}
    p0::Point{D, F}
    T::SMatrix{D, D, F, L}

    function BarycentricCoordTransform(p0::Point{D, F}, T::SMatrix{D, D, F, L}) where {D, F, L}
        new{D, F, L}(p0, T)
    end
end

function BarycentricCoordTransform(tri::Triangle)
    v = vertices(tri)
    p0 = v[1]
    T = [v[2] - p0 v[3] - p0]
    D = length(p0.coords)
    F = eltype(p0.coords)
    L = length(T)
    BarycentricCoordTransform(p0, T)
end

function tobarycoord(pt::Point{D}, bct::BarycentricCoordTransform{D}) where D
    bct.T \ (pt - bct.p0)
end

function tonatcoord(pt::Point{D}, bct::BarycentricCoordTransform{D}) where D
    bct.T * pt + bct.p0
end

function interpolation_wts(pt::Point{D}, mesh::SimpleMesh{D}, basis::PiecewiseLinearBasis{D}) where D
    el_idx = findfirst(Ref(pt) .∈ mesh)
    node_idx = mesh.topology.connec[el_idx].indices
    tri = mesh[el_idx]
    wts = _interp_wts(pt, tri)

    (cols = [node_idx...], wts = wts)
end

function _interp_wts(pt::Point, tri::Triangle)
    bct = BarycentricCoordTransform(tri)
    wt23 = tobarycoord(pt, bct)
    wt1 = 1 - sum(wt23)
    [wt1; wt23]
end

function interpolation_operator(pts::Vector{<:Point}, mesh::SimpleMesh, basis::PiecewiseLinearBasis)
    Is = Vector{Int}()
    Js = Vector{Int}()
    Vs = Vector{Float64}()
    js = zeros(3)
    vs = zeros(3)
    for (i, pt) in enumerate(pts)
        js[:], vs[:] = interpolation_wts(pt, mesh, basis)
        perm = sortperm(js)
        append!(Is, repeat([i], 3))
        append!(Js, js[perm])
        append!(Vs, vs[perm])
    end
    dropzeros!(sparse(Is, Js, Vs))
end

function interpolate(pt::Point{2}, md::MeshData, sym, basis::PiecewiseLinearBasis{2})
    js, vs = interpolation_wts(pt, domain(md), basis)
    vals = values(md, 0)[sym]
    sum(vs .* vals[js])
end

end


text: Dict{String,Any}("black" => Dict{String,Any}("allow_prereleases" => Dict{String,Any}("value" => "true","type" => "bool"),"python" => Dict{String,Any}("value" => ">3.6","type" => "string"),"version" => Dict{String,Any}("value" => ">=18.9b0","type" => "string")))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
plan_flows_fft(a::Array, args...; kwargs...) = plan_fft(a, args...; kwargs...)
plan_flows_rfft(a::Array, args...; kwargs...) = plan_rfft(a, args...; kwargs...)

"""
A placeholder grid object for `0D` problems (in other words, systems of ODEs).
"""
struct ZeroDGrid{T, A, Alias} <: AbstractGrid{T, A, Alias} end

"""
    struct OneDGrid{T<:AbstractFloat, Tk, Tx, Tfft, Trfft, Talias} <: AbstractGrid{T, Tk, Talias}

A one-dimensional `grid`.

$(TYPEDFIELDS)
"""
struct OneDGrid{T<:AbstractFloat, Tk, Tx, Tfft, Trfft, Talias} <: AbstractGrid{T, Tk, Talias}
    "number of points in x"
                nx :: Int
    "number of wavenumbers in x"
                nk :: Int
    "number of positive wavenumbers in x (real Fourier transforms)"
               nkr :: Int
    "grid spacing in x"
                dx :: T
    "domain extent in x"
                Lx :: T
    "range with x-grid-points"
                 x :: Tx
    "array with x-wavenumbers"
                 k :: Tk
    "array with positive x-wavenumbers (real Fourier transforms)"
                kr :: Tk
    "array with inverse squared k-wavenumbers, 1/k²"
            invksq :: Tk
    "array with inverse squared kr-wavenumbers, 1/kr²"
           invkrsq :: Tk
    "the FFT plan for complex-valued fields"
           fftplan :: Tfft
    "the FFT plan for real-valued fields"
          rfftplan :: Trfft
    "the fraction of wavenumbers that are aliased (e.g., 1/3 for quadradic nonlinearities)"
  aliased_fraction :: T
    "range of the indices of aliased x-wavenumbers"
            kalias :: Talias
    "range of the indices of aliased positive x-wavenumbers (real Fourier transforms)"
           kralias :: Talias
end

"""
    OneDGrid(nx, Lx; x0=-Lx/2, nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE,
                      T=Float64, aliased_fraction=1/3, ArrayType=Array)

Constructs a OneDGrid object with size `Lx`, resolution `nx`, and leftmost position `x0`.
FFT plans are generated for `nthreads` CPUs using FFTW flag `effort`. The float type is `T`
and the array types is `ArrayType`. The `aliased_fraction` keyword determines the highest
wavenubers that are being zero-ed out by `dealias!` function; 1/3 is the nominal value for
quadratic nonlinearities.
"""
function OneDGrid(nx, Lx; x0=-Lx/2, nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE,
                  T=Float64, aliased_fraction=1/3, ArrayType=Array)

  dx = Lx/nx

   nk = nx
  nkr = Int(nx/2 + 1)

  # Physical grid
  x = range(T(x0), step=T(dx), length=nx)

  # Wavenubmer grid
   k = ArrayType{T}( fftfreq(nx, 2π/Lx*nx))
  kr = ArrayType{T}(rfftfreq(nx, 2π/Lx*nx))

   invksq = @. 1 / k^2
  invkrsq = @. 1 / kr^2
  CUDA.@allowscalar  invksq[1] = 0
  CUDA.@allowscalar invkrsq[1] = 0

  FFTW.set_num_threads(nthreads)
   fftplan = plan_flows_fft(ArrayType{Complex{T}, 1}(undef, nx), flags=effort)
  rfftplan = plan_flows_rfft(ArrayType{T, 1}(undef, nx), flags=effort)

  kalias, kralias = getaliasedwavenumbers(nk, nkr, aliased_fraction)

  Tx = typeof(x)
  Tk = typeof(k)
  Tfft = typeof(fftplan)
  Trfft = typeof(rfftplan)
  Talias = typeof(kalias)

  return OneDGrid{T, Tk, Tx, Tfft, Trfft, Talias}(nx, nk, nkr, dx, Lx, x, k, kr,
                                      invksq, invkrsq, fftplan, rfftplan,
                                      aliased_fraction, kalias, kralias)
end

"""
    struct TwoDGrid{T<:AbstractFloat, Tk, Tx, Tfft, Trfft, Talias} <: AbstractGrid{T, Tk, Talias}

A two-dimensional `grid`.

$(TYPEDFIELDS)
"""
struct TwoDGrid{T<:AbstractFloat, Tk, Tx, Tfft, Trfft, Talias} <: AbstractGrid{T, Tk, Talias}
    "number of points in x"
               nx :: Int
    "number of points in y"
               ny :: Int
    "number of wavenumbers in x"
               nk :: Int
    "number of wavenumbers in y"
               nl :: Int
    "number of positive wavenumers in x (real Fourier transforms)"
               nkr :: Int
    "grid spacing in x"
                dx :: T
    "grid spacing in y"
                dy :: T
    "domain extent in x"
                Lx :: T
    "domain extent in y"
                Ly :: T
    "range with x-grid-points"
                 x :: Tx
    "range with y-grid-points"
                 y :: Tx
    "array with x-wavenumbers"
                 k :: Tk
    "array with y-wavenumbers"
                 l :: Tk
    "array with positive x-wavenumbers (real Fourier transforms)"
                kr :: Tk
    "array with squared total wavenumbers, k²+l²"
               Ksq :: Tk
    "array with inverse squared total wavenumbers, 1/(k²+l²)"
            invKsq :: Tk
    "array with squared total wavenumbers for real Fourier transforms, kr²+l²"
              Krsq :: Tk
    "array with inverse squared total wavenumbers for real Fourier transforms, 1/(kr²+l²)"
           invKrsq :: Tk
    "the FFT plan for complex-valued fields"
           fftplan :: Tfft
    "the FFT plan for real-valued fields"
          rfftplan :: Trfft
    "the fraction of wavenumbers that are aliased (e.g., 1/3 for quadradic nonlinearities)"
  aliased_fraction :: T
    "range of the indices of aliased x-wavenumbers"
            kalias :: Talias
    "range of the indices of aliased positive x-wavenumbers (real Fourier transforms)"
           kralias :: Talias
    "range of the indices of aliased y-wavenumbers"
            lalias :: Talias
end

"""
    TwoDGrid(nx, Lx, ny=nx, Ly=Lx; x0=-Lx/2, y0=-Ly/2, nthreads=Sys.CPU_THREADS,
                      effort=FFTW.MEASURE, T=Float64, aliased_fraction=1/3, ArrayType=Array)

Constructs a TwoDGrid object with size `Lx`, `Ly`, resolution `nx`, `ny`, and leftmost
positions `x0`, `y0`. FFT plans are generated for `nthreads` CPUs using FFTW flag `effort`.
The float type is `T` and the array types is `ArrayType`. The `aliased_fraction` keyword
determines the highest wavenubers that are being zero-ed out by `dealias!` function; 1/3 is
the nominal value for quadratic nonlinearities.
"""
function TwoDGrid(nx, Lx, ny=nx, Ly=Lx; x0=-Lx/2, y0=-Ly/2, nthreads=Sys.CPU_THREADS,
                  effort=FFTW.MEASURE, T=Float64, aliased_fraction=1/3, ArrayType=Array)

  dx = Lx/nx
  dy = Ly/ny

   nk = nx
   nl = ny
  nkr = Int(nx/2 + 1)

  # Physical grid
  x = range(T(x0), step=T(dx), length=nx)
  y = range(T(y0), step=T(dy), length=ny)

  # Wavenubmer grid
   k = ArrayType{T}(reshape( fftfreq(nx, 2π/Lx*nx), (nk, 1)))
   l = ArrayType{T}(reshape( fftfreq(ny, 2π/Ly*ny), (1, nl)))
  kr = ArrayType{T}(reshape(rfftfreq(nx, 2π/Lx*nx), (nkr, 1)))

     Ksq = @. k^2 + l^2
  invKsq = @. 1 / Ksq
  CUDA.@allowscalar invKsq[1, 1] = 0

     Krsq = @. kr^2 + l^2
  invKrsq = @. 1 / Krsq
  CUDA.@allowscalar invKrsq[1, 1] = 0

  # FFT plans
  FFTW.set_num_threads(nthreads)
  fftplan = plan_flows_fft(ArrayType{Complex{T}, 2}(undef, nx, ny), flags=effort)
  rfftplan = plan_flows_rfft(ArrayType{T, 2}(undef, nx, ny), flags=effort)

  kalias, kralias = getaliasedwavenumbers(nk, nkr, aliased_fraction)
  lalias, _ = getaliasedwavenumbers(nl, nl, aliased_fraction)

  Tx = typeof(x)
  Tk = typeof(k)
  Tfft = typeof(fftplan)
  Trfft = typeof(rfftplan)
  Talias = typeof(kalias)

  return TwoDGrid{T, Tk, Tx, Tfft, Trfft, Talias}(nx, ny, nk, nl, nkr, dx, dy, Lx, Ly, x, y, k, l, kr,
                                          Ksq, invKsq, Krsq, invKrsq, fftplan, rfftplan,
                                          aliased_fraction, kalias, kralias, lalias)
end

"""
    struct ThreeDGrid{T<:AbstractFloat, Tk, Tx, Tfft, Trfft, Talias} <: AbstractGrid{T, Tk, Talias}

A three-dimensional `grid`.

$(TYPEDFIELDS)
"""
struct ThreeDGrid{T<:AbstractFloat, Tk, Tx, Tfft, Trfft, Talias} <: AbstractGrid{T, Tk, Talias}
    "number of points in x"
               nx :: Int
    "number of points in y"
               ny :: Int
    "number of points in z"
               nz :: Int
    "number of wavenumbers in x"
               nk :: Int
    "number of wavenumbers in y"
               nl :: Int
    "number of wavenumbers in z"
               nm :: Int
    "number of positive wavenumers in x (real Fourier transforms)"
               nkr :: Int
    "grid spacing in x"
                dx :: T
    "grid spacing in y"
                dy :: T
    "grid spacing in z"
                dz :: T
    "domain extent in x"
                Lx :: T
    "domain extent in y"
                Ly :: T
    "domain extent in z"
                Lz :: T
    "range with x-grid-points"
                 x :: Tx
    "range with y-grid-points"
                 y :: Tx
    "range with z-grid-points"
                 z :: Tx
    "array with x-wavenumbers"
                 k :: Tk
    "array with y-wavenumbers"
                 l :: Tk
    "array with z-wavenumbers"
                 m :: Tk
    "array with positive x-wavenumbers (real Fourier transforms)"
                kr :: Tk
    "array with squared total wavenumbers, k²+l²+m²"
               Ksq :: Tk
    "array with inverse squared total wavenumbers, 1/(k²+l²+m²)"
            invKsq :: Tk
    "array with squared total wavenumbers for real Fourier transforms, kr²+l²+m²"
              Krsq :: Tk
    "array with inverse squared total wavenumbers for real Fourier transforms, 1/(kr²+l²+m²)"
           invKrsq :: Tk
    "the FFT plan for complex-valued fields"
           fftplan :: Tfft
    "the FFT plan for real-valued fields"
          rfftplan :: Trfft
    "the fraction of wavenumbers that are aliased (e.g., 1/3 for quadradic nonlinearities)"
  aliased_fraction :: T
    "range of the indices of aliased x-wavenumbers"
            kalias :: Talias
    "range of the indices of aliased positive x-wavenumbers (real Fourier transforms)"
           kralias :: Talias
    "range of the indices of aliased y-wavenumbers"
            lalias :: Talias
    "range of the indices of aliased y-wavenumbers"
            malias :: Talias
end

"""
    ThreeDGrid(nx, Lx, ny=nx, Ly=Lx, nz=nx, Lz=Lx; x0=-Lx/2, y0=-Ly/2, z0=-Lz/2,
    nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE, T=Float64, aliased_fraction=1/3, ArrayType=Array)

 Constructs a TwoDGrid object with size `Lx`, `Ly`, `Lz`, resolution `nx`, `ny`, `nz` and
 leftmost positions `x0`, `y0`, `z0`. FFT plans are generated for `nthreads` CPUs using FFTW
 flag `effort`. The float type is `T` and the array types is `ArrayType`. The `aliased_fraction`
 keyword determines the highest wavenubers that are being zero-ed out by `dealias!` function;
 1/3 is the nominal value for quadratic nonlinearities.
"""
function ThreeDGrid(nx, Lx, ny=nx, Ly=Lx, nz=nx, Lz=Lx; x0=-Lx/2, y0=-Ly/2, z0=-Lz/2,
                  nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE, T=Float64, aliased_fraction=1/3, ArrayType=Array)

  dx = Lx/nx
  dy = Ly/ny
  dz = Lz/nz

  nk = nx
  nl = ny
  nm = nz
  nkr = Int(nx/2 + 1)

  # Physical grid
  x = range(T(x0), step=T(dx), length=nx)
  y = range(T(y0), step=T(dy), length=ny)
  z = range(T(z0), step=T(dz), length=nz)

  # Wavenubmer grid
   k = ArrayType{T}(reshape( fftfreq(nx, 2π/Lx*nx), (nk, 1, 1)))
   l = ArrayType{T}(reshape( fftfreq(ny, 2π/Ly*ny), (1, nl, 1)))
   m = ArrayType{T}(reshape( fftfreq(nz, 2π/Lz*nz), (1, 1, nm)))
  kr = ArrayType{T}(reshape(rfftfreq(nx, 2π/Lx*nx), (nkr, 1, 1)))

     Ksq = @. k^2 + l^2 + m^2
  invKsq = @. 1 / Ksq
  CUDA.@allowscalar invKsq[1, 1, 1] = 0

     Krsq = @. kr^2 + l^2 + m^2
  invKrsq = @. 1 / Krsq
  CUDA.@allowscalar invKrsq[1, 1, 1] = 0

  # FFT plans
  FFTW.set_num_threads(nthreads)
  fftplan = plan_flows_fft(ArrayType{Complex{T}, 3}(undef, nx, ny, nz), flags=effort)
  rfftplan = plan_flows_rfft(ArrayType{T, 3}(undef, nx, ny, nz), flags=effort)

  kalias, kralias = getaliasedwavenumbers(nk, nkr, aliased_fraction)
  lalias, _ = getaliasedwavenumbers(nl, Int(nl/2+1), aliased_fraction)
  malias, _ = getaliasedwavenumbers(nm, Int(nm/2+1), aliased_fraction)

  Tx = typeof(x)
  Tk = typeof(k)
  Tfft = typeof(fftplan)
  Trfft = typeof(rfftplan)
  Talias = typeof(kalias)

  return ThreeDGrid{T, Tk, Tx, Tfft, Trfft, Talias}(nx, ny, nz, nk, nl, nm, nkr,
                                        dx, dy, dz, Lx, Ly, Lz, x, y, z, k, l, m, kr,
                                        Ksq, invKsq, Krsq, invKrsq, fftplan, rfftplan,
                                        aliased_fraction, kalias, kralias, lalias, malias)
end

Base.eltype(grid::OneDGrid) = eltype(grid.x)
Base.eltype(grid::TwoDGrid) = eltype(grid.x)
Base.eltype(grid::ThreeDGrid) = eltype(grid.x)

OneDGrid(dev::CPU, args...; kwargs...) = OneDGrid(args...; ArrayType=Array, kwargs...)
TwoDGrid(dev::CPU, args...; kwargs...) = TwoDGrid(args...; ArrayType=Array, kwargs...)
ThreeDGrid(dev::CPU, args...; kwargs...) = ThreeDGrid(args...; ArrayType=Array, kwargs...)

"""
    gridpoints(grid::TwoDGrid)
    gridpoints(grid::ThreeDGrid)

Returns the collocation points of the `grid` in 2D or 3D arrays `X, Y` (and `Z`).
"""
function gridpoints(grid::TwoDGrid{T, A}) where {T, A}
  X = [ grid.x[i₁] for i₁=1:grid.nx, i₂=1:grid.ny ]
  Y = [ grid.y[i₂] for i₁=1:grid.nx, i₂=1:grid.ny ]

  return A(X), A(Y)
end

function gridpoints(grid::ThreeDGrid{T, A}) where {T, A}
  X = [ grid.x[i₁] for i₁=1:grid.nx, i₂=1:grid.ny, i₃=1:grid.nz ]
  Y = [ grid.y[i₂] for i₁=1:grid.nx, i₂=1:grid.ny, i₃=1:grid.nz ]
  Z = [ grid.z[i₃] for i₁=1:grid.nx, i₂=1:grid.ny, i₃=1:grid.nz ]

  return A(X), A(Y), A(Z)
end

"""
    getaliasedwavenumbers(nk, nkr, aliased_fraction)

Returns the top `aliased_fraction` highest wavenumbers, both for and real FFTs, `kalias` and
`kralias` respectively. For example, `aliased_fraction=1/3` should return the indices of the
top-most 1/6-th (in absolute value) for both positive and negative wavenumbers (i.e., 1/3 total)
that should be set to zero after performing an FFT.
"""
function getaliasedwavenumbers(nk, nkr, aliased_fraction)
  L = (1 - aliased_fraction)/2 # e.g., (1 - 1/3) / 2 + 1 = 1/3.
  R = (1 + aliased_fraction)/2 # e.g., (1 + 1/3) / 2 - 1 = 2/3.

  iL = floor(Int, L * nk) + 1
  iR =  ceil(Int, R * nk)

  aliased_fraction < 1 || error("`aliased_fraction` must be less than 1")

   kalias = (aliased_fraction > 0) ? (iL:iR) : nothing
  kralias = (aliased_fraction > 0) ? (iL:nkr) : nothing

  return kalias, kralias
end

"""
    dealias!(fh, grid)

Dealias array `fh` on the `grid` based on the `grids`'s `aliased_factor`.
"""
function dealias!(fh, grid)
   _dealias!(fh, grid)

   return nothing
end

dealias!(fh, grid::AbstractGrid{T, A, Nothing}) where {T, A} = nothing

function _dealias!(fh, grid::OneDGrid)
  kalias = size(fh, 1) == grid.nkr ? grid.kralias : grid.kalias
  _dealias!(fh, grid, kalias)

  return nothing
end

function _dealias!(fh, grid::OneDGrid, kalias)
  @views @. fh[kalias, :] = 0

  return nothing
end

function _dealias!(fh, grid::TwoDGrid)
  kalias = size(fh, 1) == grid.nkr ? grid.kralias : grid.kalias
  _dealias!(fh, grid, kalias)

  return nothing
end

function _dealias!(fh, grid::TwoDGrid, kalias)
  @views @. fh[kalias, :, :] = 0
  @views @. fh[:, grid.lalias, :] = 0

  return nothing
end

function _dealias!(fh, grid::ThreeDGrid)
  kalias = size(fh, 1) == grid.nkr ? grid.kralias : grid.kalias
  _dealias!(fh, grid, kalias)

  return nothing
end

function _dealias!(fh, grid::ThreeDGrid, kalias)
  @views @. fh[kalias, :, :, :] = 0
  @views @. fh[:, grid.lalias, :, :] = 0
  @views @. fh[:, :, grid.malias, :] = 0

  return nothing
end

"""
    makefilter(K; order=4, innerK=0.65, outerK=1)

Returns a filter acting on the non-dimensional wavenumber `K` that decays exponentially
for `K`>`innerK`, thus removing high-wavenumber content from a spectrum it is multiplied
with. The decay rate is determined by order and `outerK` determines the outer wavenumber
at which the filter is smaller than Float64 machine precision.
"""
function makefilter(K::Array; order=4, innerK=0.65, outerK=1)
  TK = typeof(K)
  K = Array(K)
  decay = 15*log(10) / (outerK - innerK)^order # decay rate for filtering function
  filter = @. exp(-decay*(K - innerK)^order)
  filter[K .< innerK] .= 1

  return TK(filter)
end

function makefilter(g::OneDGrid; realvars=true, kwargs...)
  K = realvars ? g.kr*g.dx/π : @.(abs(g.k*g.dx/π))

  return makefilter(K; kwargs...)
end

function makefilter(g::TwoDGrid; realvars=true, kwargs...)
  K = realvars ?
      @.(sqrt((g.kr*g.dx/π)^2 + (g.l*g.dy/π)^2)) : @.(sqrt((g.k*g.dx/π)^2 + (g.l*g.dy/π)^2))

  return makefilter(K; kwargs...)
end

function makefilter(g::ThreeDGrid; realvars=true, kwargs...)
  K = realvars ?
      @.(sqrt((g.kr*g.dx/π)^2 + (g.l*g.dy/π)^2 + (g.m*g.dz/π)^2)) : @.(sqrt((g.k*g.dx/π)^2 + (g.l*g.dy/π)^2 + (g.m*g.dz/π)^2))

  return makefilter(K; kwargs...)
end

makefilter(g, T, sz; kwargs...) = ones(T, sz) .* makefilter(g; realvars=sz[1]==g.nkr, kwargs...)
makefilter(eq; kwargs...) = makefilter(eq.grid, fltype(eq.T), eq.dims; kwargs...)

"""
    griddevice(grid)

Returns the device on which the `grid` lives on.
"""
griddevice(grid::AbstractGrid{T, A}) where {T, A} = A<:Array ? "CPU" : "GPU"

show(io::IO, g::OneDGrid{T}) where T =
     print(io, "OneDimensionalGrid\n",
               "  ├─────────── Device: ", griddevice(g), '\n',
               "  ├──────── FloatType: $T", '\n',
               "  ├────────── size Lx: ", g.Lx, '\n',
               "  ├──── resolution nx: ", g.nx, '\n',
               "  ├── grid spacing dx: ", g.dx, '\n',
               "  ├─────────── domain: x ∈ [$(g.x[1]), $(g.x[end])]", '\n',
               "  └─ aliased fraction: ", g.aliased_fraction)

show(io::IO, g::TwoDGrid{T}) where T =
     print(io, "TwoDimensionalGrid\n",
               "  ├───────────────── Device: ", griddevice(g), '\n',
               "  ├────────────── FloatType: $T", '\n',
               "  ├────────── size (Lx, Ly): ", (g.Lx, g.Ly), '\n',
               "  ├──── resolution (nx, ny): ", (g.nx, g.ny), '\n',
               "  ├── grid spacing (dx, dy): ", (g.dx, g.dy), '\n',
               "  ├───────────────── domain: x ∈ [$(g.x[1]), $(g.x[end])]", '\n',
               "  |                          y ∈ [$(g.y[1]), $(g.y[end])]", '\n',
               "  └─ aliased fraction: ", g.aliased_fraction)

show(io::IO, g::ThreeDGrid{T}) where T =
     print(io, "ThreeDimensionalGrid\n",
               "  ├───────────────────── Device: ", griddevice(g), '\n',
               "  ├────────────────── FloatType: $T", '\n',
               "  ├────────── size (Lx, Ly, Lz): ", (g.Lx, g.Ly, g.Ly), '\n',
               "  ├──── resolution (nx, ny, nz): ", (g.nx, g.ny, g.nz), '\n',
               "  ├── grid spacing (dx, dy, dz): ", (g.dx, g.dy, g.dz), '\n',
               "  ├────────────────────  domain: x ∈ [$(g.x[1]), $(g.x[end])]", '\n',
               "  |                              y ∈ [$(g.y[1]), $(g.y[end])]", '\n',
               "  |                              z ∈ [$(g.z[1]), $(g.z[end])]", '\n',
               "  └─ aliased fraction: ", g.aliased_fraction)


text:
mutable struct DelayRejectionJumpAggregation{T,S,F1,F2,RNG} <: AbstractDSSAJumpAggregator
    next_jump::Int
    prev_jump::Int
    next_jump_time::T
    end_time::T
    cur_rates::Vector{T} # cur_rates here is the cumsum of rates
    sum_rate::T
    ma_jumps::S
    rates::F1
    affects!::F2
    save_positions::Tuple{Bool,Bool}
    rng::RNG
    next_delay::Union{Nothing,Vector{Int}}
    num_next_delay::Union{Nothing,Vector{Int}}
    time_to_next_jump::T
    dt_delay::T
end
#3
function DelayRejectionJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, maj::S, rs::F1, affs!::F2, sps::Tuple{Bool,Bool}, rng::RNG; kwargs...) where {T,S,F1,F2,RNG}
    nd = nothing
    nnd = nothing
    ttnj = zero(et)
    dt_delay = zero(et)
    #4
    DelayRejectionJumpAggregation{T,S,F1,F2,RNG}(nj, nj, njt, et, crs, sr, maj, rs, affs!, sps, rng, nd, nnd, ttnj, dt_delay)
end
#1
function aggregate(aggregator::DelayRejection, u, p, t, end_time, constant_jumps,
    ma_jumps, save_positions, rng; kwargs...)

    # handle constant jumps using function wrappers
    rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps)
    build_jump_aggregation(DelayRejectionJumpAggregation, u, p, t, end_time, ma_jumps, rates, affects!, save_positions, rng; kwargs...)
end

function initialize!(p::DelayRejectionJumpAggregation, integrator, u, params, t)
    find_next_delay_dt!(p, integrator)
    generate_jumps!(p, integrator, u, params, t)
    nothing
end

@inline function execute_jumps!(p::DelayRejectionJumpAggregation, integrator, u, params, t)
    # update_state!(p, integrator, u)
    update_state_delay!(p, integrator, u, t)
    nothing
end

function generate_jumps!(p::DelayRejectionJumpAggregation, integrator, u, params, t)
    time_to_next_jump!(p, integrator, u, params, t)
    if p.next_delay == nothing
        @inbounds p.next_jump = searchsortedfirst(p.cur_rates, rand(p.rng) * p.sum_rate) # cur_rates is the cumsum of rates
    else
        p.next_jump = 0
    end
    nothing
end

@fastmath function time_to_next_jump!(p::DelayRejectionJumpAggregation, integrator, u, params, t)
    prev_rate = zero(t)
    new_rate = zero(t)
    cur_rates = p.cur_rates
    # mass action rates
    majumps = p.ma_jumps
    idx = get_num_majumps(majumps)
    @inbounds for i in 1:idx
        new_rate = evalrxrate(u, i, majumps)
        cur_rates[i] = new_rate + prev_rate
        prev_rate = cur_rates[i]
    end
    # constant jump rates
    idx += 1
    rates = p.rates
    @inbounds for i in 1:length(p.rates)
        new_rate = rates[i](u, params, t)
        cur_rates[idx] = new_rate + prev_rate
        prev_rate = cur_rates[idx]
        idx += 1
    end

    @inbounds sum_rate = cur_rates[end]
    dt_reaction = randexp(p.rng) / sum_rate
    # consider delay dt
    dt_delay_generation!(p, integrator)
    compare_delay!(p, integrator.de_chan, p.dt_delay, dt_reaction, t)
    p.sum_rate = sum_rate
    nothing
end


text:
## Define the model
## parameters
γ₁, β₁, λ₁, θ₁, κ = 11.0, 1.0, 12.48, 1.6, 2.6
γ₂, β₂, λ₂, θ₂, xstar = 22.0, 1.0, 16.25, 1.0, 1.6

## state space is {11,10,01,00}
T = [
    -(γ₁ + γ₂) γ₂ γ₁ 0;
    β₂ -(γ₁ + β₂) 0 γ₁;
    β₁ 0 -(γ₂ + β₁) γ₂;
    0 β₁ β₂ -(β₂ + β₂);
    ]

C = [λ₁ - θ₁; λ₁ - θ₁; -θ₁; -θ₁]

r₁₁(x) = (λ₂ - κ) * (x .== 0) +
    (λ₂ - θ₂) * ((x .> 0) .& (x .<= xstar)) +
    λ₂ * (x .> xstar)
r₁₀(x) = -κ * (x .== 0) +
    -θ₂ * ((x .> 0) .& (x .<= xstar))
r₀₁(x) = (λ₂ - κ) * (x .== 0) +
    (λ₂ - θ₂) * ((x .> 0) .& (x .<= xstar)) +
    λ₂ * (x .> xstar)
r₀₀(x) = -κ * (x .== 0) +
    -θ₂ * ((x .> 0) .& (x .<= xstar))

R₁₁(x) = x .* (λ₂ - θ₂) .* ((x .> 0) .& (x .< xstar)) .+
    ((x .- xstar) .* λ₂ .+ xstar .* (λ₂ - θ₂)) .* (x .>= xstar)
R₁₀(x) = x .* -θ₂ * ((x .> 0) .& (x .< xstar)) + xstar .* -θ₂ * (x .>= xstar)
R₀₁(x) = x .* (λ₂ - θ₂) .* ((x .> 0) .& (x .< xstar)) .+
    ((x .- xstar) .* λ₂ .+ xstar .* (λ₂ - θ₂)) .* (x .>= xstar)
R₀₀(x) = x .* -θ₂ * ((x .> 0) .& (x .< xstar)) + xstar .* -θ₂ * (x .>= xstar)

r = (
    r = function (x)
        [r₁₁(x) r₁₀(x) r₀₁(x) r₀₀(x)]
    end,
    R = function (x)
        [R₁₁(x) R₁₀(x) R₀₁(x) R₀₀(x)]
    end,
)

## make models
simBounds = [0 Inf; -Inf Inf] # bounds for simulation only
@testset begin
    simModel1 = SFFM.Model(T, C, r, Bounds = simBounds)
    @test sum(simModel1.T,dims=2) ≈ zeros(4)
    @test simModel1.r.r([0;1]) ≈ [13.65  -2.6  13.65  -2.6; 15.25  -1.0  15.25  -1.0]
    @test simModel1.r.R(0) ≈ [0 0 0 0]
end
simModel = SFFM.Model( T, C, r, Bounds = simBounds)

approxBounds = [0 48; -Inf Inf] # bounds for approximation only
model = SFFM.Model(T, C, r, Bounds = approxBounds)

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
if Sys.islinux()
    py"""
    import tensorflow as tf
    libFwiOp = tf.load_op_library('./build/libFwiOp.so')
    @tf.custom_gradient
    def fwi_op(cp,cs,den,stf,gpu_id,shot_ids,para_fname):
        res = libFwiOp.fwi_op(cp,cs,den,stf,gpu_id,shot_ids,para_fname)
        def grad(dy):
            return libFwiOp.fwi_op_grad(dy, tf.constant(1.0,dtype=tf.float64),cp,cs,den,stf,gpu_id,shot_ids,para_fname)
        return res, grad
    def fwi_obs_op(cp,cs,den,stf,gpu_id,shot_ids,para_fname):
        res = libFwiOp.fwi_obs_op(cp,cs,den,stf,gpu_id,shot_ids,para_fname)
        return res
    """
    elseif Sys.isapple()
    py"""
    import tensorflow as tf
    libFwiOp = tf.load_op_library('./build/libFwiOp.dylib')
    @tf.custom_gradient
    def fwi_op(cp,cs,den):
        res = libFwiOp.fwi_op(cp,cs,den,stf,gpu_id,shot_ids,para_fname)
        def grad(dy):
            return libFwiOp.fwi_op_grad(dy,tf.constant(1.0,dtype=tf.float64),cp,cs,den,stf,gpu_id,shot_ids,para_fname)
        return res, grad
    def fwi_obs_op(cp,cs,den,stf,gpu_id,shot_ids,para_fname):
        res = libFwiOp.fwi_obs_op(cp,cs,den,stf,gpu_id,shot_ids,para_fname)
        return res
    """
    elseif Sys.iswindows()
    py"""
    import tensorflow as tf
    libFwiOp = tf.load_op_library('./build/libFwiOp.dll')
    @tf.custom_gradient
    def fwi_op(cp,cs,den):
        res = libFwiOp.fwi_op(cp,cs,den,stf,gpu_id,shot_ids,para_fname)
        def grad(dy):
            return libFwiOp.fwi_op_grad(dy,tf.constant(1.0,dtype=tf.float64),cp,cs,den,stf,gpu_id,shot_ids,para_fname)
        return res, grad
    def fwi_obs_op(cp,cs,den,stf,gpu_id,shot_ids,para_fname):
        res = libFwiOp.fwi_obs_op(cp,cs,den,stf,gpu_id,shot_ids,para_fname)
        return res
    """
    end

    fwi_op = py"fwi_op"
    fwi_obs_op = py"fwi_obs_op"


text:
function track_src!(S::GphysData, j::Int64, nx::Array{Int64,1}, last_src::Array{Int64,1})
  n = length(nx)

  # Check existing channels for changes
  for i in 1:n
    if length(S.x[i]) > nx[i]
      last_src[i] = j
      nx[i] = length(S.x[i])
    end
  end

  # Add new channels
  if n < S.n
    δn = S.n - n
    append!(nx, zeros(Int64, δn))
    append!(last_src, zeros(Int64, δn))
    for i in n+1:S.n
      nx[i] = length(S.x[i])
      last_src[i] = j
    end
  end
  return nothing
end


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

name = "cilantro"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/kzampog/cilantro.git", "d7e654be0deb262aa50773793a0080652c79e473")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cilantro/

mkdir build && cd build/

cmake .. \
-DCMAKE_INSTALL_PREFIX=${prefix} \
-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
-DCMAKE_BUILD_TYPE=Release \
-DBUILD_SHARED_LIBS=ON \
-DBUILD_EXAMPLES=OFF \
-DENABLE_NATIVE_BUILD_OPTIMIZATIONS=OFF \
-DENABLE_NON_DETERMINISTIC_PARALLELISM=OFF

make -j${nproc}
make install

if [[ "${target}" == *-mingw* ]]; then
#cmake install only grabs the .dll.a and leaves the actual .dll behind, manually move it
mv *.dll ${libdir}
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line

 platforms = expand_cxxstring_abis(supported_platforms(; experimental = true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libcilantro", :libcilantro)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="Eigen_jll", uuid="bc6bbf8a-a594-5541-9c57-10b0d0312c70"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
#uses CXX17 cmake standard, 5 seems to be lowest version number I can get to compile
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"5")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

using Documenter, FEMBasis

cp("../README.md", "src/index.md"; remove_destination=true)

makedocs(modules=[FEMBasis],
         format = :html,
         checkdocs = :all,
         sitename = "FEMBasis.jl",
         pages = [
                  "Introduction" => "index.md",
                  "API documentation" => "api.md",
                 ])


text:
using Revise
using JDF, CSV, DataFrames
@time a = gf()

# 2G file
# @time a = CSV.read(
#     "C:/data/Performance_All/Performance_2010Q3.txt",
#     delim = '|',
#     header = false
# );

b = type_compress(Array(a[!, :ORIGIN]));

io = iow()
@time metadata = compress_then_write(b, io)
close(io)

using Revise
using JDF, CSV, DataFrames
# using JLSO
# metadata = JLSO.load("C:/data/metatmp")["data"]

io = ior()
@time oo = column_loader(metadata.type, io, metadata);
close(io)

all(b.==oo)

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
module MyQuantum
using LinearAlgebra

include("base.jl")
include("test.jl")
# greet() = print("Hello World!")

export is_unitary, is_hermitian, is_normed, conj_transpose, outer
export Id, X, Y, Z, H, S, T, Rn, Rx, Ry, Rz, Cgate, bell_states
export ket_0, ket_1, ket_plus, ket_minus, commutator

export test_Cgates

end # module

max_stars_repo_name: anandijain/MyQuantum

max_issues_repo_name: anandijain/MyQuantum

max_forks_repo_name: anandijain/MyQuantum


text:
params = SimulatedGaussMeterParams()
gauss = SimulatedGaussMeter(deviceID="simulated_gaussmeter", params=params, dependencies=Dict{String, Union{Device, Missing}}())

@test gauss isa SimulatedGaussMeter
@test getXValue(gauss) == 1.0u"mT"
@test getYValue(gauss) == 2.0u"mT"
@test getZValue(gauss) == 3.0u"mT"
@test getXYZValues(gauss) == [1.0, 2.0, 3.0]u"mT"