text:
#=

  Euler #13 in Julia.

  Problem 13
  """
  Work out the first ten digits of the sum of the following
  one-hundred 50-digit numbers.
    37107287533902102798797998220837590246510135740250
    ....
    20849603980134001723930671666823555245252804609722
    53503534226472524250874054075591789781264330331690")
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

ns =[37107287533902102798797998220837590246510135740250,
     46376937677490009712648124896970078050417018260538,
     74324986199524741059474233309513058123726617309629,
     91942213363574161572522430563301811072406154908250,
     23067588207539346171171980310421047513778063246676,
     89261670696623633820136378418383684178734361726757,
     28112879812849979408065481931592621691275889832738,
     44274228917432520321923589422876796487670272189318,
     47451445736001306439091167216856844588711603153276,
     70386486105843025439939619828917593665686757934951,
     62176457141856560629502157223196586755079324193331,
     64906352462741904929101432445813822663347944758178,
     92575867718337217661963751590579239728245598838407,
     58203565325359399008402633568948830189458628227828,
     80181199384826282014278194139940567587151170094390,
     35398664372827112653829987240784473053190104293586,
     86515506006295864861532075273371959191420517255829,
     71693888707715466499115593487603532921714970056938,
     54370070576826684624621495650076471787294438377604,
     53282654108756828443191190634694037855217779295145,
     36123272525000296071075082563815656710885258350721,
     45876576172410976447339110607218265236877223636045,
     17423706905851860660448207621209813287860733969412,
     81142660418086830619328460811191061556940512689692,
     51934325451728388641918047049293215058642563049483,
     62467221648435076201727918039944693004732956340691,
     15732444386908125794514089057706229429197107928209,
     55037687525678773091862540744969844508330393682126,
     18336384825330154686196124348767681297534375946515,
     80386287592878490201521685554828717201219257766954,
     78182833757993103614740356856449095527097864797581,
     16726320100436897842553539920931837441497806860984,
     48403098129077791799088218795327364475675590848030,
     87086987551392711854517078544161852424320693150332,
     59959406895756536782107074926966537676326235447210,
     69793950679652694742597709739166693763042633987085,
     41052684708299085211399427365734116182760315001271,
     65378607361501080857009149939512557028198746004375,
     35829035317434717326932123578154982629742552737307,
     94953759765105305946966067683156574377167401875275,
     88902802571733229619176668713819931811048770190271,
     25267680276078003013678680992525463401061632866526,
     36270218540497705585629946580636237993140746255962,
     24074486908231174977792365466257246923322810917141,
     91430288197103288597806669760892938638285025333403,
     34413065578016127815921815005561868836468420090470,
     23053081172816430487623791969842487255036638784583,
     11487696932154902810424020138335124462181441773470,
     63783299490636259666498587618221225225512486764533,
     67720186971698544312419572409913959008952310058822,
     95548255300263520781532296796249481641953868218774,
     76085327132285723110424803456124867697064507995236,
     37774242535411291684276865538926205024910326572967,
     23701913275725675285653248258265463092207058596522,
     29798860272258331913126375147341994889534765745501,
     18495701454879288984856827726077713721403798879715,
     38298203783031473527721580348144513491373226651381,
     34829543829199918180278916522431027392251122869539,
     40957953066405232632538044100059654939159879593635,
     29746152185502371307642255121183693803580388584903,
     41698116222072977186158236678424689157993532961922,
     62467957194401269043877107275048102390895523597457,
     23189706772547915061505504953922979530901129967519,
     86188088225875314529584099251203829009407770775672,
     11306739708304724483816533873502340845647058077308,
     82959174767140363198008187129011875491310547126581,
     97623331044818386269515456334926366572897563400500,
     42846280183517070527831839425882145521227251250327,
     55121603546981200581762165212827652751691296897789,
     32238195734329339946437501907836945765883352399886,
     75506164965184775180738168837861091527357929701337,
     62177842752192623401942399639168044983993173312731,
     32924185707147349566916674687634660915035914677504,
     99518671430235219628894890102423325116913619626622,
     73267460800591547471830798392868535206946944540724,
     76841822524674417161514036427982273348055556214818,
     97142617910342598647204516893989422179826088076852,
     87783646182799346313767754307809363333018982642090,
     10848802521674670883215120185883543223812876952786,
     71329612474782464538636993009049310363619763878039,
     62184073572399794223406235393808339651327408011116,
     66627891981488087797941876876144230030984490851411,
     60661826293682836764744779239180335110989069790714,
     85786944089552990653640447425576083659976645795096,
     66024396409905389607120198219976047599490197230297,
     64913982680032973156037120041377903785566085089252,
     16730939319872750275468906903707539413042652315011,
     94809377245048795150954100921645863754710598436791,
     78639167021187492431995700641917969777599028300699,
     15368713711936614952811305876380278410754449733078,
     40789923115535562561142322423255033685442488917353,
     44889911501440648020369068063960672322193204149535,
     41503128880339536053299340368006977710650566631954,
     81234880673210146739058568557934581403627822703280,
     82616570773948327592232845941706525094512325230608,
     22918802058777319719839450180888072429661980811197,
     77158542502016545090413245809786882778948721859617,
     72107838435069186155435662884062257473692284509516,
     20849603980134001723930671666823555245252804609722,
     53503534226472524250874054075591789781264330331690
]

# 0.01485372s
function euler13a()
    return parse(Int,string(sum(ns))[1:10])
end

run_euler(euler13a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using ModelingToolkit
using Test

@parameters t
@variables x(t) y(t) z(t)

null_op = 0*t
@test isequal(simplify(null_op), 0)

one_op = 1*t
@test isequal(simplify(one_op), t)

identity_op = Num(Term(identity,[x.val]))
@test isequal(simplify(identity_op), x)

minus_op = -x
@test isequal(simplify(minus_op), -1x)
simplify(minus_op)

@variables x

@test toexpr(expand_derivatives(Differential(x)((x-2)^2))) == :($(*)(2, $(+)(-2, x)))
@test toexpr(expand_derivatives(Differential(x)((x-2)^3))) == :($(*)(3, $(^)($(+)(-2, x), 2)))
@test toexpr(simplify(x+2+3)) == :($(+)(5, x))

d1 = Differential(x)((-2 + x)^2)
d2 = Differential(x)(d1)
d3 = Differential(x)(d2)

@test toexpr(expand_derivatives(d3)) == :(0)
@test toexpr(simplify(x^0)) == :(1)

@test ModelingToolkit.substitute(2x + y == 1, Dict(x => 0.0, y => 0.0)) === false
@test ModelingToolkit.substitute(2x + y == 1, Dict(x => 0.0, y => 1.0)) === true

ext: jl


text:
# Hessian Eigenmaps (HLLE)
# ---------------------------
# Hessian eigenmaps: Locally linear embedding techniques for high-dimensional data,
# D. Donoho and C. Grimes, Proc Natl Acad Sci U S A. 2003 May 13; 100(10): 5591–5596

#### HLLE type
immutable HLLE{T <: Real} <: SpectralResult
    k::Int
    λ::AbstractVector{T}
    proj::Projection{T}

    HLLE{T}(k::Int, λ::AbstractVector{T}, proj::Projection{T}) = new(k, λ, proj)
end

## properties
outdim(M::HLLE) = size(M.proj, 1)
projection(M::HLLE) = M.proj

eigvals(M::HLLE) = M.λ
neighbors(M::HLLE) = M.k

## show & dump
function show(io::IO, M::HLLE)
    print(io, "Hessian Eigenmaps(outdim = $(outdim(M)), neighbors = $(neighbors(M)))")
end

function dump(io::IO, M::HLLE)
    show(io, M)
    println(io, "eigenvalues: ")
    Base.showarray(io, M.λ', header=false, repr=false)
    println(io)
    println(io, "projection:")
    Base.showarray(io, M.proj, header=false, repr=false)
end

## interface functions
function transform{T<:Real}(::Type{HLLE}, X::DenseMatrix{T}; d::Int=2, k::Int=12)
    n = size(X, 2)

    # Identify neighbors
    D, I = find_nn(X, k)

    # Obtain tangent coordinates and develop Hessian estimator
    hs = Int(d*(d+1)/2)
    W = spzeros(hs*n,n)
    for i=1:n
        # re-center points in neighborhood
        μ = mean(X[:,I[:,i]],2)
        N = X[:,I[:,i]] .- μ
        # calculate tangent coordinates
        #tc = svdfact(N')[:U][:,1:d]
        tc = svdfact(N)[:V][:,1:d]

        # Develop Hessian estimator
        Yi = [ones(k) tc zeros(k,hs)]
        for ii=1:d
            Yi[:,d+ii+1] = tc[:,ii].^2
        end
        yi = 2(1+d)
        for (ii,jj) in combinations(1:d,2)
            Yi[:, yi] = tc[:, ii] .* tc[:, jj]
            yi += 1
        end
        F = qrfact(Yi)
        H = full(F[:Q])[:,d+2:end]'
        W[(i-1)*hs+(1:hs),I[:,i]] = H
    end

    # decomposition
    λ, V = decompose(W'*W, d)
    return HLLE{T}(k, λ, V' .* sqrt(n))
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

default_normalization(::Type{Continuous}) =
  x -> (quantile(x, 0.75) - quantile(x, 0.25))
default_normalization(::Type{<:Compositional}) =
  x -> maximum(norm.(x))

# fallback to no normalization
default_normalization(::Type) = nothing

function normalize(tables...)
  rtables   = Tables.rowtable.(tables)
  longtable = Tables.columntable(reduce(vcat, rtables))
  colnames  = Tables.columnnames(longtable)
  scitypes  = schema(longtable).scitypes

  constants  = map(zip(colnames, scitypes)) do (c, s)
    x = Tables.getcolumn(longtable, c)
    k = default_normalization(s)
    isnothing(k) ? k : k(x)
  end

  colconst   = zip(colnames, constants)

  map(tables) do table
    # perform normalization
    colvalues = map(colconst) do (c, k)
      x = Tables.getcolumn(table, c)
      isnothing(k) ? x : x ./ k
    end

    # return same table type
    ctor = Tables.materializer(table)
    ctor((; zip(colnames, colvalues)...))
  end
end


text:
module TestContainers
using Base.Test
using DataArrays

dv = @data ones(3)
push!(dv, 3.0)
push!(dv, NA)

@assert isequal(dv, (@data [1.0, 1.0, 1.0, 3.0, NA]))

a, b = pop!(dv), pop!(dv)
@assert isna(a)
@assert b == 3.0

unshift!(dv, 3.0)
unshift!(dv, NA)

@assert isequal(dv, (@data [NA, 3.0, 1.0, 1.0, 1.0]))

a, b = shift!(dv), shift!(dv)
@assert isna(a)
@assert b == 3.0

## SPLICE
function test_splice(dv, spliceout, splicein...)
    dv = copy(dv)
    v = Any[x for x in dv]
    expval = dv[spliceout]
    retval = splice!(dv, spliceout, splicein...)
    isequal(retval, expval) ||
        error("splice!$(tuple(dv, spliceout, splicein...)) gave incorrect return value $retval (expected $expval)")
    splice!(v, spliceout, splicein...)
    for i = 1:length(v)
        isequal(dv[i], v[i]) ||
            error("splice!$(tuple(dv, spliceout, splicein...)) gave incorrect output")
    end
end

function test_deleteat(dv, spliceout)
    dv = copy(dv)
    v = Any[x for x in dv]
    @test deleteat!(dv, spliceout) === dv
    deleteat!(v, spliceout)
    for i = 1:length(v)
        isequal(dv[i], v[i]) ||
            error("deleteat!($(tuple(dv, spliceout)) gave incorrect output")
    end
end

dv1 = @data [1.0, 2.0, NA, 2.0, NA, 3.0]
for dv in (dv1, convert(DataVector{Number}, dv1), convert(PooledDataArray, dv1))
    for spliceout in (2, 3, 2:3, 5:6)
        test_splice(dv, spliceout)
        test_deleteat(dv, spliceout)
        for splicein in ([], [3], @data([3]), @pdata([3]),
                         [3, 4, 5], [3., 4., 5.], @data([3, NA, 4]),
                         @pdata([3, NA, 4]), @data([NA, 3.0, 4.0]),
                         @pdata([NA, 3.0, 4.0]))
            test_splice(dv, spliceout, splicein)
        end
    end
end

## sizehint
sizehint!(@data([1.0, 2.0]), 5)
sizehint!(@pdata([1.0, 2.0]), 5)

end


text:
using LCMCore
using CaesarLCMTypes
using Rotations, CoordinateTransformations
using MeshCat
using TransformUtils
using AprilTags

vis = Visualizer()
open(vis)

setobject!(vis["home"], Triad(1.0))

detector = AprilTagDetector()
# setobject!(vis["tag01"], Triad(0.2))
# settransform!(vis, Translation(-0.5, -0.5, 0))

function tagPicture(vis::Visualizer,msg::apriltag_t, drawdict)
    # m=homography_to_pose(Array{Float64,2}(msg.homography),300,300,320,240)
    m = msg.pose
    tkey = "tag$(msg.id)"
    if !haskey(drawdict, tkey)
        setobject!(vis[tkey], Triad(0.2))
        drawdict[tkey] = 1
    end
    wTc=LinearMap(RotX(pi/2))
    cTt=Translation(m[1:3,4]...) ∘ LinearMap(Quat(m[1:3,1:3]))
    settransform!(vis[tkey], wTc ∘ cTt)
    #settransform!(vis[tkey], cTt)
    nothing
end

function typed_callback(channel::String, msg::apriltag_t, drawdict)
    tagPicture(vis,msg, drawdict)
end

drawdict = Dict{String, Int}()

lcm=LCM()
subscribe(lcm, "TEST_CHANNEL", (m,c)->typed_callback(m,c,drawdict), apriltag_t)

while true
    handle(lcm)

end


text:
"""
    RationalQuadraticKernel(; α::Real=2.0)

Rational-quadratic kernel with shape parameter `α`.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the rational-quadratic kernel with shape parameter
``\\alpha > 0`` is defined as
```math
k(x, x'; \\alpha) = \\bigg(1 + \\frac{\\|x - x'\\|_2^2}{2\\alpha}\\bigg)^{-\\alpha}.
```

The [`SqExponentialKernel`](@ref) is recovered in the limit as ``\\alpha \\to \\infty``.

See also: [`GammaRationalQuadraticKernel`](@ref)
"""
struct RationalQuadraticKernel{Tα<:Real} <: SimpleKernel
    α::Vector{Tα}
    function RationalQuadraticKernel(; alpha::T=2.0, α::T=alpha) where {T}
        @check_args(RationalQuadraticKernel, α, α > zero(T), "α > 0")
        return new{T}([α])
    end
end

@functor RationalQuadraticKernel

function kappa(κ::RationalQuadraticKernel, d²::T) where {T<:Real}
    return (one(T) + d² / (2 * first(κ.α)))^(-first(κ.α))
end

metric(::RationalQuadraticKernel) = SqEuclidean()

function Base.show(io::IO, κ::RationalQuadraticKernel)
    return print(io, "Rational Quadratic Kernel (α = $(first(κ.α)))")
end

"""
    GammaRationalQuadraticKernel(; α::Real=2.0, γ::Real=2.0)

γ-rational-quadratic kernel with shape parameters `α` and `γ`.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the γ-rational-quadratic kernel with shape
parameters ``\\alpha > 0`` and ``\\gamma \\in (0, 2]`` is defined as
```math
k(x, x'; \\alpha, \\gamma) = \\bigg(1 + \\frac{\\|x - x'\\|_2^{\\gamma}}{\\alpha}\\bigg)^{-\\alpha}.
```

The [`GammaExponentialKernel`](@ref) is recovered in the limit as ``\\alpha \\to \\infty``.

See also: [`RationalQuadraticKernel`](@ref)
"""
struct GammaRationalQuadraticKernel{Tα<:Real,Tγ<:Real} <: SimpleKernel
    α::Vector{Tα}
    γ::Vector{Tγ}
    function GammaRationalQuadraticKernel(;
        alpha::Real=2.0, gamma::Real=2.0, α::Real=alpha, γ::Real=gamma
    )
        @check_args(GammaRationalQuadraticKernel, α, α > zero(α), "α > 0")
        @check_args(GammaRationalQuadraticKernel, γ, zero(γ) < γ ≤ 2, "γ ∈ (0, 2]")
        return new{typeof(α),typeof(γ)}([α], [γ])
    end
end

@functor GammaRationalQuadraticKernel

function kappa(κ::GammaRationalQuadraticKernel, d::Real)
    return (one(d) + d^first(κ.γ) / first(κ.α))^(-first(κ.α))
end

metric(::GammaRationalQuadraticKernel) = Euclidean()

function Base.show(io::IO, κ::GammaRationalQuadraticKernel)
    return print(
        io, "Gamma Rational Quadratic Kernel (α = $(first(κ.α)), γ = $(first(κ.γ)))"
    )
end


text:
function check_distance_matrix(D::Matrix{<:Real})
    size(D, 1) > 1 || throw(DomainError(D, "Too few nodes"))
    if any(D .< 0)
        throw(DomainError(D, "Negative entries in distance matrix"))
    end
    if any(diag(D) .!= 0.0)
        throw(DomainError(D, "Diagonal entry not 0"))
    end
end

"""
    nj!(D::Matrix{Float64}, names::AbstractVector{String}=String[];
        force_nonnegative_edges::Bool=false)

Construct a phylogenetic tree from the input distance matrix and
vector of names (as strings), using the Neighbour-Joinging algorithm
([Satou & Nei 1987](https://doi.org/10.1093/oxfordjournals.molbev.a040454)).
The order of the `names` argument should match that of the row (and
column) of the matrix `D`.
With `force_nonnegative_edges` being `true`, any negative edge length
is changed to 0.0 (with a message).

Warning: `D` is modified.
"""
function nj!(D::Matrix{Float64}, names::AbstractVector{String}=String[];
             force_nonnegative_edges::Bool=false)

    check_distance_matrix(D)
    n = size(D, 1)              # number of species

    # when no names arg is supplied
    if isempty(names)
        names = string.(1:n)
    end

    # create empty network with n unconnected leaf nodes
    nodes = map(function(i)
                node = Node(i, true)
                node.name = names[i]
                return node
                end,
                1:n)

    net = HybridNetwork(nodes, Edge[])
    # an array of Node s.t. active_nodes[i] would correspond to the
    # ith entry in distance matrix D at each iteration
    active_nodes = copy(nodes)

    neglenp = 0  # number of negative edge lengths

    while n > 2
        # compute Q matrix and find min
        # warning: here i and j are indeces for D, not actual id's
        sums = sum(D, dims = 1)
        cur = Inf
        min_index = (0, 0)
        for i = 1:n
            for j = i:n
                if j != i
                    qij = (n-2) * D[i,j] - sums[i] - sums[j]
                    if qij < cur
                        cur = qij
                        min_index = (i, j)
                    end
                end
            end
        end

        # connect the nodes, compute the length of the edge
        (i, j) = min_index
        dik = D[i,j] / 2 + (sums[i] - sums[j]) / (2 * (n - 2))
        djk = D[i,j] - dik
        # force negative lengths to 0, if any
        if dik < 0.0
            neglenp += 1
            if force_nonnegative_edges dik = 0.0; end
        end
        if djk < 0.0
            neglenp += 1
            if force_nonnegative_edges djk = 0.0; end
        end

        # create new edges and node, update tree
        edgenum = net.numEdges
        eik = Edge(edgenum + 1, dik) # edge length must be Float64 for Edge()
        ejk = Edge(edgenum + 2, djk)
        node_k = Node(net.numNodes+1, false, false, [eik, ejk]) # new node
        node_i = active_nodes[i]
        node_j = active_nodes[j]
        setNode!(eik, Node[node_i, node_k])
        setNode!(ejk, Node[node_j, node_k])
        setEdge!(node_i, eik)
        setEdge!(node_j, ejk)
        pushEdge!(net, eik)
        pushEdge!(net, ejk)
        pushNode!(net, node_k)

        # update map and D
        # replace D[l, i] with D[l, k], delete D[ , j]
        for l in 1:n
            if !(l in [i j])
                D[l, i] = (D[l,i] + D[j,l] - D[i,j]) / 2
                D[i, l] = D[l, i]
            end
        end
        newidx = filter(u->u!=j, 1:n) # index 1:n\{j}
        D = view(D, newidx, newidx)

        # update active_nodes
        active_nodes[i] = node_k
        active_nodes = view(active_nodes, newidx)

        n = n - 1
    end

    # base case
    if D[1,2] < 0.0
        neglenp += 1
        if force_nonnegative_edges D[1,2] = 0.0; end
    end
    node1 = active_nodes[1]
    node2 = active_nodes[2]
    newedge = Edge(net.numEdges+1, D[1,2])
    setNode!(newedge, [node1, node2]) # isChild1 = true by default: nodes are [child, parent]
    if node1.number == net.numNodes
        newedge.isChild1 = false # to direct the edge from the last created node to the other node
    end
    setEdge!(node1, newedge)
    setEdge!(node2, newedge)
    pushEdge!(net, newedge)

    net.root = net.numNodes # to place the root at the last created node, which is internal

    # report on number of negative branches
    if neglenp > 0
        infostr = (force_nonnegative_edges ?
                   "$neglenp branches had negative lengths, reset to 0" :
                   "$neglenp branches have negative lengths" )
        @info infostr
    end

    return net
end

"""
    nj(D::DataFrame; force_nonnegative_edges::Bool=false)

Construct a tree from a distance matrix by neighbor joining, where
`D` is a `DataFrame` of the distance matrix, with taxon names taken
from the header of the data frame.
The rows are assumed to correspond to tips in the tree in the same order
as they do in columns.
With `force_nonnegative_edges` being `true`, any negative edge length
is changed to 0.0 (with a message).

For the algorithm, see
[Satou & Nei 1987](https://doi.org/10.1093/oxfordjournals.molbev.a040454).

See [`nj!`](@ref) for using a matrix as input.
"""
function nj(D::DataFrame; force_nonnegative_edges::Bool=false)
    nj!(Matrix{Float64}(D), names(D); force_nonnegative_edges=force_nonnegative_edges)
end


text:
# diffs.jl

import MIRT
using LinearMapsAA: LinearMapAA
using Random: seed!
using BenchmarkTools: @btime

# old 2D versions for comparison with the new N-D versions

"""
    d = diff2d_forw_old(X)

2D finite differences along both dimensions, for anisotropic TV regularization.
Performs the same operations as
d = ``[(I_N \\otimes D_M); (D_N \\otimes I_M)] X[:]``
where ``D_N`` denotes the ``N-1 \\times N`` 1D finite difference matrix
and ``\\otimes`` denotes the Kronecker product,
but does it efficiently
without using `spdiagm` (or any `SparseArrays` function).

in
- `X` `M × N` array (typically a 2D image).
It cannot be a Vector! (But it can be a `M×1` or `1×N` 2D array.)

out
- `d` vector of length `N*(M-1) + (N-1)*M`
"""
function diff2d_forw_old(x::AbstractMatrix{<:Number})
    return [diff(x,dims=1)[:]; diff(x,dims=2)[:]]
end

"""
    z = diff2d_adj_old(d, M, N; out2d=false)

Adjoint of 2D finite differences along both dimensions.
Performs the same operations as
``z = [(I_N \\otimes D_M); (D_N \\otimes I_M)]' * d``
where `D_N` denotes the `N-1 × N` 1D finite difference matrix
and `\\otimes` denotes the Kronecker product, but does it efficiently
without using `spdiagm` (or any `SparseArrays` function).

in
- `d` vector of length `N*(M-1) + (N-1)*M`
- `M,N` desired output size

option
- `out2d` if true then return `M × N` array, else `M*N` vector

out
- `z` `M*N` vector or `M × N` array (typically a 2D image)

"""
function diff2d_adj_old(d::AbstractVector{<:Number}, M::Int, N::Int ; out2d=false)

    length(d) != N*(M-1) + (N-1)*M && throw("length(d)")

    dx = d[1:(N*(M-1))]
    dx = reshape(dx, M-1, N)
    dy = d[1+(N*(M-1)):end]
    dy = reshape(dy, M, N-1)
    zx = [-transpose(dx[1,:]);
        (@views dx[1:end-1,:] - dx[2:end,:]);
        transpose(dx[end,:])]
    zy = [-dy[:,1] (@views dy[:,1:end-1] - dy[:,2:end]) dy[:,end]]
    z = zx + zy # M by N

    return out2d ? z : z[:] # array or vector?
end

"""
    T = diff2d_map_old(M::Int, N::Int)
"""
function diff2d_map_old(M::Int, N::Int)
    return LinearMapAA(
        x -> diff2d_forw_old(reshape(x,M,N)),
        d -> diff2d_adj_old(d, M, N),
        (N*(M-1)+M*(N-1), N*M), (name="diff2_map",),
    )
end

# timing runs
function diff_map_time( ; M::Int=2^7, N::Int=2^7+1)
    seed!(0)
    x = randn(M * N)
    T2d = diff2d_map_old(M, N)
    Tnd = MIRT.diff_map((M, N))

    # check consistency
    d2d = T2d * x
    dnd = Tnd * x
    @assert d2d == dnd

    d = d2d
    y2d = T2d' * d
    ynd = Tnd' * d
    @assert y2d == ynd[:]

    # time diff_forw
    @btime $T2d * $x
    @btime $Tnd * $x

    # time diff_adj
    @btime $T2d' * $d
    @btime $Tnd' * $d
    nothing
end

#= ir71 with 1.4.2
  98.381 μs (16 allocations: 762.73 KiB)
  58.344 μs (18 allocations: 508.77 KiB)
  159.237 μs (39 allocations: 1.00 MiB)
  102.423 μs (115 allocations: 384.05 KiB)
=#

diff_map_time()


text:
# =========================================================================
#    ENVIRONMENT SETTINGS
# =========================================================================

#OPTIMIZER SETTING
const default_optimizer = Ref{Any}(nothing)
set_optimizer!(O) = default_optimizer[] = O
get_optimizer() = default_optimizer[]
reset_optimizer() = default_optimizer[] = nothing

#SLACK SURP COEFF SETTING
const slack_coeff = Ref{Any}(nothing)
const surp_coeff = Ref{Any}(nothing)
set_slack_coeff!(O) = slack_coeff[] = O
set_surp_coeff!(O) = surp_coeff[] = O
su_C() = surp_coeff[]
sl_C() = slack_coeff[]

#SILENCE SETTING
const verbosity = Ref{Any}(nothing)
silence!(O) = verbosity[] = O
silent() = verbosity[]


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie
CairoMakie.activate!() #HIDE
let
    fig = Figure(resolution = (700,450))
    ax1 = Axis(fig[1, 1], yticklabelcolor = :black, rightspinevisible = false)
    ax2 = Axis(fig[1, 1], yticklabelcolor = :orangered, yaxisposition = :right,
        rightspinecolor = :orangered, ytickcolor = :orangered)
    lines!(ax1, 0..10, x -> x, color = :black)
    lines!(ax2, 0..10, x -> exp(-x), color = :orangered)
    hidespines!(ax2, :l,:b, :t)
    hidexdecorations!(ax2)
    fig
    save(joinpath(@__DIR__, "output", "twinAxis.png"), current_figure(), px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using ProgressMeter
using MATLAB
using Printf
include("PSS_function.jl")

H=range(.11,stop=0.15,length=3);
max_hit=100
@showprogress 1 "Computing..." for Energy in H

mat"""
figure; hold on;
set(gcf, 'Position',  [0, 0, 2000, 1500])
"""
# set(gca,'XTick',[])
# set(gca,'YTick',[])
# set(gca,'xticklabel',[])
# set(gca,'yticklabel',[])
# set(gca,'XColor','none')
# set(gca,'YColor','none')
# set(gca,'XColor','none')
# set(gca,'YColor','none')
# location="~/Desktop/PSS_MOVIE/"
location="~/MATLAB-Drive/PSS_LEAP_FROG_CENTER/"
h=replace(@sprintf("%.7f",Energy),"."=>"_")
file_name=location*h*".fig"
println(file_name)

# n_iter_Q=100;
# Q_start=-1
# Q_end=1
# n_iter_P=401;
# P_start=-2
# P_end=2
# t_end=1000

n_iter_Q=10;#50
Q_start=-.25
Q_end=.25
n_iter_P=11;#51
P_start=-.015
P_end=.015
t_end=1e3
#create an empty list to store dH
H_differences=zeros(0)
#julia's version of linspace
ArrP=range(P_start,stop=P_end,length=n_iter_P)
ArrQ=range(Q_start,stop=Q_end,length=n_iter_Q)

#defining colors for PSS
COLOR=["b", "g", "c", "m", "y", "k", "r"]
# COLOR=['b','g','c','m','y','k','r']
# axis([ -.003,.003,-.08, .08])
mat"axis([ -.03,.03,-.4,.4 ])"
#
# mat"axis([ -3,3,-1.5, 1.5])"
# @showprogress 1 "Computing..."
for j=1:n_iter_P
    for k=1:n_iter_Q
        Q,P,dH=PSS_function(ArrQ[k], ArrP[j], Energy, t_end, max_hit)
        if dH<5e-5 #make sure dH is ok

            #push new dH to list
            # push!(H_differences,dH)
            #choose new color
            current_color=COLOR[mod(k,length(COLOR))+1]

            #take all symmetric combinations
            Q_PSS=vcat(Q)
            P_PSS=vcat(P)
            # Qy=A[2,:];
            # Px=A[4,:];
            mat"plot($P_PSS,$Q_PSS,'.','MarkerSize',.1,'color',$current_color); hold on;"
            # mat"plot($P_PSS,$Q_PSS,'k.','MarkerSize',.1); hold on;"
        end
    end
end
mat"savefig($file_name)"
mat"close"
end


text:
using StrainRecon
using Base.Test

@testset "entropyM" begin

M  = [1 1 1; 1 1 0; 1 0 1; 1 0 0; 0 1 1; 0 1 0; 0 0 1;]
(m,n) = size(M)
R     = StrainRecon.sdiag(ones(m)*100);
atrue = [.6; 0.3; 0.1]
Data  = M*atrue

param = getStrainReconParam(Data,n,R=R)

entropy = entropyM(param,atrue)

@test entropy >= 0 && entropy <= m*n
end


text:
init_cache(model, method::ExplicitDiffusion) =
    ExplicitDiffusionCache(copy(model.u0), copy(model.u0), copy(model.parameters.tspan[1]))
function init_cache(model, method::CrankNickolsonDiffusion)
    u     = copy(model.u0)
    uprev = copy(model.u0)
    rhs   = copy(model.u0)

    a1, a2 = init_CN(model)

    t = copy(model.parameters.tspan[1])

    return CNDiffusionCache(u, uprev, rhs, a1, a2, t);
end
function init_cache(model, method::ADIDiffusion)
    u     = copy(model.u0)
    uhalf = copy(model.u0)
    uprev = copy(model.u0)
    rhs   = copy(model.u0)

    a1x, a1y, a2x, a2y = init_ADI(model)

    t = copy(model.parameters.tspan[1])

    return ADIDiffusionCache(u, uhalf, uprev, rhs, a1x, a1y, a2x, a2y, t);
end


text:
using Test
using SafeTestsets

include("diff/diff_test.jl")

@safetestset "Factorization Test" begin
    include("factorization_test.jl")
end

include("initialization/initialization_test.jl")

@safetestset "Shape Tests" begin
    include("shape_test.jl")
end

@safetestset "Dynamics Tests" begin
    include("dynamics/dynamics_test.jl")
end

@safetestset "Example Tests" begin
    include("example_test.jl")
end


text:
__precompile__(false)

using Revise

import Documenter

using BiTWeather

Documenter.makedocs(
    sitename="BiTWeather.jl",
    modules=[BiTWeather],
    clean=true,
    format=Documenter.HTML(
        prettyurls=false
    ),
    pages = [
        "Introduction" => "index.md",
        "Chill Models" => [
            "chill_CumulativeChillHour.md",
            "chill_CumulativeChillUnit.md",
            "chill_CumulativeChillPortion.md",
            "chill_MeanTemperature.md"
        ],
        "API" => [
            "api_types.md",
            "api_functions.md",
            "api_constants.md"
        ],
        "Background" => "background.md"
    ]
)

max_stars_repo_name: pebender/BiTWeather

max_issues_repo_name: pebender/BiTWeather

max_forks_repo_name: pebender/BiTWeather


text:
"""
$(SIGNATURES)

Sample particle velocity from Maxwellian. Unconstraint methods use Box-Muller sampling,
while constraint one depends on Distributions.jl.
"""
function sample_maxwell(λ::Real, v = 0.0::Real)
    rd1 = rand()
    rd2 = rand()

    return sqrt(-log(rd1) / λ) * sin(2.0 * π * rd2) + v
end

"""
$(SIGNATURES)
"""
function sample_maxwell(prim::T) where {T<:AA{<:Real,1}}

    u = prim[2]
    if length(prim) == 3
        v = 0.0
        w = 0.0
    elseif length(prim) == 4
        v = prim[3]
        w = 0.0
    elseif length(prim) == 5
        v = prim[3]
        w = prim[4]
    end

    rd1 = rand(3)
    rd2 = rand(3)

    return @. sqrt(-log(rd1) / prim[end]) * sin(2.0 * π * rd2) + [u, v, w]

end

"""
$(SIGNATURES)

Truncated sampling with Distributions.jl
"""
function sample_maxwell(λ::Real, l, u)
    d = truncated(Normal(0.0, sqrt(1.0 / λ)), l, u)
    return rand(d)
end

"""
$(SIGNATURES)
"""
function sample_maxwell(λ::Real, v, l, u)
    d = truncated(Normal(0.0, sqrt(1.0 / λ)), l, u)
    return rand(d) + v
end

"""
$(SIGNATURES)
"""
function sample_maxwell(prim::T, l, u) where {T<:AA{<:Real,1}}
    u = prim[2]
    if length(prim) == 3
        v = 0.0
        w = 0.0
    elseif length(prim) == 4
        v = prim[3]
        w = 0.0
    elseif length(prim) == 5
        v = prim[3]
        w = prim[4]
    end

    d = truncated(Normal(0.0, sqrt(1.0 / prim[end])), l, u)
    return rand(d, 3) + [u, v, w]
end

"""
$(SIGNATURES)

Calculate the time interval until next collision
"""
next_collision_time(τ) = -τ * log(rand())

"""
$(SIGNATURES)

Calculate traveling time to edges
"""
function boundary_time(x, v::Real, xL, xR)
    if v < 0.0
        tb = (xL - x) / v
    elseif v > 0.0
        tb = (xR - x) / v
    else
        tb = 1e8
    end

    return tb
end

"""
$(SIGNATURES)
"""
function boundary_time(x, v::T, xL, xR) where {T<:AA{<:Real,1}}
    if v[1] < 0.0
        tb = (xL - x) / v[1]
    elseif v[1] > 0.0
        tb = (xR - x) / v[1]
    else
        tb = 1e8
    end

    return tb
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL; GL = ViewerGL

#//////////////////////////////////////////////////////////////////////
# ----- solid ridge construction --------------------------------------
#//////////////////////////////////////////////////////////////////////

# parametric solid generation -----------------------------------------
# ridge = extrusion solid ---------------------------------------------
# whose section is a circle approximation -- with #sides ≥ 3 ----------

#-- topological algorithm (distribution of product over skeletons)-----
n = 10
circle1 = Lar.circle()([n])
interval1 = Lar.qn(1)([3])

circle0 = (circle1[1], [[k] for k=1:size(circle1[1],2)])
interval0 = (interval1[1], [[k] for k=1:size(interval1[1],2)])

cyl1a = Lar.larModelProduct(circle0,interval1)
cyl1b = Lar.larModelProduct(circle1,interval0)
cyl1 = ( cyl1a[1],append!(cyl1a[2], cyl1b[2]) )

GL.VIEW([ GL.GLGrid(cyl1[1], cyl1[2],GL.COLORS[1],1) ]);

cyl2a = Lar.larModelProduct(circle1,interval1)
circle2 = (circle1[1],[collect(1:n)])
cyl2b = Lar.larModelProduct(circle2,interval0)
cyl2 = ( cyl2a[1],append!(cyl2a[2], cyl2b[2]) )

#-- LAR model building -----------------------------------------------
V,EV = cyl1
V,FV = cyl2
rod = V,FV,EV

#///////// single rod instance /////////////////////////////

cop_EV = convert(Lar.ChainOp, Lar.coboundary_0( EV::Lar.Cells ));
cop_FE = Lar.coboundary_1( V, FV::Lar.Cells, EV::Lar.Cells );
W = convert(Lar.Points, V');
#------ E^3 arrangement ------------------------------------
V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement( W, cop_EV, cop_FE)

#------ visualization --------------------------------------
W = convert(Lar.Points, V')
V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF)
GL.VIEW(GL.GLExplode(V,FVs,1.,1.,1.,99,1));

#/////// triple instance ///////////////////////////////////

#------ model generation --------------------------------------
tris = Lar.Struct([
	rod,
	Lar.Struct([Lar.r(pi/2,0,0), rod ]) ,
	Lar.Struct([Lar.r(0,pi/2,0), rod ])
])
V,FV,EV = Lar.struct2lar( tris )
GL.VIEW([ GL.GLGrid( V, EV, GL.COLORS[1],1) ]);

#------ E^3 arrangement ------------------------------------
cop_EV = convert(Lar.ChainOp, Lar.coboundary_0(EV::Lar.Cells));
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');
V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement( W, cop_EV, cop_FE)

#------ visualization --------------------------------------
W = convert(Lar.Points, V')
V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF)
GL.VIEW(GL.GLExplode(V,FVs,1.,1.,1.,99,1));

#///////// Exploded views of 2-sells and boundaries of 2-cells ///

GL.VIEW(GL.GLExplode(V,FVs,1.5,1.5,1.5,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));

#////////// Algebra atoms (generated 3-cells) //////////////////

meshes = GL.GLExplode(V,CVs[2:end],4,4,4,99,0.5);
GL.VIEW( push!( meshes, GL.GLFrame) );


text:
const N_CYCLES = 6

function count_neighbors(grid::Array{Int,3}, row::Int, col::Int, layer::Int)::Int
    count::Int = 0
    for x in row - 1:row + 1
        for y in col - 1:col + 1
            for z in layer - 1:layer + 1
                if x == row && y == col && z == layer
                    continue
                else
                    count += grid[x, y, z]
                end
            end
        end
    end
    return count
end

function run(s)
    offset = N_CYCLES + 1
    grid = zeros(Int, (8 + 2 * offset, 8 + 2 * offset, 1 + 2 * offset))
    for (row, line) in enumerate(eachline(IOBuffer(s)))
        for (col, c) in enumerate(line)
            if c == '#'
                grid[row + offset, col + offset, 1 + offset] = 1
            end
        end
    end
    old_grid = copy(grid)
    n_neighbors::Int = 0
    for cycle in 1:N_CYCLES
        old_grid .= grid
        for row in offset + 1 - cycle:offset + 8 + cycle
            for col in offset + 1 - cycle:offset + 8 + cycle
                for layer in offset + 1 - cycle:offset + 1 + cycle
                    n_neighbors = count_neighbors(old_grid, row, col, layer)
                    if (old_grid[row, col, layer] != 0) && (n_neighbors != 2) && (n_neighbors != 3)
                        grid[row, col, layer] = 0
                    elseif (old_grid[row, col, layer] == 0) && (n_neighbors == 3)
                        grid[row, col, layer] = 1
                    end
                end
            end
        end
    end
    return sum(grid)
end

#########################################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
module helper
#Solves for steady state distribution of workforce across age stages
function solveStages(var, fixedParam)
    NC_f    = var[1];
    YC_f    = var[2];
    PL_f    = var[3];
    D_f     = var[4];

    NC_m    = var[5];
    YC_m    = var[6];
    PL_m    = var[7];
    D_m     = var[8];

    # 4 equations for women
    eq1f = NC_f + YC_f + PL_f + D_f - 0.5;
    eq2f = fixedParam.chi*(NC_f+YC_f) - (fixedParam.etaf+fixedParam.gamma)*PL_f; # flows into and out of PL
    eq3f = fixedParam.etaf*PL_f - (fixedParam.chi+fixedParam.gamma)*YC_f; # flows into and out of YC
    eq4f = fixedParam.gamma*(NC_f+YC_f+PL_f) - fixedParam.phi*D_f;

    # 4 equations for men
    eq1m = NC_m + YC_m + PL_m + D_m - 0.5;
    eq2m = fixedParam.chi*(NC_m+YC_m) - (fixedParam.etam+fixedParam.gamma)*PL_m; # flows into and out of PL
    eq3m = fixedParam.etam*PL_m - (fixedParam.chi+fixedParam.gamma)*YC_m; # flows into and out of YC
    eq4m = fixedParam.gamma*(NC_m+YC_m+PL_m) - fixedParam.phi*D_m;

    return [eq1f, eq2f, eq3f, eq4f, eq1m, eq2m, eq3m, eq4m];
end
export solveStages
end


text:
function solve()
    n, s, d = [parse(Int, x) for x in split(readline())]
    for _=1:n
        x, y = [parse(Int, x) for x in split(readline())]
        if x < s && y > d
            return "Yes"
        end
    end
    "No"
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Documenter, TextStylometry

makedocs(
    format = :html,
    sitename = "Text Stylometry",
    modules = [TextStylometry],
    pages = [
        "index.md",
        "Corpus" => "corpus.md",
        "Features" => "features.md",
        "Measures" => "measures.md",
        "Bootstrap measures" => "bootstrap.md",
    ]
)

deploydocs(
    repo   = "github.com/borh/TextStylometry.jl.git",
    target = "build",
    julia  = "0.6",
    deps   = nothing, # Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    make   = nothing
)


text:
using Documenter
using RefractiveIndexDatabase

makedocs(
    sitename = "RefractiveIndexDatabase.jl",
    authors="Adam Fekete <adam.fekete@unamur.be> and contributors",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        canonical="https://fekad.github.io/RefractiveIndexDatabase.jl"
    ),
    modules = [RefractiveIndexDatabase],
    pages=[
        "Home" => "index.md",
        "Usage" => "usage.md",
        "Formualas" => "formulas.md",
    ]
)

deploydocs(
    repo = "github.com/fekad/RefractiveIndexDatabase.jl",
    devbranch = "main"
)


text:
using MIDI
using Base.Test

include("variablelength.jl")
include("sysexevent.jl")
include("note.jl")
include("midievent.jl")
include("metaevent.jl")
include("miditrack.jl")


text:
function pisum()
sum = 0.0
for j = 1:500
sum = 0.0
for k = 1:10000
sum += 1.0/(k*k)
end
end
sum
end

max_stars_repo_path: Module 2/Chapter01/p1.jl

max_issues_repo_path: Module 2/Chapter01/p1.jl

max_forks_repo_path: Module 2/Chapter01/p1.jl


text:
#####
##### Benchmark.jl
##### Script to run a new benchmark on a previous training session
##### (requires training with option `save_intermediate=true`)
#####

using AlphaZero

include("games.jl")
GAME = get(ENV, "GAME", "double-dummy-full-state")
SelectedGame = GAME_MODULE[GAME]
using .SelectedGame: Game, Training

SESSION_DIR = "sessions/$GAME"

alphazero = Benchmark.Full(Training.arena.mcts)
baselines = [Benchmark.Solver(ϵ=0), Training.minmax_baseline]

make_duel(player, baseline) =
  Benchmark.Duel(
    player,
    baseline,
    num_games=1,
    reset_every=1,
    flip_probability=0.0,
    color_policy=BASELINE_WHITE)

benchmark = [make_duel(alphazero, Benchmark.Solver(ϵ=0))]
label = "full"

AlphaZero.UserInterface.run_new_benchmark(
  Game, Training.Network{Game}, SESSION_DIR,
  label, benchmark,
  params=Training.params, itcmax=nothing)

max_stars_repo_name: btheodorou/Bridge-NN-Research

max_issues_repo_name: btheodorou/Bridge-NN-Research

max_forks_repo_name: btheodorou/Bridge-NN-Research


text:
__precompile__()

module FLAC

using FileIO
using FLAC_jll

export StreamMetaData,
       InfoMetaData,
       VorbisCommentMetaData,
       PaddingMetaData,
       ApplicationMetaData,
       SeekTableMetaData,
       CueSheetMetaData,
       StreamDecoderPtr,
       StreamEncoderPtr,
       initfile!,
       FLACDecoder,
       seek,
       read,
       size,
       length

import Base: read, seek, size, length

include("metadata.jl")
include("format.jl")
include("decoder.jl")
include("encoder.jl")

function __init__()
    init_decoder_cfunctions()
    init_encoder_cfunctions()
end

end # module


text:
"""
    IsotropicDiffusivity{N, K}

Holds viscosity and diffusivities for models with prescribed isotropic diffusivities.
"""
struct IsotropicDiffusivity{N, K} <: AbstractIsotropicDiffusivity
    ν :: N
    κ :: K
end

"""
    IsotropicDiffusivity(; ν=ν₀, κ=κ₀)

Returns parameters for an isotropic diffusivity model with viscosity `ν`
and thermal diffusivities `κ` for each tracer field in `tracers`
`ν` and the fields of `κ` may be constants or functions of `(x, y, z, t)`, and
may represent molecular diffusivities in cases that all flow
features are explicitly resovled, or turbulent eddy diffusivities that model the effect of
unresolved, subgrid-scale turbulence.
`κ` may be a `NamedTuple` with fields corresponding
to each tracer, or a single number to be a applied to all tracers.

By default, a molecular viscosity of `ν₀ = 1.05×10⁻⁶` m² s⁻¹ and a molecular thermal
diffusivity of `κ₀ = 1.46×10⁻⁷` m² s⁻¹ is used for each tracer. These molecular values are
the approximate viscosity and thermal diffusivity for seawater at 20°C and 35 psu,
according to Sharqawy et al., "Thermophysical properties of seawater: A review of existing
correlations and data" (2010).
"""
function IsotropicDiffusivity(FT=Float64; ν=ν₀, κ=κ₀)
    if ν isa Number && κ isa Number
        κ = convert_diffusivity(FT, κ)
        return IsotropicDiffusivity(FT(ν), κ)
    else
        return IsotropicDiffusivity(ν, κ)
    end
end

function with_tracers(tracers, closure::IsotropicDiffusivity)
    κ = tracer_diffusivities(tracers, closure.κ)
    return IsotropicDiffusivity(closure.ν, κ)
end

calculate_diffusivities!(K, arch, grid, closure::IsotropicDiffusivity, args...) = nothing

@inline function ∇_κ_∇c(i, j, k, grid, clock, closure::IsotropicDiffusivity,
                        c, ::Val{tracer_index}, args...) where tracer_index

    @inbounds κ = closure.κ[tracer_index]

    return ∂ⱼκᵢⱼ∂ᵢc(i, j, k, grid, clock, κ, c)
end

@inline ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, clock, closure::IsotropicDiffusivity, U, args...) =
    ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, clock, closure.ν, U.u)

@inline ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, clock, closure::IsotropicDiffusivity, U, args...) =
    ∂ⱼνᵢⱼ∂ᵢv(i, j, k, grid, clock, closure.ν, U.v)

@inline ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, clock, closure::IsotropicDiffusivity, U, args...) =
    ∂ⱼνᵢⱼ∂ᵢw(i, j, k, grid, clock, closure.ν, U.w)

Base.show(io::IO, closure::IsotropicDiffusivity) =
    print(io, "IsotropicDiffusivity: ν=$(closure.ν), κ=$(closure.κ)")


text:
wiener_randn!() = nothing

function __init__()
    @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
        wiener_randn!(rng::AbstractRNG,rand_vec::CuArrays.CuArray) = randn!(rand_vec)
    end

    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
        wiener_randn!(rng::AbstractRNG,rand_vec::CUDA.CuArray) = randn!(rand_vec)
    end

    @require ReverseDiff="37e2e3b7-166d-5795-8a7a-e32c996b4267" begin
        @inline function DiffEqNoiseProcess.wiener_randn(rng::Random.AbstractRNG,proto::ReverseDiff.TrackedArray)
          ReverseDiff.track(convert.(eltype(proto.value),randn(rng,size(proto))))
        end
        @inline function DiffEqNoiseProcess.wiener_randn!(rng::AbstractRNG,rand_vec::Array{<:ReverseDiff.TrackedReal})
            rand_vec .= ReverseDiff.track.(randn.((rng,),typeof.(DiffEqBase.value.(rand_vec))))
        end
        @inline function DiffEqNoiseProcess.wiener_randn!(rng::AbstractRNG,rand_vec::AbstractArray{<:ReverseDiff.TrackedReal})
            rand_vec .= ReverseDiff.track.(randn.((rng,),typeof.(DiffEqBase.value.(rand_vec))))
        end
    end
end


text:
#=
diffs.jl
finite differences
2019-03-06 Jeff Fessler, University of Michigan
2020-06 N-D version by Steven Whitaker
=#

export diff_map, diff_forw, diff_adj

using LinearMapsAA: LinearMapAA

# size after apply `diff` along `dim`
diff_size(N::Dims, dim::Int) = (N[1:dim-1]..., N[dim] - 1, N[dim+1:end]...)

# corresponding length
diff_length(N::Dims, dim::Int) = prod(diff_size(N, dim))

# check argument validity
# dims should be Int or AbstractVector{Int} or Dims
# but using duck typing for simplicity
function diff_check(N::Dims, dims)
    (length(unique(dims)) != length(dims)) &&
        throw(ArgumentError("non-unique dims $dims"))
    map(dim -> diff_check(N, dim), dims) # check each dim
end

function diff_check(N::Dims, dim::Int)
    !(1 <= dim <= length(N)) && throw(ArgumentError("dim range $dim"))
    (N[dim] == 1) && throw(ArgumentError("invalid size $N for dim $dim"))
end

"""
    d = diff_forw(X ; dims = 1:ndims(X))

Finite differences along one or more dimensions of an array,
e.g., for anisotropic TV regularization.

By default performs the same operations as
``d = [(I_{N_d} \\otimes \\cdots \\otimes D_{N_1}); \\dots; (D_{N_d} \\otimes \\cdots \\otimes I_{N_1})] vec(X)``
where ``D_N`` denotes the `N-1 × N` 1D finite difference matrix
and `⊗` denotes the Kronecker product, but does it efficiently
without using `spdiagm` (or any `SparseArrays` function).

Input dimension `N` must exceed `1` for each dimension specified by `dims`.

in
- `X` `N_1 × ... × N_d` array (typically an N-D image).

option
- `dims` dimension(s) for performing finite differences; default `1:ndims(X)`
must have unique elements and be a nonempty subset of `1:ndims(X)`

out
- `d` vector of default length `N_d*...*(N_1-1) + ... + (N_d-1)*...*N_1`
"""
function diff_forw(x::AbstractArray{<:Number,D} ; dims = 1:D) where {D}
    diff_check(size(x), dims)
    return reduce(vcat, vec(diff(x, dims = d)) for d in dims)
end

"""
    Z = diff_adj(dx, N::Dims{D} ; dims = 1:D)

Adjoint of finite differences of arrays along one or more dimensions.
By default performs the same operations as
``vec(Z) = [(I_{N_d} \\otimes \\cdots \\otimes D_{N_1}); \\dots; (D_{N_d} \\otimes \\cdots \\otimes I_{N_1})]' * d``
where `D_N` denotes the `N-1 × N` 1D finite difference matrix
and `⊗` denotes the Kronecker product,
but does it efficiently without using `spdiagm` (or any `SparseArrays` function).

in
- `dx` vector of typical length `N_d*...*(N_1-1) + ... + (N_d-1)*...*N_1`
- `N::Dims` desired output size

option
- `dims` dimension(s) for performing adjoint finite differences; default `1:ndims(X)`

out
- `Z` `N_1 × ... × N_d` array by default

"""
function diff_adj(d::AbstractVector{<:Number}, N::Dims{D} ; dims = 1:D) where {D}

    length(d) != sum(diff_length(N,dim) for dim in dims) && throw("length(d)")

    z = zeros(eltype(d), N)
    for (i, dim) in enumerate(dims)
        if i == 1
            di = @view(d[1:diff_length(N,dim)])
        else
            start = 1 + sum(diff_length(N,n) for n in dims[1:i-1])
            len = diff_length(N,dim)
            di = @view(d[start:start+len-1])
        end
        di = reshape(di, diff_size(N,dim))
        slice1 = selectdim(z, dim, 1)
        slice1 .-= selectdim(di, dim, 1)
        slicen = selectdim(z, dim, 2:N[dim]-1)
        slicen .+= selectdim(di, dim, 1:N[dim]-2) - selectdim(di, dim, 2:N[dim]-1)
        sliceN = selectdim(z, dim, N[dim])
        sliceN .+= selectdim(di, dim, N[dim] - 1)
    end

    return z
end

# backward compatibility, undocumented because deprecated
function diff2d_forw(x::AbstractMatrix{<:Number})
    isinteractive() && @warn("diff2d_forw is deprecated; use diff_forw")
    return diff_forw(x)
end

function diff2d_adj(d::AbstractVector{<:Number}, M::Int, N::Int ; out2d=false)
    isinteractive() && @warn("diff2d_adj is deprecated; use diff_adj")
    tmp = diff_adj(d, (M, N))
    return out2d ? tmp : vec(tmp)
end

function diff2d_map(M::Int, N::Int)
    isinteractive() && @warn("diff2d_map is deprecated; use diff_map")
    return diff_map((M, N))
end

"""
    T = diff_map(N::Dims{D} ; dims = 1:D)

in
- `N::Dims` image size

out
- `T` `LinearMapAA` object for computing finite differences via `T*x`
"""
function diff_map(N::Dims{D} ; dims = 1:D) where {D}
    diff_check(N, dims)
    return LinearMapAA(
        x -> diff_forw(reshape(x, N), dims=dims),
        d -> vec(diff_adj(d, N, dims=dims)),
        (sum(diff_length(N,dim) for dim in dims), prod(N)),
        (name="diff_map", dims=dims),
    )
end


text:
module MarketBasket

export  Transactions,
        Cluster,
        getindex,
        length,
        items,
        refine!

include("clustering.jl")

end # module


text:
#
# create Lienar FEM Hat functions with ReLU
#
# centred at x=0, width = 1.0 => spacing = 0.5
#
using Flux, NNlib, DiffEqFlux
using Plots
import SEM
#--------------------------------------#
x = SEM.linspace(-2,2,1000)

fm1 = @. 2.0relu(x+0.5)
fm0 = @. 4.0relu(x-0.0)
fp1 = @. 2.0relu(x-0.5)

hat = fm1 - fm0 + fp1

p = plot()
p = plot!(title="FEM Hat Function",xlims=[-1,1],ylims=[0,1.5])
p = plot!(x,hat,width=4,color=:black,style=:solid,legend=false)
savefig(p,"reluHat.png")

# ReLU
p=plot(title="ReLU"); p=plot!(x,relu.(x),width=4,color="black",legend=false,style=:solid); savefig(p,"relu.png")


text:
export MtlLibrary, MtlLibraryFromFile, MtlLibraryFromData

const MTLLibrary = Ptr{MtLibrary}

"""
    MtlDevice(i::Integer)

Get a handle to a compute device.
"""
mutable struct MtlLibrary
    handle::MTLLibrary
    device::MtlDevice
end

Base.unsafe_convert(::Type{MTLLibrary}, lib::MtlLibrary) = lib.handle

Base.:(==)(a::MtlLibrary, b::MtlLibrary) = a.handle == b.handle
Base.hash(lib::MtlLibrary, h::UInt) = hash(lib.handle, h)

function MtlLibrary(device::MtlDevice, src::String, opts::MtlCompileOptions=MtlCompileOptions())
    handle = @mtlthrows _errptr mtNewLibraryWithSource(device, src, opts, _errptr)

    obj = MtlLibrary(handle, device)
    finalizer(unsafe_destroy!, obj)

    return obj
end

function MtlLibraryFromFile(device::MtlDevice, path::String)
    handle = @mtlthrows _errptr mtNewLibraryWithFile(device, path, _errptr)

    obj = MtlLibrary(handle, device)
    finalizer(unsafe_destroy!, obj)

    return obj
end

function MtlLibraryFromData(device::MtlDevice, data)
    GC.@preserve data begin
        handle = @mtlthrows _errptr mtNewLibraryWithData(device, pointer(data), sizeof(data), _errptr)
    end

    obj = MtlLibrary(handle, device)
    finalizer(unsafe_destroy!, obj)

    return obj
end

function unsafe_destroy!(lib::MtlLibrary)
    mtRelease(lib.handle)
end

## properties

Base.propertynames(::MtlLibrary) = (:device, :label, :functionNames)

function Base.getproperty(lib::MtlLibrary, f::Symbol)
    if f === :label
        ptr = mtLibraryLabel(lib)
        ptr == C_NULL ? nothing : unsafe_string(ptr)
    elseif f === :functionNames
        count = Ref{Csize_t}(0)
        mtLibraryFunctionNames(lib, count, C_NULL)
        names = Vector{Cstring}(undef, count[])
        mtLibraryFunctionNames(lib, count, names)
        unsafe_string.(names)
    else
        getfield(lib, f)
    end
end

## display

function Base.show(io::IO, lib::MtlLibrary)
    print(io, "MtlLibrary($(lib.device))")
end

function Base.show(io::IO, ::MIME"text/plain", lib::MtlLibrary)
    println(io, "MtlLibrary:")
    println(io, " device: ", lib.device)
    print(io,   " label:  ", lib.label)
end


text:
export isbase32

function isbase32(str::AbstractString)::Bool
    base32Reg = r"^[A-Z2-7]+=*$"

    if length(str) % 8 === 0 && contains(str, base32Reg)
        return true
    end

    return false
end


text:
type BSONArray
    _wrap_::Ptr{Void}
    _ref_::Any

    BSONArray() = begin
        _wrap_ = ccall(
            (:bson_new, libbson),
            Ptr{Void}, ()
            )
        bsonArray = new(_wrap_, Union{})
        finalizer(bsonArray, destroy)
        return bsonArray
    end

    BSONArray(vec::Vector) = begin
        bsonArray = BSONArray()
        for ele in vec
            append(bsonArray, ele)
        end
        return bsonArray
    end

    BSONArray(data::Ptr{UInt8}, length::Integer, _ref_::Any) = begin
        buffer = Array(UInt8, 128)
        ccall(
            (:bson_init_static, libbson),
            Bool, (Ptr{Void}, Ptr{UInt8}, UInt32),
            buffer, data, length
            ) || error("bson_init_static: failure")
        b = Compat.unsafe_convert(Ptr{Void}, buffer)
        new(b, (_ref_, buffer))
    end

    BSONArray(_wrap_::Ptr{Void}, _ref_::Any) = new(_wrap_, _ref_)
end
export BSONArray

function convert(::Type{AbstractString}, bsonArray::BSONArray)
    cstr = ccall(
        (:bson_array_as_json, libbson),
        Ptr{UInt8}, (Ptr{Void}, Ptr{UInt8}),
        bsonArray._wrap_,
        C_NULL
        )
    cstr == C_NULL && return ""
    result = unsafe_string(cstr)
    ccall(
        (:bson_free, libbson),
        Void, (Ptr{Void},),
        cstr
        )
    return result
end
export convert

string(bsonArray::BSONArray) = convert(AbstractString, bsonArray)

show(io::IO, bsonArray::BSONArray) = print(io, "BSONArray($(convert(AbstractString, bsonArray)))")
export show

length(bsonArray::BSONArray) =
    ccall(
        (:bson_count_keys, libbson),
        UInt32, (Ptr{Void},),
        bsonArray._wrap_
        )

function append(bsonArray::BSONArray, val::Bool)
    keyCStr = string(length(bsonArray))
    ccall(
        (:bson_append_bool, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint, Bool),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr),
        val
        ) || error("libBSON: overflow")
end
function append(bsonArray::BSONArray, val::Real)
    keyCStr = string(length(bsonArray))
    ccall(
        (:bson_append_double, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint, Cdouble),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr),
        val
        ) || error("libBSON: overflow")
end
using Base.Dates: datetime2unix
function append(bsonArray::BSONArray, val::Union{Date,DateTime})
    keyCStr = string(length(bsonArray))
    ts = (typeof(val) == Date ? datetime2unix(DateTime(val)) : datetime2unix(val)) * 1000
    ccall(
        (:bson_append_date_time, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint, Clong),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr),
        ts
        ) || error("libBSON: overflow")
end
function append(bsonArray::BSONArray, val::BSONArray)
    keyCStr = string(length(bsonArray))
    childBuffer = Array(UInt8, 128)
    ccall(
        (:bson_append_array_begin, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint, Ptr{Void}),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr),
        childBuffer
        ) || error("bson_append_array_begin: failure")
    childBSONArray = BSONArray(Compat.unsafe_convert(Ptr{Void}, childBuffer), childBuffer)
    for element in val
        append(childBSONArray, element)
    end
    ccall(
        (:bson_append_array_end, libbson),
        Bool, (Ptr{Void}, Ptr{Void}),
        bsonArray._wrap_,
        childBuffer
        ) || error("bson_append_array_end: failure")
end
function append(bsonArray::BSONArray, val::Union{Int8, UInt8, Int16, UInt16, Int32, UInt32})
    keyCStr = string(length(bsonArray))
    ccall(
        (:bson_append_int32, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint, Int32),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr),
        val
        ) || error("libBSON: overflow")
end
function append(bsonArray::BSONArray, val::Union{Int64, UInt64})
    keyCStr = string(length(bsonArray))
    ccall(
        (:bson_append_int64, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint, Int64),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr),
        val
        ) || error("libBSON: overflow")
end
function append(bsonArray::BSONArray, val::BSONOID)
    keyCStr = string(length(bsonArray))
    ccall(
        (:bson_append_oid, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint, Ptr{UInt8}),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr),
        val._wrap_
        ) || error("libBSON: overflow")
end
function append(bsonArray::BSONArray, val::AbstractString)
    keyCStr = string(length(bsonArray))
    valUTF8 = String(val)
    ccall(
        (:bson_append_utf8, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint, Ptr{UInt8}, Cint),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr),
        valUTF8,
        sizeof(valUTF8)
        ) || error("libBSON: overflow")
end
function append(bsonArray::BSONArray, val::Void)
    append_null(bsonArray)
end
function append(bsonArray::BSONArray, val::Symbol)
    if val == :null
        append_null(bsonArray)
    elseif val == :minkey
        append_minkey(bsonArray)
    elseif val == :maxkey
        append_maxkey(bsonArray)
    else
        append(bsonArray, string(val))
    end
end
function append(bsonArray::BSONArray, val::Associative)
    keyCStr = string(length(bsonArray))
    childBuffer = Array(UInt8, 128)
    ccall(
        (:bson_append_document_begin, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint, Ptr{Void}),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr),
        childBuffer
        ) || error("bson_append_document_begin: failure")
    childBSONObject = BSONObject(Compat.unsafe_convert(Ptr{Void}, childBuffer), childBuffer)
    for (k, v) in val
        append(childBSONObject, k, v)
    end
    ccall(
        (:bson_append_document_end, libbson),
        Bool, (Ptr{Void}, Ptr{Void}),
        bsonArray._wrap_,
        childBuffer
        ) || error("bson_append_document_end: failure")
end
function append(bsonArray::BSONArray, val::Vector)
    keyCStr = string(length(bsonArray))
    childBuffer = Array(UInt8, 128)
    ccall(
        (:bson_append_array_begin, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint, Ptr{Void}),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr),
        childBuffer
        ) || error("bson_append_array_begin: failure")
    childBSONArray = BSONArray(Compat.unsafe_convert(Ptr{Void}, childBuffer), childBuffer)
    for element in val
        append(childBSONArray, element)
    end
    ccall(
        (:bson_append_array_end, libbson),
        Bool, (Ptr{Void}, Ptr{Void}),
        bsonArray._wrap_,
        childBuffer
        ) || error("bson_append_array_end: failure")
end
export append

function append_null(bsonArray::BSONArray)
    keyCStr = string(length(bsonArray))
    ccall(
        (:bson_append_null, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr)
        ) || error("libBSON: overflow")
end
export append_null

function append_minkey(bsonArray::BSONArray)
    keyCStr = string(length(bsonArray))
    ccall(
        (:bson_append_minkey, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr)
        ) || error("libBSON: overflow")
end
export append_minkey

function append_maxkey(bsonArray::BSONArray)
    keyCStr = string(length(bsonArray))
    ccall(
        (:bson_append_maxkey, libbson),
        Bool, (Ptr{Void}, Ptr{UInt8}, Cint),
        bsonArray._wrap_,
        keyCStr,
        length(keyCStr)
        ) || error("libBSON: overflow")
end
export append_maxkey

vector(bsonArray::BSONArray) = begin
    a = Any[]
    for v in bsonArray
        if isa(v, BSONObject)
            push!(a, dict(v))
        elseif isa(v, BSONArray)
            push!(a, vector(v))
        else
            push!(a, v)
        end
    end
    return a
end
export vector

# Private

function destroy(bsonArray::BSONArray)
    ccall(
        (:bson_destroy, libbson),
        Void, (Ptr{Void},),
        bsonArray._wrap_
        )
end


text:
# Julia wrapper for header: ida.h
# Automatically generated using Clang.jl

function IDACreate()
    ccall((:IDACreate, libsundials_ida), IDAMemPtr, ())
end

function IDAInit(ida_mem, res::IDAResFn, t0::realtype, yy0::N_Vector, yp0::N_Vector)
    ccall((:IDAInit, libsundials_ida), Cint, (IDAMemPtr, IDAResFn, realtype, N_Vector, N_Vector), ida_mem, res, t0, yy0, yp0)
end

function IDAInit(ida_mem, res::IDAResFn, t0::realtype, yy0, yp0)
    __yy0 = convert(NVector, yy0)
    __yp0 = convert(NVector, yp0)
    IDAInit(ida_mem, res, t0, convert(N_vector, __yy0), convert(N_vector, __yp0))
end

function IDASStolerances(ida_mem, reltol::realtype, abstol::realtype)
    ccall((:IDASStolerances, libsundials_ida), Cint, (IDAMemPtr, realtype, realtype), ida_mem, reltol, abstol)
end

function IDASVtolerances(ida_mem, reltol::realtype, abstol::N_vector)
    ccall((:IDASVtolerances, libsundials_ida), Cint, (IDAMemPtr, realtype, N_Vector), ida_mem, reltol, abstol)
end

function IDASVtolerances(ida_mem, reltol::realtype, abstol)
    __abstol = Convert(NVector, abstol)
    IDASVtolerances(ida_mem, reltol::realtype, convert(N_vector, __abstol))
end

function IDAWFtolerances(ida_mem, efun::IDAEwtFn)
    ccall((:IDAWFtolerances, libsundials_ida), Cint, (IDAMemPtr, IDAEwtFn), ida_mem, efun)
end

function IDACalcIC(ida_mem, icopt::Cint, tout1::realtype)
    ccall((:IDACalcIC, libsundials_ida), Cint, (IDAMemPtr, Cint, realtype), ida_mem, icopt, tout1)
end

function IDACalcIC(ida_mem, icopt, tout1::realtype)
    IDACalcIC(ida_mem, convert(Cint, icopt), tout1)
end

function IDASetNonlinConvCoefIC(ida_mem, epiccon::realtype)
    ccall((:IDASetNonlinConvCoefIC, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, epiccon)
end

function IDASetMaxNumStepsIC(ida_mem, maxnh::Cint)
    ccall((:IDASetMaxNumStepsIC, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxnh)
end

function IDASetMaxNumStepsIC(ida_mem, maxnh)
    IDASetMaxNumStepsIC(ida_mem, convert(Cint, maxnh))
end

function IDASetMaxNumJacsIC(ida_mem, maxnj::Cint)
    ccall((:IDASetMaxNumJacsIC, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxnj)
end

function IDASetMaxNumJacsIC(ida_mem, maxnj)
    IDASetMaxNumJacsIC(ida_mem, convert(Cint, maxnj))
end

function IDASetMaxNumItersIC(ida_mem, maxnit::Cint)
    ccall((:IDASetMaxNumItersIC, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxnit)
end

function IDASetMaxNumItersIC(ida_mem, maxnit)
    IDASetMaxNumItersIC(ida_mem, convert(Cint, maxnit))
end

function IDASetLineSearchOffIC(ida_mem, lsoff::Cint)
    ccall((:IDASetLineSearchOffIC, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, lsoff)
end

function IDASetLineSearchOffIC(ida_mem, lsoff)
    IDASetLineSearchOffIC(ida_mem, convert(Cint, lsoff))
end

function IDASetStepToleranceIC(ida_mem, steptol::realtype)
    ccall((:IDASetStepToleranceIC, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, steptol)
end

function IDASetMaxBacksIC(ida_mem, maxbacks::Cint)
    ccall((:IDASetMaxBacksIC, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxbacks)
end

function IDASetMaxBacksIC(ida_mem, maxbacks)
    IDASetMaxBacksIC(ida_mem, convert(Cint, maxbacks))
end

function IDASetErrHandlerFn(ida_mem, ehfun, eh_data)
    ccall((:IDASetErrHandlerFn, libsundials_ida), Cint, (IDAMemPtr, IDAErrHandlerFn, Ptr{Cvoid}), ida_mem, ehfun, eh_data)
end

function IDASetErrFile(ida_mem, errfp)
    ccall((:IDASetErrFile, libsundials_ida), Cint, (IDAMemPtr, Ptr{FILE}), ida_mem, errfp)
end

function IDASetUserData(ida_mem, user_data)
    ccall((:IDASetUserData, libsundials_ida), Cint, (IDAMemPtr, Any), ida_mem, user_data)
end

function IDASetMaxOrd(ida_mem, maxord::Cint)
    ccall((:IDASetMaxOrd, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxord)
end

function IDASetMaxOrd(ida_mem, maxord)
    IDASetMaxOrd(ida_mem, convert(Cint, maxord))
end

function IDASetMaxNumSteps(ida_mem, mxsteps::Clong)
    ccall((:IDASetMaxNumSteps, libsundials_ida), Cint, (IDAMemPtr, Clong), ida_mem, mxsteps)
end

function IDASetMaxNumSteps(ida_mem, mxsteps)
    IDASetMaxNumSteps(ida_mem, convert(Cint, mxsteps))
end

function IDASetInitStep(ida_mem, hin::realtype)
    ccall((:IDASetInitStep, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, hin)
end

function IDASetMaxStep(ida_mem, hmax::realtype)
    ccall((:IDASetMaxStep, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, hmax)
end

function IDASetStopTime(ida_mem, tstop::realtype)
    ccall((:IDASetStopTime, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, tstop)
end

function IDASetNonlinConvCoef(ida_mem, epcon::realtype)
    ccall((:IDASetNonlinConvCoef, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, epcon)
end

function IDASetMaxErrTestFails(ida_mem, maxnef::Cint)
    ccall((:IDASetMaxErrTestFails, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxnef)
end

function IDASetMaxErrTestFails(ida_mem, maxnef)
    IDASetMaxErrTestFails(ida_mem, convert(Cint, maxnef))
end

function IDASetMaxNonlinIters(ida_mem, maxcor::Cint)
    ccall((:IDASetMaxNonlinIters, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxcor)
end

function IDASetMaxNonlinIters(ida_mem, maxcor)
    IDASetMaxNonlinIters(ida_mem, convert(Cint,maxcor))
end

function IDASetMaxConvFails(ida_mem, maxncf::Cint)
    ccall((:IDASetMaxConvFails, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxncf)
end

function IDASetMaxConvFails(ida_mem, maxncf)
    IDASetMaxConvFails(ida_mem, convert(Cint, maxncf))
end

function IDASetSuppressAlg(ida_mem, suppressalg::Cint)
    ccall((:IDASetSuppressAlg, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, suppressalg)
end

function IDASetSuppressAlg(ida_mem, suppressalg)
    IDASetSuppressAlg(ida_mem, convert(Cint, suppressalg))
end

function IDASetId(ida_mem, id::N_Vector)
    ccall((:IDASetId, libsundials_ida), Cint, (IDAMemPtr, N_Vector), ida_mem, id)
end

function IDASetId(ida_mem, id)
    __id = convert(NVector, id)
    IDASetId(ida_mem, convert(N_Vector, __id))
end

function IDASetConstraints(ida_mem, constraints::N_Vector)
    ccall((:IDASetConstraints, libsundials_ida), Cint, (IDAMemPtr, N_Vector), ida_mem, constraints)
end

function IDASetConstraints(ida_mem, constraints)
    __constraints = convert(NVector, constraints)
    IDASetConstraints(ida_mem, convert(N_Vector, constraints))
end

function IDASetNonlinearSolver(ida_mem, NLS::SUNNonlinearSolver)
    ccall((:IDASetNonlinearSolver, libsundials_ida), Cint, (IDAMemPtr, SUNNonlinearSolver), ida_mem, NLS)
end

function IDARootInit(ida_mem, nrtfn::Cint, g::IDARootFn)
    ccall((:IDARootInit, libsundials_ida), Cint, (IDAMemPtr, Cint, IDARootFn), ida_mem, nrtfn, g)
end

function IDARootInit(ida_mem, nrtfn, g::IDARootFn)
    IDARootInit(ida_mem, convert(Cint, nrtfn), g)
end

function IDASetRootDirection(ida_mem, rootdir)
    ccall((:IDASetRootDirection, libsundials_ida), Cint, (IDAMemPtr, Ptr{Cint}), ida_mem, rootdir)
end

function IDASetNoInactiveRootWarn(ida_mem)
    ccall((:IDASetNoInactiveRootWarn, libsundials_ida), Cint, (IDAMemPtr,), ida_mem)
end

function IDASolve(ida_mem, tout::realtype, tret, yret::N_Vector, ypret::N_Vector, itask::Cint)
    ccall((:IDASolve, libsundials_ida), Cint, (IDAMemPtr, realtype, Ptr{realtype}, N_Vector, N_Vector, Cint), ida_mem, tout, tret, yret, ypret, itask)
end

function IDASolve(ida_mem, tout::realtype, tret, yret, ypret, itask)
    __yret = convert(NVector, yret)
    __ypret = convert(NVector, ypret)
    IDASolve(ida_mem, tout::realtype, tret, convert(N_Vector, __yret), convert(N_Vector, __ypret), convert(Cint, itask))
end

function IDAComputeY(ida_mem, ycor::N_Vector, y::N_Vector)
    ccall((:IDAComputeY, libsundials_ida), Cint, (IDAMemPtr, N_Vector, N_Vector), ida_mem, ycor, y)
end

function IDAComputeY(ida_mem, ycor, y)
    __ycor = convert(NVector, ycor)
    __y = convert(NVector, y)
    IDAComputeY(ida_mem, convert(N_Vector, __ycor), convert(N_Vector, __y))
end

function IDAComputeYp(ida_mem, ycor::N_Vector, yp::N_Vector)
    ccall((:IDAComputeYp, libsundials_ida), Cint, (IDAMemPtr, N_Vector, N_Vector), ida_mem, ycor, yp)
end

function IDAComputeYp(ida_mem, ycor, yp)
    __ycor = convert(NVector, ycor)
    __yp = convert(NVector, yp)
    IDAComputeY(ida_mem, convert(N_Vector, __ycor), convert(N_Vector, __y))
end

function IDAGetDky(ida_mem, t::realtype, k::Cint, dky::N_Vector)
    ccall((:IDAGetDky, libsundials_idas), Cint, (IDAMemPtr, realtype, Cint, N_Vector), ida_mem, t, k, dky)
end

function IDAGetDky(ida_mem, t, k, dky)
    __dky = convert(NVector, dky)
    IDAGetDky(ida_mem, t, convert(Cint, k), convert(N_Vector, __dky))
end

function IDAGetWorkSpace(ida_mem, lenrw, leniw)
    ccall((:IDAGetWorkSpace, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}), ida_mem, lenrw, leniw)
end

function IDAGetNumSteps(ida_mem, nsteps)
    ccall((:IDAGetNumSteps, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nsteps)
end

function IDAGetNumResEvals(ida_mem, nrevals)
    ccall((:IDAGetNumResEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nrevals)
end

function IDAGetNumLinSolvSetups(ida_mem, nlinsetups)
    ccall((:IDAGetNumLinSolvSetups, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nlinsetups)
end

function IDAGetNumErrTestFails(ida_mem, netfails)
    ccall((:IDAGetNumErrTestFails, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, netfails)
end

function IDAGetNumBacktrackOps(ida_mem, nbacktr)
    ccall((:IDAGetNumBacktrackOps, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nbacktr)
end

function IDAGetConsistentIC(ida_mem, yy0_mod::N_Vector, yp0_mod::N_Vector)
    ccall((:IDAGetConsistentIC, libsundials_idas), Cint, (IDAMemPtr, N_Vector, N_Vector), ida_mem, yy0_mod, yp0_mod)
end

function IDAGetConsistentIC(ida_mem, yy0_mod, yp0_mod)
    __yy0_mod = convert(NVector, yy0_mod)
    __yp0_mod = convert(NVector, yp0_mod)
    IDAGetConsistentIC(ida_mem, convert(N_Vector, __yy0_mod), convert(N_Vector, __yp0_mod))
end

function IDAGetLastOrder(ida_mem, klast)
    ccall((:IDAGetLastOrder, libsundials_ida), Cint, (IDAMemPtr, Ptr{Cint}), ida_mem, klast)
end

function IDAGetCurrentOrder(ida_mem, kcur)
    ccall((:IDAGetCurrentOrder, libsundials_ida), Cint, (IDAMemPtr, Ptr{Cint}), ida_mem, kcur)
end

function IDAGetCurrentCj(ida_mem, cj)
    ccall((:IDAGetCurrentCj, libsundials_ida), Cint, (IDAMemPtr, Ptr{realtype}), ida_mem, cj)
end

function IDAGetCurrentY(ida_mem, ycur)
    ccall((:IDAGetCurrentY, libsundials_ida), Cint, (IDAMemPtr, Ptr{N_Vector}), ida_mem, ycur)
end

function IDAGetCurrentYp(ida_mem, ypcur)
    ccall((:IDAGetCurrentYp, libsundials_ida), Cint, (IDAMemPtr, Ptr{N_Vector}), ida_mem, ypcur)
end

function IDAGetActualInitStep(ida_mem, hinused)
    ccall((:IDAGetActualInitStep, libsundials_ida), Cint, (IDAMemPtr, Ptr{realtype}), ida_mem, hinused)
end

function IDAGetLastStep(ida_mem, hlast)
    ccall((:IDAGetLastStep, libsundials_ida), Cint, (IDAMemPtr, Ptr{realtype}), ida_mem, hlast)
end

function IDAGetCurrentStep(ida_mem, hcur)
    ccall((:IDAGetCurrentStep, libsundials_ida), Cint, (IDAMemPtr, Ptr{realtype}), ida_mem, hcur)
end

function IDAGetCurrentTime(ida_mem, tcur)
    ccall((:IDAGetCurrentTime, libsundials_ida), Cint, (IDAMemPtr, Ptr{realtype}), ida_mem, tcur)
end

function IDAGetTolScaleFactor(ida_mem, tolsfact)
    ccall((:IDAGetTolScaleFactor, libsundials_ida), Cint, (IDAMemPtr, Ptr{realtype}), ida_mem, tolsfact)
end

function IDAGetErrWeights(ida_mem, eweight::N_Vector)
    ccall((:IDAGetErrWeights, libsundials_idas), Cint, (IDAMemPtr, N_Vector), ida_mem, eweight)
end

function IDAGetErrWeights(ida_mem, eweight)
    __eweight = convert(NVector, eweight)
    IDAGetErrWeights(ida_mem, convert(N_Vector, __eweight))
end

function IDAGetEstLocalErrors(ida_mem, ele::N_Vector)
    ccall((:IDAGetEstLocalErrors, libsundials_idas), Cint, (IDAMemPtr, N_Vector), ida_mem, ele)
end

function IDAGetEstLocalErrors(ida_mem, ele)
    __ele = convert(NVector, ele)
    IDAGetEstLocalErrors(ida_mem, convert(N_Vector, __ele))
end

function IDAGetNumGEvals(ida_mem, ngevals)
    ccall((:IDAGetNumGEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, ngevals)
end

function IDAGetRootInfo(ida_mem, rootsfound)
    ccall((:IDAGetRootInfo, libsundials_ida), Cint, (IDAMemPtr, Ptr{Cint}), ida_mem, rootsfound)
end

function IDAGetIntegratorStats(ida_mem, nsteps, nrevals, nlinsetups, netfails, qlast, qcur, hinused, hlast, hcur, tcur)
    ccall((:IDAGetIntegratorStats, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Cint}, Ptr{Cint}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}), ida_mem, nsteps, nrevals, nlinsetups, netfails, qlast, qcur, hinused, hlast, hcur, tcur)
end

function IDAGetNumNonlinSolvIters(ida_mem, nniters)
    ccall((:IDAGetNumNonlinSolvIters, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nniters)
end

function IDAGetNumNonlinSolvConvFails(ida_mem, nncfails)
    ccall((:IDAGetNumNonlinSolvConvFails, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nncfails)
end

function IDAGetNonlinSolvStats(ida_mem, nniters, nncfails)
    ccall((:IDAGetNonlinSolvStats, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}), ida_mem, nniters, nncfails)
end

function IDAGetReturnFlagName(flag::Clong)
    ccall((:IDAGetReturnFlagName, libsundials_ida), Cstring, (Clong,), flag)
end

function IDAGetReturnFlagName(flag)
    IDAGetReturnFlagName(convert(Clong, flag))
end

function IDAFree(ida_mem)
    ccall((:IDAFree, libsundials_ida), Cvoid, (Ref{IDAMemPtr},), ida_mem)
end
# Julia wrapper for header: ida_bbdpre.h
# Automatically generated using Clang.jl

function IDABBDPrecInit(ida_mem, Nlocal::sunindextype, mudq::sunindextype, mldq::sunindextype, mukeep::sunindextype, mlkeep::sunindextype, dq_rel_yy::realtype, Gres::IDABBDLocalFn, Gcomm::IDABBDCommFn)
    ccall((:IDABBDPrecInit, libsundials_ida), Cint, (IDAMemPtr, sunindextype, sunindextype, sunindextype, sunindextype, sunindextype, realtype, IDABBDLocalFn, IDABBDCommFn), ida_mem, Nlocal, mudq, mldq, mukeep, mlkeep, dq_rel_yy, Gres, Gcomm)
end

function IDABBDPrecReInit(ida_mem, mudq::sunindextype, mldq::sunindextype, dq_rel_yy::realtype)
    ccall((:IDABBDPrecReInit, libsundials_ida), Cint, (IDAMemPtr, sunindextype, sunindextype, realtype), ida_mem, mudq, mldq, dq_rel_yy)
end

function IDABBDPrecGetWorkSpace(ida_mem, lenrwBBDP, leniwBBDP)
    ccall((:IDABBDPrecGetWorkSpace, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}), ida_mem, lenrwBBDP, leniwBBDP)
end

function IDABBDPrecGetNumGfnEvals(ida_mem, ngevalsBBDP)
    ccall((:IDABBDPrecGetNumGfnEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, ngevalsBBDP)
end
# Julia wrapper for header: ida_direct.h
# Automatically generated using Clang.jl

function IDADlsSetLinearSolver(ida_mem, LS::SUNLinearSolver, A::SUNMatrix)
    ccall((:IDADlsSetLinearSolver, libsundials_idas), Cint, (IDAMemPtr, SUNLinearSolver, SUNMatrix), ida_mem, LS, A)
end

function IDADlsSetJacFn(ida_mem, jac::IDADlsJacFn)
    ccall((:IDADlsSetJacFn, libsundials_idas), Cint, (IDAMemPtr, IDADlsJacFn), ida_mem, jac)
end

function IDADlsGetWorkSpace(ida_mem, lenrwLS, leniwLS)
    ccall((:IDADlsGetWorkSpace, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}), ida_mem, lenrwLS, leniwLS)
end

function IDADlsGetNumJacEvals(ida_mem, njevals)
    ccall((:IDADlsGetNumJacEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, njevals)
end

function IDADlsGetNumResEvals(ida_mem, nrevalsLS)
    ccall((:IDADlsGetNumResEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nrevalsLS)
end

function IDADlsGetLastFlag(ida_mem, flag)
    ccall((:IDADlsGetLastFlag, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, flag)
end

function IDADlsGetReturnFlagName(flag)
    ccall((:IDADlsGetReturnFlagName, libsundials_ida), Cstring, (Clong,), flag)
end
# Julia wrapper for header: ida_ls.h
# Automatically generated using Clang.jl

function IDASetLinearSolver(ida_mem, LS::SUNLinearSolver, A::SUNMatrix)
    ccall((:IDASetLinearSolver, libsundials_ida), Cint, (IDAMemPtr, SUNLinearSolver, SUNMatrix), ida_mem, LS, A)
end

function IDASetJacFn(ida_mem, jac::IDADlsJacFn)
    ccall((:IDASetJacFn, libsundials_ida), Cint, (IDAMemPtr, IDALsJacFn), ida_mem, jac)
end

function IDASetPreconditioner(ida_mem, pset::IDALsPrecSetupFn, psolve::IDALsPrecSolveFn)
    ccall((:IDASetPreconditioner, libsundials_ida), Cint, (IDAMemPtr, IDALsPrecSetupFn, IDALsPrecSolveFn), ida_mem, pset, psolve)
end

function IDASetJacTimes(ida_mem, jtsetup::IDALsJacTimesSetupFn, jtimes::IDALsJacTimesVecFn)
    ccall((:IDASetJacTimes, libsundials_ida), Cint, (IDAMemPtr, IDALsJacTimesSetupFn, IDALsJacTimesVecFn), ida_mem, jtsetup, jtimes)
end

function IDASetEpsLin(ida_mem, eplifac::realtype)
    ccall((:IDASetEpsLin, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, eplifac)
end

function IDASetIncrementFactor(ida_mem, dqincfac::realtype)
    ccall((:IDASetIncrementFactor, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, dqincfac)
end

function IDAGetLinWorkSpace(ida_mem, lenrwLS, leniwLS)
    ccall((:IDAGetLinWorkSpace, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}), ida_mem, lenrwLS, leniwLS)
end

function IDAGetNumJacEvals(ida_mem, njevals)
    ccall((:IDAGetNumJacEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, njevals)
end

function IDAGetNumPrecEvals(ida_mem, npevals)
    ccall((:IDAGetNumPrecEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, npevals)
end

function IDAGetNumPrecSolves(ida_mem, npsolves)
    ccall((:IDAGetNumPrecSolves, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, npsolves)
end

function IDAGetNumLinIters(ida_mem, nliters)
    ccall((:IDAGetNumLinIters, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nliters)
end

function IDAGetNumLinConvFails(ida_mem, nlcfails)
    ccall((:IDAGetNumLinConvFails, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nlcfails)
end

function IDAGetNumJTSetupEvals(ida_mem, njtsetups)
    ccall((:IDAGetNumJTSetupEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, njtsetups)
end

function IDAGetNumJtimesEvals(ida_mem, njvevals)
    ccall((:IDAGetNumJtimesEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, njvevals)
end

function IDAGetNumLinResEvals(ida_mem, nrevalsLS)
    ccall((:IDAGetNumLinResEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nrevalsLS)
end

function IDAGetLastLinFlag(ida_mem, flag)
    ccall((:IDAGetLastLinFlag, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, flag)
end

function IDAGetLinReturnFlagName(flag::Clong)
    ccall((:IDAGetLinReturnFlagName, libsundials_ida), Cstring, (Clong,), flag)
end

function IDAGetLinReturnFlagName(flag)
    IDAGetLinReturnFlagName(convert(Clong,flag))
end

# Julia wrapper for header: ida_spils.h
# Automatically generated using Clang.jl

function IDASpilsSetLinearSolver(ida_mem, LS::SUNLinearSolver)
    ccall((:IDASpilsSetLinearSolver, libsundials_ida), Cint, (IDAMemPtr, SUNLinearSolver), ida_mem, LS)
end

function IDASpilsSetPreconditioner(ida_mem, pset::IDASpilsPrecSetupFn, psolve::IDASpilsPrecSolveFn)
    ccall((:IDASpilsSetPreconditioner, libsundials_ida), Cint, (IDAMemPtr, IDASpilsPrecSetupFn, IDASpilsPrecSolveFn), ida_mem, pset, psolve)
end

function IDASpilsSetJacTimes(ida_mem, jtsetup::IDASpilsJacTimesSetupFn, jtimes::IDASpilsJacTimesVecFn)
    ccall((:IDASpilsSetJacTimes, libsundials_ida), Cint, (IDAMemPtr, IDASpilsJacTimesSetupFn, IDASpilsJacTimesVecFn), ida_mem, jtsetup, jtimes)
end

function IDASpilsSetEpsLin(ida_mem, eplifac::realtype)
    ccall((:IDASpilsSetEpsLin, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, eplifac)
end

function IDASpilsSetIncrementFactor(ida_mem, dqincfac::realtype)
    ccall((:IDASpilsSetIncrementFactor, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, dqincfac)
end

function IDASpilsGetWorkSpace(ida_mem, lenrwLS, leniwLS)
    ccall((:IDASpilsGetWorkSpace, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}), ida_mem, lenrwLS, leniwLS)
end

function IDASpilsGetNumPrecEvals(ida_mem, npevals)
    ccall((:IDASpilsGetNumPrecEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, npevals)
end

function IDASpilsGetNumPrecSolves(ida_mem, npsolves)
    ccall((:IDASpilsGetNumPrecSolves, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, npsolves)
end

function IDASpilsGetNumLinIters(ida_mem, nliters)
    ccall((:IDASpilsGetNumLinIters, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nliters)
end

function IDASpilsGetNumConvFails(ida_mem, nlcfails)
    ccall((:IDASpilsGetNumConvFails, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nlcfails)
end

function IDASpilsGetNumJTSetupEvals(ida_mem, njtsetups)
    ccall((:IDASpilsGetNumJTSetupEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, njtsetups)
end

function IDASpilsGetNumJtimesEvals(ida_mem, njvevals)
    ccall((:IDASpilsGetNumJtimesEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, njvevals)
end

function IDASpilsGetNumResEvals(ida_mem, nrevalsLS)
    ccall((:IDASpilsGetNumResEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nrevalsLS)
end

function IDASpilsGetLastFlag(ida_mem, flag)
    ccall((:IDASpilsGetLastFlag, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, flag)
end

function IDASpilsGetReturnFlagName(flag::Clong)
    ccall((:IDASpilsGetReturnFlagName, libsundials_ida), Cstring, (Clong,), flag)
end

function IDASpilsGetReturnFlagName(flag)
    IDASpilsGetReturnFlagName(convert(Clong, flag))
end


text:
using PythonCall, Test, Dates, Aqua

Aqua.test_all(PythonCall)

@testset "PythonCall.jl" begin
    @testset "abstract" begin
        include("abstract.jl")
    end
    @testset "concrete" begin
        include("concrete.jl")
    end
    @testset "jlwrap" begin
        include("jlwrap.jl")
    end
end


text:
module FieldOperationsTests

end # module


text:
using Weave
using Flux
using DataFrames, CSV
using Binning
using Random:shuffle!
using BSON: @save, @load
using Plots
include("./trainModel.jl")

filename = "./evalTemplate.jl"

weave(filename; doctype="md2html", out_path=:pwd)
weave(filename; doctype="github", out_path=:pwd)


text:
__precompile__()
module SingularIntegralEquations
    using Base, BandedMatrices, ApproxFun, DualNumbers, RecipesBase, DomainSets,
            LinearAlgebra, Random, SpecialFunctions, LowRankApprox, InteractiveUtils

export cauchy, cauchyintegral, stieltjes, logkernel,
       stieltjesintegral, hilbert, pseudohilbert, pseudocauchy,
       stieltjesjacobimoment, logjacobimoment, singularintegral

import Base: values, getindex, setindex!, *, +, -, ==, <, <=, >,
                >=, /, ^, \, ∪, transpose, convert, Array, Vector, Matrix,
                AbstractVector, AbstractMatrix, AbstractArray, in

import Base.Broadcast: broadcasted, DefaultArrayStyle

import LinearAlgebra: ldiv!, mul!, rank, cond, qr

import ApproxFun
import DomainSets: UnionDomain, TypedEndpointsInterval

import ApproxFun: bandwidths, blockbandwidths, SpaceOperator, bilinearform, linebilinearform,dotu, blocklengths,
                  plan_transform,plan_itransform,transform,itransform,transform!,itransform!,
                  rangespace, domainspace, promotespaces, InterlaceOperator, coefficientmatrix,
                  canonicalspace, domain, space, Space, promotedomainspace, promoterangespace, AnyDomain, CalculusOperator,
                  ConcreteDefiniteIntegral, ConcreteDefiniteLineIntegral,
                  SumSpace,PiecewiseSpace, interlace,Multiplication, VectorSpace, ArraySpace,
                  BandedMatrix,ChebyshevDirichlet,PolynomialSpace,AbstractProductSpace,evaluate,order,
                  UnsetSpace, MultivariateFun, BivariateFun,linesum,complexlength,
                  Fun, ProductFun, LowRankFun, mappoint, JacobiZ, prectype,
                  real, UnivariateSpace, RealUnivariateSpace, setdomain, eps, choosedomainspace, isapproxinteger,
                  ConstantSpace,ReOperator,DirectSumSpace, ArraySpace, ZeroSpace,
                  LowRankPertOperator, LaurentDirichlet, setcanonicaldomain, SubSpace,
                  IntervalCurve, PeriodicCurve, reverseorientation, @wrapper, mobius,
                  defaultgetindex, WeightSpace, pochhammer, spacescompatible, ∞, LowRankMatrix, SubOperator,
                  Block, BlockBandedMatrix, BandedBlockBandedMatrix, DFunction, Infinity,
                  component, ncomponents, factor, nfactors, components, factors, rangetype,
                  VFun, Point, dynamic, pieces, npieces, piece, cfstype, isreal, IntervalOrSegmentDomain,
                  IntervalOrSegment, canonicaldomain

import ApproxFun: testbandedoperator

import DualNumbers: dual

import LowRankApprox: refactorsvd!

export ⁺, ⁻

"""
`Directed` represents a number that is a limit from either left (s=true) or right (s=false)
For functions with branch cuts, it is assumed that the value is on the branch cut,
Therefore not requiring tolerances.  This will naturally give the analytic continuation.
"""
struct Directed{s,T} <: Number
    x::T
    Directed{s,T}(x::T) where {s,T} = new{s,T}(x)
    Directed{s,T}(x::Number) where {s,T} = new{s,T}(T(x))
end

Directed{s}(x) where {s} = Directed{s,eltype(x)}(x)

convert(::Type{Directed{s,T}},x::Directed{s}) where {s,T} = Directed{s,T}(T(x.x))
convert(::Type{Directed{s,T}},x::T) where {s,T} = Directed{s,T}(x)
convert(::Type{Directed{s,T}},x::Real) where {s,T} = Directed{s,T}(T(x))
convert(::Type{Directed{s,T}},x::Complex) where {s,T} = Directed{s,T}(T(x))

const ⁺ = Directed{true}(true)
const ⁻ = Directed{false}(true)

orientationsign(::Type{Directed{true}}) = 1
orientationsign(::Type{Directed{false}}) = -1
orientation(::Type{Directed{s}}) where {s} = s
orientation(::Directed{s}) where {s} = s

# removes direction from a number
undirected(x::Number) = x
undirected(x::Fun) = x
undirected(x::Directed) = undirected(x.x)  # x might also have directeion
reverseorientation(x::Number) = x
reverseorientation(x::Directed{s}) where {s} = Directed{!s}(reverseorientation(x.x))

for OP in (:*,:+,:-,:/)
    @eval begin
        $OP(a::Directed{s}) where {s} = Directed{s}($OP(a.x))
        $OP(a::Directed{s},b::Directed{s}) where {s} = Directed{s}($OP(a.x,b.x))
        $OP(a::Directed{s},b::Number) where {s} = Directed{s}($OP(a.x,b))
        $OP(a::Number,b::Directed{s}) where {s} = Directed{s}($OP(a,b.x))
    end
end

real(::Type{Directed{s,T}}) where {s,T} = real(T)

# abs, real and imag delete orientation.
for OP in (:(Base.isfinite),:(Base.isinf),:(Base.abs),:(Base.real),:(Base.imag),:(Base.angle))
    @eval $OP(a::Directed) = $OP(a.x)
end

# branchcuts of log, sqrt, etc. are oriented from (0,-∞)
Base.log(x::Directed{true}) = log(-x.x) - π*im
Base.log(x::Directed{false}) = log(-x.x) + π*im
Base.log1p(x::Directed) = log(1+x)
Base.sqrt(x::Directed{true}) = real(x.x) ≥ 0 ? sqrt(complex(x.x)) : -im*sqrt(-x.x)
Base.sqrt(x::Directed{false}) = real(x.x) ≥ 0 ? sqrt(complex(x.x)) : im*sqrt(-x.x)
^(x::Directed{true},a::Integer) = x.x^a
^(x::Directed{false},a::Integer) = x.x^a
^(x::Directed{true},a::Number) = exp(-a*π*im)*(-x.x)^a
^(x::Directed{false},a::Number) = exp(a*π*im)*(-x.x)^a

dual(a::Directed{s},b) where {s} = Directed{s}(dual(undirected(a),b))

for OP in (:stieltjes,:stieltjesintegral,:pseudostieltjes)
    @eval $OP(f::Fun) = $OP(space(f),coefficients(f))
    @eval $OP(f::Fun,z) = $OP(space(f),coefficients(f),z)
end

stieltjes(s::SumSpace, c::AbstractVector, z) =
    sum(stieltjes.(components(Fun(s,c)), z))

hilbert(f) = Hilbert()*f
hilbert(S,f,z) = hilbert(Fun(S,f))(z)
hilbert(f::Fun,z) = hilbert(space(f),coefficients(f),z)

logkernel(f::Fun,z) = logkernel(space(f),coefficients(f),z)
in(x::Directed, d::Domain) = x.x ∈ d
in(x::Directed, d::TypedEndpointsInterval{:closed,:closed}) = x.x ∈ d

cauchy(f...) = stieltjes(f...)*(im/(2π))
pseudocauchy(f...) = pseudostieltjes(f...)*(im/(2π))
cauchyintegral(u...) = stieltjesintegral(u...)*(im/(2π))

function singularintegral(k::Integer,s::Space,f,z)
    k == 0 && return logkernel(s,f,z)
    k == 1 && isreal(domain(s)) && return stieltjes(s,f,z)/π
    k == 1 && domain(s) isa Segment && return stieltjes(s,f/sign(domain(s)),z)/π
    error("Not implemented")
end

singularintegral(k::Integer,s::PiecewiseSpace,f,z) =
    mapreduce(g -> singularintegral(k, g, z), +, pieces(Fun(s,f)))

function csingularintegral(k::Integer,s::Space,f,z)
    k == 0 && error("Not defined")
    k == 1 && return stieltjes(s,f,z)/π
    return epsilon(csingularintegral(k-1,s,f,dual(z,1)))
end

singularintegral(k::Integer,f::Fun,z) = singularintegral(k,space(f),coefficients(f),z)
csingularintegral(k::Integer,f::Fun,z) = csingularintegral(k,space(f),coefficients(f),z)

# Modifier spaces

stieltjes(sp::SubSpace,v,z) = stieltjes(sp.space,coefficients(v,sp,sp.space),z)

include("Operators/Operators.jl")
include("FundamentalSolutions/FundamentalSolutions.jl")
include("HypergeometricFunctions/HypergeometricFunctions.jl")

include("JacobiQWeight.jl")
include("JacobiQ.jl")
include("stieltjesmoment.jl")

include("circlecauchy.jl")
include("intervalcauchy.jl")
include("singfuncauchy.jl")

include("vectorcauchy.jl")

include("GreensFun/GreensFun.jl")

include("periodicline.jl")
include("arc.jl")
include("curve.jl")
include("asymptotics.jl")

# if isdir(Pkg.dir("TikzGraphs"))
#     include("introspect.jl")
# end

include("Extras/Extras.jl")

using Test

function testsieoperators(S::Space)
    testbandedoperator(SingularIntegral(S,0))
    testbandedoperator(SingularIntegral(S,1))
    testbandedoperator(Hilbert(S))
    p=ApproxFun.checkpoints(S)[1] # random point on contour
    x=Fun(domain(S))
    z=2.12312231+1.433453443534im # random point not on contour

    for k=1:5
        f=Fun(S,[zeros(k-1);1])
        @test (SingularIntegral(S,0)*f)(p) ≈ logkernel(f,p)
        @test (Hilbert(S,1)*f)(p) ≈ hilbert(f,p)
    end
end

function testsieeval(S::Space;posdirection=im)
    p=ApproxFun.checkpoints(S)[1] # random point on contour
    x=Fun(domain(S))
    z=2.12312231+1.433453443534im # random point not on contour

    for k=1:5
        f=Fun(S,[zeros(k-1);1])
        @test abs(sum(f/(z-x))-stieltjes(f,z)) ≤ 100eps()
        @test stieltjes(f,(p)⁺) ≈ stieltjes(f,p+eps()*posdirection)
        @test stieltjes(f,(p)⁻) ≈ stieltjes(f,p-eps()*posdirection)
        @test cauchy(f,(p)⁺)-cauchy(f,(p)⁻) ≈ f(p)
        @test im*(cauchy(f,(p)⁺)+cauchy(f,(p)⁻)) ≈ hilbert(f,p)

        @test abs(linesum(f*logabs(x-z))/π-logkernel(f,z)) ≤ 100eps()
        @test logkernel(f,p) ≈ logkernel(f,p+eps()*posdirection)
    end
end

function testsies(S::Space;posdirection=im)
    testsieoperators(S)
    testsieeval(S;posdirection=posdirection)
end

function testsies(d::IntervalOrSegmentDomain; posdirection=im)
    testsies(JacobiWeight(-0.5,-0.5,Chebyshev(d)))
    testsies(JacobiWeight(0.5,0.5,Ultraspherical(1,d)))
    testsies(JacobiWeight(-0.5,-0.5,ChebyshevDirichlet{1,1}(d)))
    testsieeval(Legendre(d);posdirection=posdirection)
end

end #module


text:
using QContinuousMeasurement
using HDF5
using DelimitedFiles
using TimerOutputs
using LinearAlgebra

function random_rholike(N)
    blocksizes = QContinuousMeasurement.block_sizes(N)
    blocks = map(x -> x + x', map(x-> rand(ComplexF64, x, x), blocksizes))
    r = BlockDiagonal(blocks)
    return r / tr(r)
end

outsteps = 1
Tfinal = 1.0
dt = 0.01
outpoints = 10
Ntraj = 3

if outpoints > 0
    try
        outsteps = Int(round(Tfinal / dt / outpoints, digits=3))
    catch InexactError
        @warn "The requested $outpoints output points does not divide
        the total time steps. Using the full time output."
    end
end

outsteps = Int(round(Tfinal / dt / outpoints, digits=3))

modelparams = ModelParameters(Nj=40, kind=1.0, Gamma=1.0, omega=0.01, Tfinal=1.0, dt=dt)

Ntime = Int(floor(Tfinal/dt)) # Number of timesteps
t = (1 : Ntime) * dt
t = t[outsteps:outsteps:end]

# Precompile
model_p = InitializeModel(modelparams)
initial_state_p = blockdiag_css(modelparams.Nj)
dy = measure_current(initial_state_p, model_p)
#updatekraus!(model_p, dy)
tr_ρ, tr_τ = updatestate!(initial_state_p, model_p, dy)
QFI(initial_state_p)
expectation_value!(initial_state_p, model_p.Jx)

TimerOutputs.enable_debug_timings(QContinuousMeasurement)
reset_timer!()
@timeit "Model creation" model = InitializeModel(modelparams)
@timeit "State creation" initial_state = blockdiag_css(modelparams.Nj)

@info "Starting trajectories simulation..."
for ktraj = 1 : Ntraj
    @timeit "State copy" state = State(copy(initial_state.ρ))

    # Output variables
    jx = similar(t)
    jy = similar(t)
    jz = similar(t)

    jx2 = similar(t)
    jy2 = similar(t)
    jz2 = similar(t)

    FisherT = similar(t)
    QFisherT = similar(t)

    jto = 1 # Counter for the output
    for jt = 1 : Ntime
        @timeit "Current" dy = measure_current(state, model)
        #@timeit to "Kraus" updatekraus!(model, dy)
        @timeit "Dynamics" tr_ρ, tr_τ = updatestate!(state, model, dy)

        if jt % outsteps == 0
        @timeit "Output" begin
            jx[jto] = expectation_value!(state, model.Jx)
            jy[jto] = expectation_value!(state, model.Jy)
            jz[jto] = expectation_value!(state, model.Jz)

            jx2[jto] = expectation_value!(state, model.Jx2)
            jy2[jto] = expectation_value!(state, model.Jy2)
            jz2[jto] = expectation_value!(state, model.Jz2)

            # We evaluate the classical FI for the continuous measurement
            FisherT[jto] = real(tr_τ^2)

            # We evaluate the QFI for a final strong measurement done at time t
            @timeit "QFI" QFisherT[jto] = QFI(state)

            jto += 1
            end
        end
    end
end
print_timer()


text:
@testset "stddpp.jl" begin
    test_Zalg(StdDPP(SqExponentialKernel()))
end


text:
using ChebyshevQuantum
using LinearAlgebra

function go(N ; Rmax = 20, L=0)

    @time begin
        pts = ChebyshevQuantum.Interp.getCpts(N);
        D = ChebyshevQuantum.Interp.getD(pts);
    end

    @time begin
        r = (1 .+ pts) * Rmax / 2
        H1 = -0.5*(D*D )[2:N,2:N] * 4/Rmax^2
        H2 = -diagm(1 ./ r[2:N])
        H3 = L*(L+1)*diagm(1 ./ r[2:N].^2) / 2.0
        H = H1+H2+H3
    end

    @time v = eigvals( H )

    m = min(10, length(v))
    return v[1:m] * 2 .* ( (1:m) .+ L).^2
end

function go_exp(N ; Rmax = 20, L=0, a = 1.0)

    @time begin
        pts = ChebyshevQuantum.Interp.getCpts(N);
        D = ChebyshevQuantum.Interp.getD(pts);
    end

    @time begin
        r = (1 .+ pts) * Rmax / 2
        H1 = -0.5*(D*D )[2:N,2:N] * 4/Rmax^2

        #H1a = zeros(N-1,N-1)
        H1a = (D )[2:N,2:N] * 2/Rmax * a

        #        H1b = -0.5*I(N-1) * a^2
        #H1b = zeros(N-1, N-1)

        H2 = -diagm(1 ./ r[2:N])
        H3 = L*(L+1)*diagm(1 ./ r[2:N].^2) / 2.0
        H = H1+H1a+H2+H3
    end

    @time v = eigvals( H )
    v = v .- 0.5 * a^2
    m = min(10, length(v))
    return v[1:m] * 2 .* ( (1:m) .+ L).^2
end

#=function go_exp2(N ; Rmax = 20, L=0, a = 1.0)

    @time begin
        pts = ChebyshevQuantum.Interp.getCpts(N);
        D = ChebyshevQuantum.Interp.getD(pts);
    end

    @time begin
        r = (1 .+ pts) * Rmax / 2
        H1 = -0.5*(D*D )[2:N,2:N] * 4/Rmax^2

        #H1a = zeros(N-1,N-1)
        H1a = (D )[2:N,2:N] * 2/Rmax * a

        H1b = -0.5*I(N-1) * a^2
        #H1b = zeros(N-1, N-1)

        H2 = -diagm(1 ./ r[2:N])
        H3 = L*(L+1)*diagm(1 ./ r[2:N].^2) / 2.0
        H = H1+H1a+H1b+H2+H3
    end

    @time v = eigvals( H )
#    v = v .- 0.5 * a^2
    m = min(10, length(v))
    return v[1:m] * 2 .* ( (1:m) .+ L).^2
end
=#