text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    GraphMol, QueryMol,
    graphmol, querymol,
    SDFile, SMILES, SMARTS

struct GraphMol{A<:Atom,B<:Bond} <: OrderedGraph
    neighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    nodeattrs::Vector{A}
    edgeattrs::Vector{B}
    cache::Dict{Symbol,Any}
    attributes::Dict{Symbol,Any}
end

"""
    graphmol() -> GraphMol

Generate empty `GraphMol`.
"""
graphmol(::Type{A}, ::Type{B}
    ) where {A<:Atom,B<:Bond} = GraphMol{A,B}([], [], [], [], Dict(), Dict())

"""
    graphmol(atoms::Vector{Atom}, bonds::Vector{Bond}) -> GraphMol

Generate `GraphMol` that has the given atom objects and edge objects.
"""
function graphmol(
        edges, atoms::Vector{A}, bonds::Vector{B}) where {A<:Atom,B<:Bond}
    nbrmap = [Dict{Int,Int}() for i in 1:length(atoms)]
    edges = collect(edges)
    for (i, (u, v)) in enumerate(edges)
        nbrmap[u][i] = v
        nbrmap[v][i] = u
    end
    return GraphMol(
        nbrmap, edges, atoms, bonds, Dict{Symbol,Any}(), Dict{Symbol,Any}())
end

"""
    graphmol(mol::GraphMol) -> GraphMol

Copy `GraphMol`.
"""
function graphmol(graph::GraphMol{A,B}) where {A<:Atom,B<:Bond}
    newg = graphmol(A, B)
    for nbr in graph.neighbormap
        push!(newg.neighbormap, copy(nbr))
    end
    append!(newg.edges, graph.edges)
    append!(newg.nodeattrs, graph.nodeattrs)
    append!(newg.edgeattrs, graph.edgeattrs)
    merge!(newg.cache, graph.cache)
    merge!(newg.attributes, graph.attributes)
    return newg
end

"""
    graphmol(mol::SubgraphView{GraphMol}) -> GraphMol

Generate a new `GraphMol` from a substructure view.

Graph property caches and attributes are not inherited.
"""
function graphmol(view::SubgraphView)
    newg = graphmol(nodeattrtype(view), edgeattrtype(view))
    nkeys = sort(collect(nodeset(view)))
    ekeys = sort(collect(edgeset(view)))
    nmap = Dict{Int,Int}()
    for (i, n) in enumerate(nkeys)
        nmap[n] = i
        push!(newg.nodeattrs, nodeattr(view, n))
        push!(newg.neighbormap, Dict())
    end
    for (i, e) in enumerate(ekeys)
        (oldu, oldv) = getedge(view, e)
        u = nmap[oldu]
        v = nmap[oldv]
        push!(newg.edges, (u, v))
        push!(newg.edgeattrs, edgeattr(view, e))
        newg.neighbormap[u][i] = v
        newg.neighbormap[v][i] = u
    end
    return newg
end

struct QueryMol{A<:QueryAtom,B<:QueryBond} <: OrderedGraph
    neighbormap::Vector{Dict{Int,Int}}
    edges::Vector{Tuple{Int,Int}}
    nodeattrs::Vector{A}
    edgeattrs::Vector{B}
    cache::Dict{Symbol,Any}
    attributes::Dict{Symbol,Any}
    connectivity::Vector{Vector{Int}}
end

"""
    querymol() -> QueryMol

Generate empty `QueryMol`.
"""
querymol(
    ::Type{A}, ::Type{B}
) where {A<:QueryAtom,B<:QueryBond} = QueryMol{A,B}(
    [], [], [], [], Dict(), Dict(), [])

"""
    querymol(atoms::Vector{Atom}, bonds::Vector{Bond}) -> GraphMol

Generate `QueryMol` that has the given atom objects and edge objects.
"""
function querymol(edges, atoms::Vector{A}, bonds::Vector{B},
        connectivity::Vector{Vector{Int}}) where {A<:QueryAtom,B<:QueryBond}
    nbrmap = [Dict{Int,Int}() for i in 1:length(atoms)]
    edges = collect(edges)
    for (i, (u, v)) in enumerate(edges)
        nbrmap[u][i] = v
        nbrmap[v][i] = u
    end
    return QueryMol(
        nbrmap, edges, atoms, bonds,
        Dict{Symbol,Any}(), Dict{Symbol,Any}(), connectivity)
end

# Aliases

SDFile = GraphMol{SDFileAtom,SDFileBond}
SMILES = GraphMol{SmilesAtom,SmilesBond}
SMARTS = QueryMol{SmartsAtom,SmartsBond}

function Base.getindex(graph::GraphMol, sym::Symbol)
    if haskey(graph.cache, sym)
        return graph.cache[sym]
    else
        return eval(Expr(:call, sym, graph))
    end
end

function Base.getindex(view::SubgraphView, sym::Symbol)
    if haskey(view.graph.cache, sym)
        return view.graph.cache[sym]
    else
        return eval(Expr(:call, sym, view.graph))
    end
end


text:
# https://levelup.gitconnected.com/functional-one-liners-in-julia-e0ed35d4ff7b

f = <(4) # (::Base.Fix2{typeof(<), Int64}) (generic function with 1 method)
@show f
@show filter(<(4), 1:10)
@show findall(==(4), [4, 8, 4, 2, 1, 5])
@show findfirst(==(4), [4, 8, 4, 2, 1, 5])

nothing

ext: jl


text:
using HDF5, JLD

# writes x to group g under name v, with compression
function writecompressed{T<:Number}(g::HDF5Group, v, x::AbstractArray{T})
    g[v, "chunk", size(x), "compress", 9] = x
end
writecompressed(g::HDF5Group, v, x) = write(g, v, x)

# writes variables vars to group gname in file fname
macro writedata(fname, gname, vars...)
    writeexprs = Array(Expr, length(vars))
    for i in 1:length(vars)
        # remove daata first, if already exists, then write in all cases
        writeexprs[i] = Expr(:block,
            :(if $(string(vars[i])) in names(g); o_delete(g, $(string(vars[i]))); end),
            :(writecompressed(g.plain, $(string(vars[i])), $(esc(vars[i])))))
    end
    # open file "r+" if already exists, "w" otherwise
    Expr(:block,
         :(local f = jldopen($fname, isfile($fname) ? "r+" : "w")),
         Expr(:try,
            Expr(:block,
                :(local g = $gname in names(f) ? f[$gname] : g_create(f, $gname)),
                writeexprs...),
            false, false, :(close(f)))
        )
end

# reads variables vars from group gname in file fname
macro readdata(fname, gname, vars...)
    readexprs = Array(Expr, length(vars))
    for i in 1:length(vars)
        readexprs[i] = :($(esc(vars[i])) = read(g, $(string(vars[i]))))
    end
    Expr(:block,
         :(local f = jldopen($fname, "r")),
         Expr(:try,
            Expr(:block,
                :(local g = f[$gname]),
                Expr(:global, map(esc, vars)...),
                readexprs...),
            false, false, :(close(f)))
        )
end

# computes a histogram and binned statistics of additional variables
#
# v is a vector of variables to be binned into n bins. vars is a list of
# further variables (vectors the same size as v) whose statistics are computed
# for each bin.
#
# The function returns (e, counts, stats...), where e is the range of bin
# edges, counts is a vector over counts of bins, and stats is a sequence of
# n x 2 statistic matrices, which means and SEM per bin in the first and
# and second column, respectively.
function histbin(v, n::Integer, vars...)
    minv, maxv = minimum(v), maximum(v)
    w = (maxv - minv) / n
    nvars, nv = length(vars), length(v)
    counts = zeros(Integer, n)
    stats = [zeros(Float64, n, 2) for i = 1:nvars]
    # collect counts and <X> / <X^2> statistics
    for i = 1:nv
        const b = (v[i] - minv) / w
        const bi = b >= float(n) ? n : ifloor(b) + 1
        counts[bi] += 1
        for j = 1:nvars
            const varsji = vars[j][i]
            stats[j][bi,:] += [varsji (varsji * varsji)]
        end
    end
    # turn statistics into mean and SEM
    for i = 1:n
        const ci = counts[i]
        if ci == 0
            for j = 1:nvars
                stats[j][i,:] = fill(NaN, 2)
            end
        elseif ci == 1
            for j = 1:nvars
                stats[j][i,2] = NaN
            end
        else
            for j = 1:nvars
                stats[j][i,:] /= ci
                stats[j][i,2] = sqrt(max(0.0, stats[j][i,2] - stats[j][i,1]^2) / ci)
            end
        end
    end
    tuple(minv:w:maxv, counts, stats...)
end


text:
# Test that the transition density
# in a non-linear, non-homogenous, non-constant diffusivity model
# estimated by forward simulation
# agrees with the density obtained from the linearisation
# reweighted with importance weights using guided proposals.
const 𝕏 = SVector{1}
using GaussianDistributions

struct TargetSDE <: Bridge.ContinuousTimeProcess{Float64}
end
struct LinearSDE{T}  <: Bridge.ContinuousTimeProcess{Float64}
    σ::T
end
Bridge.b(s, x, P::TargetSDE) = -0.1x .+ .5sin(x[1]) .+ 0.5sin(s/4)
Bridge.b(s, x, P::LinearSDE) = Bridge.B(s, P)*x + Bridge.β(s, P)
Bridge.B(s, P::LinearSDE) = SMatrix{1}(-0.1)
Bridge.β(s, P::LinearSDE) = 𝕏(0.5sin(s/4))

Bridge.σ(s, x, P::TargetSDE) = SMatrix{1}(2.0 + 0.5cos(x[1]))
Bridge.σ(s, x, P::LinearSDE) = P.σ
Bridge.σ(s, P::LinearSDE) = P.σ
Bridge.a(s, P::LinearSDE) = P.σ^2

Bridge.constdiff(::TargetSDE) = false
Bridge.constdiff(::LinearSDE) = true

binind(r, x) = searchsortedfirst(r, x) - 1

_euler = BridgeSDEInference.EulerMaruyamaBounded()

function test_measchange()
    Random.seed!(1)
    fextra = 1.
    f = 1.0
    T = round(f*4*pi, digits=2)
    P = TargetSDE()
    v = 𝕏(pi/2)

    x0 = 𝕏(-pi/2)

    t = 0:0.01*f*fextra:T
    t = Bridge.tofs.(t, 0, T)
    W = Bridge.samplepath(t, 𝕏(0.0))

    Wnr =  Wiener{𝕏{Float64}}()

    Σ = SMatrix{1}(0.1)
    L = SMatrix{1}(1.0)
    Noise = Gaussian(𝕏(0.0), Σ)

    sample!(W, Wnr)
    _, X = solve(_euler, x0, W, P)
    v1 = X.yy[end]
    X.yy[end] = zero(v1)
    solve!(_euler, X, x0, W, P)
    @test v1 ≈ X.yy[end]

    K = 50
    vrange = range(-10,10, length=K+1)
    vints = [(vrange[i], vrange[i+1]) for i in 1:K]

    k = 1
    N = 50000

    # Forward simulation

    vs = Float64[]
    for i in 1:N
        sample!(W, Wnr)
        solve!(_euler, X, x0, W, P)
        v = L*X.yy[end] + rand(Noise)
        push!(vs, v[1])
    end

    counts = zeros(K+2)
    [counts[binind(vrange, v)+1] += 1 for v in vs]
    counts /= length(vs)

    wcounts = zeros(K)

    VProp = Uniform(-10,10)

    fpt = fill(NaN, 1)
    v = 𝕏(0.0)
    P̃ = LinearSDE(Bridge.σ(T, v, P)) # use a law with large variance
    Pᵒ = BSI.GuidPropBridge(eltype(x0), t, P, P̃, L, v, Σ)

    # Guided proposals

    for i in 1:N
        v = 𝕏(5*rand(VProp))
        while !((binind(vrange, v[1]) in 1:K))
            v = 𝕏(5*rand(VProp))
        end
        # other possibility: change proposal each step
    #    P̃ = LinearSDE(Bridge.σ(T, v, P))
        Pᵒ = BSI.GuidPropBridge(eltype(x0), t, P, P̃, L, v, Σ)

        sample!(W, Wnr)
        solve!(_euler, X, x0, W, Pᵒ)
        ll = BSI.path_log_likhd(BSI.PartObs(), [X], [Pᵒ], 1:1, fpt, skipFPT=true)
        ll += BSI.lobslikelihood(Pᵒ, x0)
        ll -= logpdf(VProp, v[1])
        wcounts[binind(vrange, v[1])] += exp(ll)/N
    end
    bias = wcounts - counts[2:end-1]

    @testset "Statistical correctness of guided proposals" begin
        @test norm(bias) < 0.05
    end
end

test_measchange()


text:
# SPDX-License-Identifier: MIT

"""
    AbstractDirichlet{D} <: Distribution

Abstract type for Dirichlet distribution implementations.
"""
abstract type AbstractDirichlet{D} <: Distribution end

#######################################################################
# Parameter of the Dirichlet distribution.

function DefaultDirichletParameter(α)
    DefaultParameter(α)
end

#######################################################################
# Dirichlet distribution

"""
    mutable struct Dirichlet{D} <: AbstractDirichlet{D}
        param::P where P <: AbstractParameter
    end

Dirichlet distribution.

# Constructors

    Dirichlet(α)

where `α` is a vector of concentrations.

# Examples
```jldoctest
julia> Dirichlet([1.0, 2.0, 3.0])
Dirichlet{3}:
  α = [1.0, 2.0, 3.0]
```
"""
mutable struct Dirichlet{D} <: AbstractDirichlet{D}
    param::P where P <: AbstractParameter
end

function Dirichlet(α::AbstractVector)
    param = DefaultDirichletParameter(α)
    P = typeof(param)
    D = length(α)
    Dirichlet{D}(param)
end

#######################################################################
# Distribution interface.

basemeasure(::AbstractDirichlet, x::AbstractVector) = -log.(x)

function lognorm(d::AbstractDirichlet, η::AbstractVector = naturalform(d.param))
    sum(loggamma_dot(η)) - loggamma(sum(η))
end

function gradlognorm(d::AbstractDirichlet)
    α = stdparam(d).α
    digamma.(α) .- digamma(sum(α))
end

function sample(d::AbstractDirichlet{D}, size) where D
    α = stdparam(d).α
    retval = similar(α, D, size)
    d_ = Dists.Dirichlet(Array(α))
    retval[:, :] = Dists.rand(d_, size)
    retval
end

# vectorization is effect less
splitgrad(d::AbstractDirichlet, μ) = μ

stats(::AbstractDirichlet, x::AbstractVector) = log.(x)

function stdparam(d::AbstractDirichlet, η::AbstractVector = naturalform(d.param))
    (α = η,)
end

max_stars_repo_name: lucasondel/ExpFamilyDistributions-jl

max_issues_repo_name: BUTSpeechFIT/ExpFamilyDistributions

max_forks_repo_name: lucasondel/ExpFamilyDistributions


text:
export Tensor,
       inds,
       store

"""
Tensor{StoreT,IndsT}

A plain old tensor (with order independent
interface and no assumption of labels)
"""
struct Tensor{ElT,N,StoreT,IndsT} <: AbstractArray{ElT,N}
  store::StoreT
  inds::IndsT
  # The resulting Tensor is a view into the input data
  function Tensor(store::StoreT,inds::IndsT) where {StoreT,IndsT}
    new{eltype(StoreT),length(IndsT),StoreT,IndsT}(store,inds)
  end
end

store(T::Tensor) = T.store
inds(T::Tensor) = T.inds
ind(T::Tensor,j::Integer) = inds(T)[j]

#
# Tools for working with Dims/Tuples
#

# dim and dims are used in the Tensor interface, overload
# base Dims here
dims(ds::Dims) = ds
dim(ds::Dims) = prod(ds)

Base.length(ds::Type{<:Dims{N}}) where {N} = N

# This may be a bad idea to overload?
# Type piracy?
Base.strides(is::Dims) = Base.size_to_strides(1, dims(is)...)
Base.copy(ds::Dims) = ds

## TODO: should this be StaticArrays.similar_type?
#Base.promote_rule(::Type{<:Dims},
#                  ::Type{Val{N}}) where {N} = Dims{N}

ValLength(::Type{Dims{N}}) where {N} = Val{N}
ValLength(::Dims{N}) where {N} = Val{N}()

# This is to help with some generic programming in the Tensor
# code (it helps to construct a Tuple(::NTuple{N,Int}) where the
# only known thing for dispatch is a concrete type such
# as Dims{4})
StaticArrays.similar_type(::Type{<:Dims},
                          ::Type{Val{N}}) where {N} = Dims{N}

unioninds(is1::Dims{N1},
          is2::Dims{N2}) where {N1,N2} = Dims{N1+N2}((is1...,is2...))

function deleteat(t::NTuple{N},pos::Int) where {N}
  return ntuple(i -> i < pos ? t[i] : t[i+1],Val(N-1))
end

function insertat(t::NTuple{N},
                  val::NTuple{M},
                  pos::Int) where {N,M}
  return ntuple(i -> i < pos ? t[i] :
                ( i > pos+M-1 ? t[i-1] :
                 val[i-pos+1] ), Val(N+M-1))
end

#
# Generic Tensor functions
#

# The size is obtained from the indices
dims(T::Tensor) = dims(inds(T))
dim(T::Tensor) = dim(inds(T))
Base.size(T::Tensor) = dims(T)

Base.copy(T::Tensor) = Tensor(copy(store(T)),copy(inds(T)))

Base.complex(T::Tensor) = Tensor(complex(store(T)),copy(inds(T)))

function Base.similar(::Type{T},dims) where {T<:Tensor{<:Any,<:Any,StoreT}} where {StoreT}
  return Tensor(similar(StoreT,dims),dims)
end
# TODO: make sure these are implemented correctly
#Base.similar(T::Type{<:Tensor},::Type{S}) where {S} = Tensor(similar(store(T),S),inds(T))
#Base.similar(T::Type{<:Tensor},::Type{S},dims) where {S} = Tensor(similar(store(T),S),dims)

Base.similar(T::Tensor) = Tensor(similar(store(T)),copy(inds(T)))
Base.similar(T::Tensor,dims) = Tensor(similar(store(T),dims),dims)
# To handle method ambiguity with AbstractArray
Base.similar(T::Tensor,dims::Dims) = Tensor(similar(store(T),dims),dims)
Base.similar(T::Tensor,::Type{S}) where {S} = Tensor(similar(store(T),S),copy(inds(T)))
Base.similar(T::Tensor,::Type{S},dims) where {S} = Tensor(similar(store(T),S),dims)
# To handle method ambiguity with AbstractArray
Base.similar(T::Tensor,::Type{S},dims::Dims) where {S} = Tensor(similar(store(T),S),dims)

#function Base.convert(::Type{Tensor{<:Number,N,StoreR,Inds}},
#                      T::Tensor{<:Number,N,<:Any,Inds}) where {N,Inds,StoreR}
#  return Tensor(convert(StoreR,store(T)),copy(inds(T)))
#end

function Base.zeros(::Type{<:Tensor{ElT,<:Any,StoreT,<:Any}},inds::IndsT) where {ElT,N,StoreT,IndsT}
  return Tensor(zeros(StoreT,dim(inds)),inds)
end

function Base.promote_rule(::Type{<:Tensor{ElT1,N1,StoreT1}},
                           ::Type{<:Tensor{ElT2,N2,StoreT2}}) where {ElT1,ElT2,N1,N2,StoreT1,StoreT2}
  return Tensor{promote_type(ElT1,ElT2),N3,promote_type(StoreT1,StoreT2)} where {N3}
end

function Base.promote_rule(::Type{Tensor{ElT1,N,StoreT1,Inds}},
                           ::Type{Tensor{ElT2,N,StoreT2,Inds}}) where {ElT1,ElT2,N,
                                                                       StoreT1,StoreT2,Inds}
  return Tensor{promote_type(ElT1,ElT2),N,promote_type(StoreT1,StoreT2),Inds}
end

#function Base.promote_rule(::Type{<:Tensor{ElT,<:Any,StoreT,<:Any}},::Type{IndsR}) where {N,ElT,StoreT,IndsR}
#  return Tensor{ElT,length(IndsR),StoreT,IndsR}
#end

function StaticArrays.similar_type(::Type{<:Tensor{ElT,<:Any,StoreT,<:Any}},indsR) where {N,ElT,StoreT}
  return Tensor{ElT,length(indsR),StoreT,typeof(indsR)}
end

Base.BroadcastStyle(::Type{T}) where {T<:Tensor} = Broadcast.ArrayStyle{T}()

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{T}},
                      ::Type{<:Any}) where {T<:Tensor}
  A = find_tensor(bc)
  return similar(A)
end

# This is used for overloading broadcast
"`A = find_tensor(As)` returns the first Tensor among the arguments."
find_tensor(bc::Base.Broadcast.Broadcasted) = find_tensor(bc.args)
find_tensor(args::Tuple) = find_tensor(find_tensor(args[1]), Base.tail(args))
find_tensor(x) = x
find_tensor(a::Tensor, rest) = a
find_tensor(::Any, rest) = find_tensor(rest)

# TODO: implement some generic fallbacks for necessary parts of the API?
#Base.getindex(A::TensorT, i::Int) where {TensorT<:Tensor} = error("getindex not yet implemented for Tensor type $TensorT")


text:
module LineDP

import Printf: @printf
import ..Pwl: clip, Range, Event, EPS, DEBUG
import GraphIdx: Ones, Vec, is_const

"""
Contains all memory needed for `line_las!`.
"""
struct LineDPMem{F}
    events::Vector{Event}
    lb::Vector{F}
end

LineDPMem{F}() where {F} =
    LineDPMem{F}([], [])

LineDPMem{F}(n::Integer) where {F} =
    resize!(LineDPMem{F}(), n)

function Base.resize!(mem::LineDPMem{F}, n::Integer)::LineDPMem{F} where {F}
    resize!(mem.events, 2n)
    resize!(mem.lb, n-1)
    mem
end

function line_las!(
    mem::LineDPMem{F},
    x::Array{F,N},
    y::Array{F,N},
    λ::Lam,
    µ::Mu,
)::Array{F,N} where {F, N, Lam, Mu}
    line_las!(mem, x, y, λ, μ, Val(!is_const(Mu)))
end

function line_las!(
    mem::LineDPMem{F},
    x::Array{F,N},
    y::Array{F,N},
    λ::Lam,
    µ::Mu,
    C::Val{check},
)::Array{F,N} where {F, N, Lam, Mu, check}
    n = length(y)
    @assert n == length(x)
    resize!(mem, n)
    local events::Vector{Event} = mem.events
    local lb::Vector{F} = mem.lb
    local ub::Vector{F} = x
    local pq = Ref(Range(n+1, n))
    local λ0::F = F(0.0)

    @static DEBUG && @info y
    @inbounds for i = 1:(n-1)
        lb[i] = clip(events, pq, +μ[i], -μ[i] * y[i] - λ0 + λ[i], Val(true), C)
        ub[i] = clip(events, pq, -μ[i], +μ[i] * y[i] - λ0 + λ[i], Val(false), C)
        λ0 = (!check || μ[i] > EPS) ? λ[i] : min(λ0, λ[i])
    end

    x[n] = clip(events, pq, +μ[n], -μ[n] * y[n] - λ0 + 0, Val(true), C)
    @static DEBUG && @printf("x[%d] = %g\n", n, x[n])
    @inbounds for i = n-1:-1:1
        x[i] = clamp(x[i+1], lb[i], ub[i])
    end

    return x
end

line_las!(x::Array{F,N}, y::Array{F,N}, λ::Lam, µ::Mu) where {F, N, Lam, Mu} =
    line_las!(LineDPMem{F}(N), x, y, λ, µ)

line_las(y::Array{F,N}, λ::Lam, µ::Mu = Ones{F}()) where {F, N, Lam, Mu} =
    line_las!(similar(y), y, λ, µ)

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
"
Author: Daniel O'Connor
Email: dt.oconnor@outlook.com
"

module SA
    include("tools.jl")
    using ..SamplingTools: hj_info, default_temp_range, generate_nbr_dict,
                            generate_nbr, discrete_aggregation, add_noise,
                            add_cross_talk, hj_info!, generate_nbr!
    using Random
    using Distributed
    using Profile
    using BenchmarkTools
    using LinearAlgebra
    using SharedArrays
    function anneal!(spin_vector::Array{Float64, 1},
                    nbs_idxs::Array{Int64, 2},
                    nbs_vals::Array{Float64, 2},
                    temp_schedule::Array{Float64, 1})
        maxnb = size(nbs_idxs, 1)
        nsteps = size(temp_schedule, 1)
        nspins = size(spin_vector, 1)
        perm = collect(1:nspins)
        istep = 1
        temp = 0.0
        ispin = 1
        spin_nb = 1
        deg = 1
        ediff =  0.0
        spin = 1
        idx = 1

        @fastmath @inbounds begin
            for istep = 1:nsteps
                temp = temp_schedule[istep]
                randperm!(perm)
                for idx = 1:nspins
                    ispin = perm[idx]
                    ediff = 0.0
                    spin = spin_vector[ispin]
                    for deg = 1:maxnb
                        spin_nb = nbs_idxs[deg, ispin]
                        if spin_nb === ispin
                            ediff += -2.0*spin*nbs_vals[deg, ispin]
                        elseif spin_nb !== 0
                            ediff += -2.0*spin*(nbs_vals[deg, ispin]*spin_vector[spin_nb])
                        end
                    end
                    if ediff <= 0.0
                        spin_vector[ispin] *= -1.0
                    elseif exp(-1.0 * ediff / temp) > rand()
                        spin_vector[ispin] *= -1.0
                    end
                end
            end
        end
        # return spin_vector
    end

    function sample(h::Dict, J::Dict; num_reads=1000, num_steps=1000,
                    add_xtalk=false, noisy=false,
                    temp_schedule=zeros(Float64, 10),
                    initial_state=Array{Float64, 1}(), return_raw=false)
        # forward anneal usage only
        nspins = length(h)
        JM, M, degs = hj_info(h, J)

        # prepare spin vectors
        if length(initial_state) == nspins
            @assert all(abs.(initial_state) .== 1.0)
            spin_vectors = [deepcopy(initial_state) for i =1:num_reads]
        else
            spin_vectors = rand([-1.0, 1.0], nspins, num_reads)
            spin_vectors = [spin_vectors[:, i] for i =1:num_reads]
        end

        # allow for custom temperature schedules, but if none given then
        # one will be calculated according to the num_steps
        if all(temp_schedule .== 0)
            hot_t, cold_t = default_temp_range(JM)
            temp_sch = collect(LinRange(hot_t, cold_t, num_steps))
        end

        # generate index and value array of problem
        nbs_idxs, nbs_vals = generate_nbr(nspins, M, maximum(degs) + 1)
        Js = copy(JM)

        if !noisy & !add_xtalk
            # standard anneal
            @inbounds for i = 1:num_reads
                anneal!(spin_vectors[i], nbs_idxs, nbs_vals, temp_schedule)
            end
        elseif noisy & add_xtalk
            # anneal with noise and x-talk, need to re-egenrate the problem
            # upon each read due to the change of parameter
            @inbounds for i = 1:num_reads
                hc, jc = add_cross_talk(add_noise(h, J)...)
                hj_info!(Js, M, degs, hc, jc)
                generate_nbr!(nbs_idxs, nbs_vals, nspins, M, maximum(degs) + 1)
                anneal!(spin_vectors[i], nbs_idxs, nbs_vals, temp_schedule)
            end
        elseif add_xtalk & !noisy
            # anneal with x-talk, need to re-egenrate the problem
            # once before the runs as problem is the same consistently
            hc, jc = add_cross_talk(h, J)
            hj_info!(Js, M, degs, hc, jc)
            generate_nbr!(nbs_idxs, nbs_vals, nspins, M, maximum(degs) + 1)
            @inbounds for i = 1:num_reads
                anneal!(spin_vectors[i], nbs_idxs, nbs_vals, temp_schedule)
            end
        else
            # anneal with noise, need to re-egenrate the problem
            # upon each read due to the change of parameters
            @inbounds for i = 1:num_reads
                hc, jc = add_noise(h, J)
                hj_info!(Js, M, degs, hc, jc)
                generate_nbr!(nbs_idxs, nbs_vals, nspins, M, maximum(degs) + 1)
                anneal!(spin_vectors[i], nbs_idxs, nbs_vals, temp_schedule)
            end
        end
        if return_raw
            return spin_vectors
        end
        return discrete_aggregation(spin_vectors, num_reads, collect(UpperTriangular(JM)))
    end

    export sample, anneal!
end

max_stars_repo_name: dtoconnor/Annealing-Samplers

max_issues_repo_name: dtoconnor/Annealing-Samplers

max_forks_repo_name: dtoconnor/Annealing-Samplers


text:
function interior_dGx_pos_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    thread_idx = threadIdx().x
    block_dim = blockDim().x

    sum_delx_sqr = 0.0
    sum_dely_sqr = 0.0
    sum_delx_dely = 0.0
    sum_delx_delf = SVector{4,Float64}(0, 0, 0, 0)
    sum_dely_delf = SVector{4,Float64}(0, 0, 0, 0)
    qtilde_i = SVector{4,Float64}(0, 0, 0, 0)
    qtilde_k = SVector{4,Float64}(0, 0, 0, 0)
    # phi_i = SVector{4,Float64}(0, 0, 0, 0)
    # phi_k = SVector{4,Float64}(0, 0, 0, 0)

    x_i = gpuGlobalDataFixedPoint[idx].x
    y_i = gpuGlobalDataFixedPoint[idx].y
    nx = gpuGlobalDataFixedPoint[idx].nx
    ny = gpuGlobalDataFixedPoint[idx].ny

    # limiter_flag = gpuConfigData[7]
    # power = gpuConfigData[6]
    # gamma = gpuConfigData[15]
    for iter in 15:24
        conn = gpuGlobalDataConn[iter, idx]
        if conn == 0
            break
        end

        delx = gpuGlobalDataFixedPoint[conn].x - x_i
        dely = gpuGlobalDataFixedPoint[conn].y - y_i
        dels = delx*ny - dely*nx
        deln = delx*nx + dely*ny
        # dist = CUDAnative.hypot(dels, deln)
        # weights = CUDAnative.pow(dist, gpuConfigData[6])
        weights = 1.0

        dels_weights = dels*weights
        deln_weights = deln*weights

        sum_delx_sqr += dels*dels_weights
        sum_dely_sqr += deln*deln_weights
        sum_delx_dely += dels*deln_weights

        # if limiter_flag == 1
            venkat_limiter_kernel(gpuGlobalDataFixedPoint, gpuGlobalDataRest, idx, gpuConfigData, delx, dely, shared, thread_idx, block_dim)
            qtilde_i = @SVector [ gpuGlobalDataRest[9, idx] - 0.5*shared[thread_idx]*(delx * gpuGlobalDataRest[13, idx] + dely * gpuGlobalDataRest[17, idx]),
                gpuGlobalDataRest[10, idx] - 0.5*shared[thread_idx + block_dim]*(delx * gpuGlobalDataRest[14, idx] + dely * gpuGlobalDataRest[18, idx]),
                gpuGlobalDataRest[11, idx] - 0.5*shared[thread_idx + block_dim * 2]*(delx * gpuGlobalDataRest[15, idx] + dely * gpuGlobalDataRest[19, idx]),
                gpuGlobalDataRest[12, idx] - 0.5*shared[thread_idx + block_dim * 3]*(delx * gpuGlobalDataRest[16, idx] + dely * gpuGlobalDataRest[20, idx])]
            venkat_limiter_kernel(gpuGlobalDataFixedPoint, gpuGlobalDataRest, conn, gpuConfigData, delx, dely, shared, thread_idx, block_dim)
            qtilde_k = @SVector [ gpuGlobalDataRest[9, conn] - 0.5*shared[thread_idx]*(delx * gpuGlobalDataRest[13, conn] + dely * gpuGlobalDataRest[17, conn]),
                gpuGlobalDataRest[10, conn] - 0.5*shared[thread_idx + block_dim]*(delx * gpuGlobalDataRest[14, conn] + dely * gpuGlobalDataRest[18, conn]),
                gpuGlobalDataRest[11, conn] - 0.5*shared[thread_idx + block_dim * 2]*(delx * gpuGlobalDataRest[15, conn] + dely * gpuGlobalDataRest[19, conn]),
                gpuGlobalDataRest[12, conn] - 0.5*shared[thread_idx + block_dim * 3]*(delx * gpuGlobalDataRest[16, conn] + dely * gpuGlobalDataRest[20, conn])]
        # end

        shared[thread_idx], shared[thread_idx + block_dim], shared[thread_idx + block_dim * 2], shared[thread_idx + block_dim * 3] = 0.0, 0.0, 0.0, 0.0

        qtilde_to_primitive_kernel(qtilde_i, gpuConfigData, shared, thread_idx, block_dim)
        flux_Gxp_kernel(nx, ny, idx, shared, +, thread_idx, block_dim)
        qtilde_to_primitive_kernel(qtilde_k, gpuConfigData, shared, thread_idx, block_dim)
        flux_Gxp_kernel(nx, ny, idx, shared, -, thread_idx, block_dim)

        # CUDAnative.synchronize()
        temp_var = @SVector [shared[thread_idx], shared[thread_idx + block_dim], shared[thread_idx + block_dim * 2], shared[thread_idx + block_dim * 3]]
        sum_delx_delf += temp_var * dels_weights
        sum_dely_delf += temp_var * deln_weights
    end

    det = sum_delx_sqr*sum_dely_sqr - sum_delx_dely*sum_delx_dely
    one_by_det = 1.0 / det

    flux_shared[thread_idx] += (sum_delx_delf[1]*sum_dely_sqr - sum_dely_delf[1]*sum_delx_dely)*one_by_det
    flux_shared[thread_idx + block_dim] += (sum_delx_delf[2]*sum_dely_sqr - sum_dely_delf[2]*sum_delx_dely)*one_by_det
    flux_shared[thread_idx + block_dim * 2] += (sum_delx_delf[3]*sum_dely_sqr - sum_dely_delf[3]*sum_delx_dely)*one_by_det
    flux_shared[thread_idx + block_dim * 3] += (sum_delx_delf[4]*sum_dely_sqr - sum_dely_delf[4]*sum_delx_dely)*one_by_det
    # if idx == 1
    #     @cuprintf("\n %f %f %f %f", gpuGlobalDataRest[5, idx],gpuGlobalDataRest[6, idx],gpuGlobalDataRest[7, idx],gpuGlobalDataRest[8, idx])
    # end
    return nothing
end

function interior_dGx_neg_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    thread_idx = threadIdx().x
    block_dim = blockDim().x

    sum_delx_sqr = 0.0
    sum_dely_sqr = 0.0
    sum_delx_dely = 0.0
    sum_delx_delf = SVector{4,Float64}(0, 0, 0, 0)
    sum_dely_delf = SVector{4,Float64}(0, 0, 0, 0)
    qtilde_i = SVector{4,Float64}(0, 0, 0, 0)
    qtilde_k = SVector{4,Float64}(0, 0, 0, 0)

    x_i = gpuGlobalDataFixedPoint[idx].x
    y_i = gpuGlobalDataFixedPoint[idx].y
    nx = gpuGlobalDataFixedPoint[idx].nx
    ny = gpuGlobalDataFixedPoint[idx].ny

    # power = gpuConfigData[6]
    # gamma = gpuConfigData[15]
    for iter in 25:34
        conn = gpuGlobalDataConn[iter, idx]
        if conn == 0
            break
        end

        delx = gpuGlobalDataFixedPoint[conn].x - x_i
        dely = gpuGlobalDataFixedPoint[conn].y - y_i
        dels = delx*ny - dely*nx
        deln = delx*nx + dely*ny
        # dist = CUDAnative.hypot(dels, deln)
        # weights = CUDAnative.pow(dist, gpuConfigData[6])
        weights = 1.0

        dels_weights = dels*weights
        deln_weights = deln*weights

        sum_delx_sqr += dels*dels_weights
        sum_dely_sqr += deln*deln_weights
        sum_delx_dely += dels*deln_weights

        # if limiter_flag == 1
            venkat_limiter_kernel(gpuGlobalDataFixedPoint, gpuGlobalDataRest, idx, gpuConfigData, delx, dely, shared, thread_idx, block_dim)
            qtilde_i = @SVector [ gpuGlobalDataRest[9, idx] - 0.5*shared[thread_idx]*(delx * gpuGlobalDataRest[13, idx] + dely * gpuGlobalDataRest[17, idx]),
                gpuGlobalDataRest[10, idx] - 0.5*shared[thread_idx + block_dim]*(delx * gpuGlobalDataRest[14, idx] + dely * gpuGlobalDataRest[18, idx]),
                gpuGlobalDataRest[11, idx] - 0.5*shared[thread_idx + block_dim * 2]*(delx * gpuGlobalDataRest[15, idx] + dely * gpuGlobalDataRest[19, idx]),
                gpuGlobalDataRest[12, idx] - 0.5*shared[thread_idx + block_dim * 3]*(delx * gpuGlobalDataRest[16, idx] + dely * gpuGlobalDataRest[20, idx])]
            venkat_limiter_kernel(gpuGlobalDataFixedPoint, gpuGlobalDataRest, conn, gpuConfigData, delx, dely, shared, thread_idx, block_dim)
            qtilde_k = @SVector [ gpuGlobalDataRest[9, conn] - 0.5*shared[thread_idx]*(delx * gpuGlobalDataRest[13, conn] + dely * gpuGlobalDataRest[17, conn]),
                gpuGlobalDataRest[10, conn] - 0.5*shared[thread_idx + block_dim]*(delx * gpuGlobalDataRest[14, conn] + dely * gpuGlobalDataRest[18, conn]),
                gpuGlobalDataRest[11, conn] - 0.5*shared[thread_idx + block_dim * 2]*(delx * gpuGlobalDataRest[15, conn] + dely * gpuGlobalDataRest[19, conn]),
                gpuGlobalDataRest[12, conn] - 0.5*shared[thread_idx + block_dim * 3]*(delx * gpuGlobalDataRest[16, conn] + dely * gpuGlobalDataRest[20, conn])]
        # end

        shared[thread_idx], shared[thread_idx + block_dim], shared[thread_idx + block_dim * 2], shared[thread_idx + block_dim * 3] = 0.0, 0.0, 0.0, 0.0

        qtilde_to_primitive_kernel(qtilde_i, gpuConfigData, shared, thread_idx, block_dim)
        flux_Gxn_kernel(nx, ny, idx, shared, +, thread_idx, block_dim)
        qtilde_to_primitive_kernel(qtilde_k, gpuConfigData, shared, thread_idx, block_dim)
        flux_Gxn_kernel(nx, ny, idx, shared, -, thread_idx, block_dim)
        # CUDAnative.synchronize()
        temp_var = @SVector [shared[thread_idx], shared[thread_idx + block_dim], shared[thread_idx + block_dim * 2], shared[thread_idx + block_dim * 3] ]
        sum_delx_delf += temp_var * dels_weights
        sum_dely_delf += temp_var * deln_weights
    end

    det = sum_delx_sqr*sum_dely_sqr - sum_delx_dely*sum_delx_dely
    one_by_det = 1.0 / det
    flux_shared[thread_idx] += (sum_delx_delf[1]*sum_dely_sqr - sum_dely_delf[1]*sum_delx_dely)*one_by_det
    flux_shared[thread_idx + block_dim] += (sum_delx_delf[2]*sum_dely_sqr - sum_dely_delf[2]*sum_delx_dely)*one_by_det
    flux_shared[thread_idx + block_dim * 2] += (sum_delx_delf[3]*sum_dely_sqr - sum_dely_delf[3]*sum_delx_dely)*one_by_det
    flux_shared[thread_idx + block_dim * 3] += (sum_delx_delf[4]*sum_dely_sqr - sum_dely_delf[4]*sum_delx_dely)*one_by_det

    return nothing
end

function interior_dGy_pos_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    thread_idx = threadIdx().x
    block_dim = blockDim().x

    sum_delx_sqr = 0.0
    sum_dely_sqr = 0.0
    sum_delx_dely = 0.0
    sum_delx_delf = SVector{4,Float64}(0, 0, 0, 0)
    sum_dely_delf = SVector{4,Float64}(0, 0, 0, 0)
    qtilde_i = SVector{4,Float64}(0, 0, 0, 0)
    qtilde_k = SVector{4,Float64}(0, 0, 0, 0)

    x_i = gpuGlobalDataFixedPoint[idx].x
    y_i = gpuGlobalDataFixedPoint[idx].y
    nx = gpuGlobalDataFixedPoint[idx].nx
    ny = gpuGlobalDataFixedPoint[idx].ny

    # power = gpuConfigData[6]
    # gamma = gpuConfigData[15]

    for iter in 35:44
        conn = gpuGlobalDataConn[iter, idx]
        if conn == 0
            break
        end

        delx = gpuGlobalDataFixedPoint[conn].x - x_i
        dely = gpuGlobalDataFixedPoint[conn].y - y_i
        dels = delx*ny - dely*nx
        deln = delx*nx + dely*ny
        # dist = CUDAnative.hypot(dels, deln)
        # weights = CUDAnative.pow(dist, gpuConfigData[6])
        weights = 1.0

        dels_weights = dels*weights
        deln_weights = deln*weights

        sum_delx_sqr += dels*dels_weights
        sum_dely_sqr += deln*deln_weights
        sum_delx_dely += dels*deln_weights

        # if limiter_flag == 1
            venkat_limiter_kernel(gpuGlobalDataFixedPoint, gpuGlobalDataRest, idx, gpuConfigData, delx, dely, shared, thread_idx, block_dim)
            qtilde_i = @SVector [ gpuGlobalDataRest[9, idx] - 0.5*shared[thread_idx]*(delx * gpuGlobalDataRest[13, idx] + dely * gpuGlobalDataRest[17, idx]),
                gpuGlobalDataRest[10, idx] - 0.5*shared[thread_idx + block_dim]*(delx * gpuGlobalDataRest[14, idx] + dely * gpuGlobalDataRest[18, idx]),
                gpuGlobalDataRest[11, idx] - 0.5*shared[thread_idx + block_dim * 2]*(delx * gpuGlobalDataRest[15, idx] + dely * gpuGlobalDataRest[19, idx]),
                gpuGlobalDataRest[12, idx] - 0.5*shared[thread_idx + block_dim * 3]*(delx * gpuGlobalDataRest[16, idx] + dely * gpuGlobalDataRest[20, idx])]
            venkat_limiter_kernel(gpuGlobalDataFixedPoint, gpuGlobalDataRest, conn, gpuConfigData, delx, dely, shared, thread_idx, block_dim)
            qtilde_k = @SVector [ gpuGlobalDataRest[9, conn] - 0.5*shared[thread_idx]*(delx * gpuGlobalDataRest[13, conn] + dely * gpuGlobalDataRest[17, conn]),
                gpuGlobalDataRest[10, conn] - 0.5*shared[thread_idx + block_dim]*(delx * gpuGlobalDataRest[14, conn] + dely * gpuGlobalDataRest[18, conn]),
                gpuGlobalDataRest[11, conn] - 0.5*shared[thread_idx + block_dim * 2]*(delx * gpuGlobalDataRest[15, conn] + dely * gpuGlobalDataRest[19, conn]),
                gpuGlobalDataRest[12, conn] - 0.5*shared[thread_idx + block_dim * 3]*(delx * gpuGlobalDataRest[16, conn] + dely * gpuGlobalDataRest[20, conn])]
        # end

        # if limiter_flag == 2
        #     @cuprintf("\n Havent written the code - die \n")
        # end
        shared[thread_idx], shared[thread_idx + block_dim], shared[thread_idx + block_dim * 2], shared[thread_idx + block_dim * 3] = 0.0, 0.0, 0.0, 0.0

        qtilde_to_primitive_kernel(qtilde_i, gpuConfigData, shared, thread_idx, block_dim)
        flux_Gyp_kernel(nx, ny, idx, shared, +, thread_idx, block_dim)
        qtilde_to_primitive_kernel(qtilde_k, gpuConfigData, shared, thread_idx, block_dim)
        flux_Gyp_kernel(nx, ny, idx, shared, -, thread_idx, block_dim)
        # CUDAnative.synchronize()
        temp_var = @SVector [shared[thread_idx], shared[thread_idx + block_dim], shared[thread_idx + block_dim * 2], shared[thread_idx + block_dim * 3] ]
        sum_delx_delf += temp_var * dels_weights
        sum_dely_delf += temp_var * deln_weights

    end

    det = sum_delx_sqr*sum_dely_sqr - sum_delx_dely*sum_delx_dely
    one_by_det = 1.0 / det
    flux_shared[thread_idx] += (sum_dely_delf[1]*sum_delx_sqr - sum_delx_delf[1]*sum_delx_dely)*one_by_det
    flux_shared[thread_idx + block_dim] += (sum_dely_delf[2]*sum_delx_sqr - sum_delx_delf[2]*sum_delx_dely)*one_by_det
    flux_shared[thread_idx + block_dim * 2] += (sum_dely_delf[3]*sum_delx_sqr - sum_delx_delf[3]*sum_delx_dely)*one_by_det
    flux_shared[thread_idx + block_dim * 3] += (sum_dely_delf[4]*sum_delx_sqr - sum_delx_delf[4]*sum_delx_dely)*one_by_det
    # if idx ==1
    #     @cuprintf("\n %f %f %f %f", gpuGlobalDataRest[5, idx],gpuGlobalDataRest[6, idx],gpuGlobalDataRest[7, idx],gpuGlobalDataRest[8, idx])
    # end
    return nothing
end

function interior_dGy_neg_kernel(gpuGlobalDataConn, gpuGlobalDataFixedPoint, gpuGlobalDataRest, gpuConfigData, shared, flux_shared)
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    thread_idx = threadIdx().x
    block_dim = blockDim().x

    sum_delx_sqr = 0.0
    sum_dely_sqr = 0.0
    sum_delx_dely = 0.0
    sum_delx_delf = SVector{4,Float64}(0, 0, 0, 0)
    sum_dely_delf = SVector{4,Float64}(0, 0, 0, 0)
    qtilde_i = SVector{4,Float64}(0, 0, 0, 0)
    qtilde_k = SVector{4,Float64}(0, 0, 0, 0)

    x_i = gpuGlobalDataFixedPoint[idx].x
    y_i = gpuGlobalDataFixedPoint[idx].y
    nx = gpuGlobalDataFixedPoint[idx].nx
    ny = gpuGlobalDataFixedPoint[idx].ny

    # power = gpuConfigData[6]
    # gamma = gpuConfigData[15]

    for iter in 45:54
        conn = gpuGlobalDataConn[iter, idx]
        if conn == 0
            break
        end

        delx = gpuGlobalDataFixedPoint[conn].x - x_i
        dely = gpuGlobalDataFixedPoint[conn].y - y_i
        dels = delx*ny - dely*nx
        deln = delx*nx + dely*ny
        # dist = CUDAnative.hypot(dels, deln)
        # weights = CUDAnative.pow(dist, gpuConfigData[6])
        weights = 1.0

        dels_weights = dels*weights
        deln_weights = deln*weights

        sum_delx_sqr += dels*dels_weights
        sum_dely_sqr += deln*deln_weights
        sum_delx_dely += dels*deln_weights

        # if limiter_flag == 1
            venkat_limiter_kernel(gpuGlobalDataFixedPoint, gpuGlobalDataRest, idx, gpuConfigData, delx, dely, shared, thread_idx, block_dim)
            qtilde_i = @SVector [ gpuGlobalDataRest[9, idx] - 0.5*shared[thread_idx]*(delx * gpuGlobalDataRest[13, idx] + dely * gpuGlobalDataRest[17, idx]),
                gpuGlobalDataRest[10, idx] - 0.5*shared[thread_idx + block_dim]*(delx * gpuGlobalDataRest[14, idx] + dely * gpuGlobalDataRest[18, idx]),
                gpuGlobalDataRest[11, idx] - 0.5*shared[thread_idx + block_dim * 2]*(delx * gpuGlobalDataRest[15, idx] + dely * gpuGlobalDataRest[19, idx]),
                gpuGlobalDataRest[12, idx] - 0.5*shared[thread_idx + block_dim * 3]*(delx * gpuGlobalDataRest[16, idx] + dely * gpuGlobalDataRest[20, idx])]
            venkat_limiter_kernel(gpuGlobalDataFixedPoint, gpuGlobalDataRest, conn, gpuConfigData, delx, dely, shared, thread_idx, block_dim)
            qtilde_k = @SVector [ gpuGlobalDataRest[9, conn] - 0.5*shared[thread_idx]*(delx * gpuGlobalDataRest[13, conn] + dely * gpuGlobalDataRest[17, conn]),
                gpuGlobalDataRest[10, conn] - 0.5*shared[thread_idx + block_dim]*(delx * gpuGlobalDataRest[14, conn] + dely * gpuGlobalDataRest[18, conn]),
                gpuGlobalDataRest[11, conn] - 0.5*shared[thread_idx + block_dim * 2]*(delx * gpuGlobalDataRest[15, conn] + dely * gpuGlobalDataRest[19, conn]),
                gpuGlobalDataRest[12, conn] - 0.5*shared[thread_idx + block_dim * 3]*(delx * gpuGlobalDataRest[16, conn] + dely * gpuGlobalDataRest[20, conn])]
        # end
        # if limiter_flag == 2
        #     @cuprintf("\n Havent written the code - die \n")
        # end
        shared[thread_idx], shared[thread_idx + block_dim], shared[thread_idx + block_dim * 2], shared[thread_idx + block_dim * 3] = 0.0, 0.0, 0.0, 0.0

        qtilde_to_primitive_kernel(qtilde_i, gpuConfigData, shared, thread_idx, block_dim)
        flux_Gyn_kernel(nx, ny, idx, shared, +, thread_idx, block_dim)
        qtilde_to_primitive_kernel(qtilde_k, gpuConfigData, shared, thread_idx, block_dim)
        flux_Gyn_kernel(nx, ny, idx, shared, -, thread_idx, block_dim)
        # CUDAnative.synchronize()
        temp_var = @SVector [shared[thread_idx], shared[thread_idx + block_dim], shared[thread_idx + block_dim * 2], shared[thread_idx + block_dim * 3] ]
        sum_delx_delf += temp_var * dels_weights
        sum_dely_delf += temp_var * deln_weights

    end

    det = sum_delx_sqr*sum_dely_sqr - sum_delx_dely*sum_delx_dely
    one_by_det = 1.0 / det
    flux_shared[thread_idx] += (sum_dely_delf[1]*sum_delx_sqr - sum_delx_delf[1]*sum_delx_dely)*one_by_det
    flux_shared[thread_idx + block_dim] += (sum_dely_delf[2]*sum_delx_sqr - sum_delx_delf[2]*sum_delx_dely)*one_by_det
    flux_shared[thread_idx + block_dim * 2] += (sum_dely_delf[3]*sum_delx_sqr - sum_delx_delf[3]*sum_delx_dely)*one_by_det
    flux_shared[thread_idx + block_dim * 3] += (sum_dely_delf[4]*sum_delx_sqr - sum_delx_delf[4]*sum_delx_dely)*one_by_det
    # if idx ==1
    #     @cuprintf("\n %f %f %f %f", gpuGlobalDataRest[5, idx],gpuGlobalDataRest[6, idx],gpuGlobalDataRest[7, idx],gpuGlobalDataRest[8, idx])
    # end
    return nothing
end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
#! /usr/bin/env julia
function incDict!(dict, key)
  value = get!(dict, key, 0) + 1
  setindex!(dict, value, key)
  return dict
end

charCounts = str -> foldl(incDict!, str; init=Dict())

counts = map(charCounts, readlines())
twos = sum(count -> 2 in values(count) ? 1 : 0, counts)
threes = sum(count -> 3 in values(count) ? 1 : 0, counts)
println(twos * threes)

max_stars_repo_name: paolobueno/advent-of-julia

max_issues_repo_name: paolobueno/advent-of-julia

max_forks_repo_name: paolobueno/advent-of-julia


text:
"""
TODO: Description of _deme_ module.
"""
module _deme_
importall common
using core, breeder, species, fitness
export Deme, contents, deme, DemeDefinition

type Deme{F}
  capacity::Int
  breeder::Breeder
  species::Species
  num_offspring::Int
  offspring::IndividualCollection{F}
  members::IndividualCollection{F}

  Deme(capacity::Int, breeder::Breeder, species::Species, num_offspring::Int) =
    new(capacity, breeder, species, num_offspring, empty_individual_collection(species), empty_individual_collection(species))
end

"""
Used to provide a definition for building a deme.
"""
type DemeDefinition
  """
  The number of individuals within this deme.
  """
  capacity::Int

  """
  The number of offspring produced by this deme at each generation.
  If set to `0`, the number of offspring will default to being equal to
  the capacity of the deme.
  """
  offspring::Int

  """
  The species to which all members of this deme belong.
  """
  species::SpeciesDefinition

  """
  The breeder used to produce the offspring for this deme.
  """
  breeder::BreederDefinition

  """
  Constructs a partial deme specification using default values.
  """
  DemeDefinition() = new(100, -1)
  DemeDefinition(
    capacity::Int,
    offspring::Int,
    species::SpeciesDefinition,
    breeder::BreederDefinition
  ) =
    new(capacity, offspring, species, breeder)
end

"""
Composes a deme from a provided specification.
"""
function compose!(d::DemeDefinition)
  sp = species.compose!(d.species)
  br = breeder.compose!(d.breeder, sp)
  if d.offspring == 0
    d.offspring = d.capacity
  end
  Deme{uses(sp.fitness)}(d.capacity, br, sp, d.offspring)
end

"""
TODO: Description of deme models.

**Properties:**

* `capacity::Int`, the number of individuals that this deme may hold.
* `offspring::Int`, the number of offspring that should be produced by this
deme at every generation.
* `species::SpeciesSpecification`, the species to which the individuals in the
deme belong.
* `breeder::BreederSpecification`, the breeder used to produce the offspring
for this deme.
"""
deme( capacity::Int,
      offspring::Int,
      species::species.SpeciesDefinition,
      breeder::breeder.BreederDefinition) =
  DemeDefinition(capacity, offspring, species, breeder)

function deme(spec::Function)
  ds = DemeDefinition()
  spec(ds)
  ds
end

"""
Produces the offspring for a given deme at each generation.
"""
breed!(d::Deme) =
  breed!(d.breeder, d.species, d.offspring, d.members, d.num_offspring)
end


text:
__precompile__(true)
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module Triangles

using Nemo, SeqUtils

export ModuleTriangles
export Triangle, ZTriangle, QTriangle, RecTriangle
export TriangularNumber, isTriangular, assertTriangular
export ShowAsΔ, ShowAsMatrix, Row, RowSums
export fromΔ, toΔ, TriangleToList

"""
* Triangle, ZTriangle, QTriangle, RecTriangle, TriangularNumber, isTriangular, assertTriangular, ShowAsΔ, ShowAsMatrix, Row, RowSums, fromΔ, toΔ, TriangleToList
"""
const ModuleTriangles = ""

#    T(0,0)                          row 0
#    T(1,0) T(1,1)                   row 1
#    T(2,0) T(2,1) T(2,2)            row 2
#    T(3,0) T(3,1) T(3,2) T(3,3)     row 3
#    col 0  col 1  col 2  col 3
#
# The two dimensional index is ``(0, 0)``-based and displayed row-wise as
# a lower triangular array. In particular note that a `Triangle` has rows
# and cols from ``0`` to dim ``-1``.

abstract type AbstractTriangle end

"""
A Triangle is a list of rows. The rows of the triangle are generated by a function `(n::Int, k::Int) → gen::fmpz` for ``0 ≤ k ≤ n`` and ``0 ≤ n <``dim.
"""
struct Triangle <: AbstractTriangle
    dim::Int
    gen::Function
    function Triangle(dim::Int, gen::Function)
        dim <= 0 && throw(ArgumentError("dim has to be a positive integer."))
        new(dim, (n::Int, k::Int) -> gen(n, k)::fmpz )
    end
end

"""
Return the first row of the triangle.
"""
function Base.iterate(T::Triangle)
    ([T.gen(0, 0)], 1)
end

"""
Return the next row of the triangle.
"""
function Base.iterate(T::Triangle, n)
    n >= T.dim && return nothing
    row = [T.gen(n, k) for k in 0:n]
    (row, n + 1)
end

Base.length(T::Triangle) = T.dim
Base.eltype(T::Triangle) = fmpz

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

"""
A recursive triangle `RecTriangle` is a subtype of `AbstractTriangle`.
The rows of the triangle are generated by a function `gen(n, k, prevrow)` defined for ``n ≥ 0`` and ``0 ≤ k ≤ n``. The function returns value of type Nemo.fmpz.

The parameter prevrow is a function which returns the values of `row(n-1)` of the triangle and 0 if ``k < 0`` or ``k > n``. The function prevrow is provided
by an instance of `RecTriangle` and must not be defined by the user.
"""
struct RecTriangle <: AbstractTriangle
    dim::Int
    A::Array{fmpz}
    gen::Function  # generating function
    function RecTriangle(dim::Int, gen::Function)
        dim <= 0 && throw(ArgumentError("dim has to be a positive integer."))
        new( dim,
             fill(ZZ(0), dim),
             (n::Int, k::Int, f::Function) -> gen(n, k, f)::fmpz)
    end
end

function Base.iterate(T::RecTriangle)
    top = T.gen(0, 0, () -> undefined)
    T.A[1] = ZZ(top)
    ([top], 1)
end

"""
Return the next row of the triangle.
"""
function Base.iterate(T::RecTriangle, n)
    @inline prevrow(A, n) = (k) -> (k<0 || k>n) ? ZZ(0) : A[k+1]
    n == T.dim && return nothing
    F(n, k) = T.gen(n, k, prevrow(T.A, n))
    row = [F(n, k) for k in 0:n]
    for k in 1:n+1 T.A[k] = row[k] end
    (row, n + 1)
end

Base.length(R::RecTriangle) = R.dim
Base.eltype(R::RecTriangle) = fmpz

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

"""
Return the``n``-th triangular number.
"""
TriangularNumber(n) = div(n * (n + 1), 2)

"""
Is ``n`` a triangular number?
"""
isTriangular(n) = n == TriangularNumber(isqrt(2n))

"""
Return the square root of ``2n`` or throw an ArgumentError if ``n`` is not a triangular number.
"""
function assertTriangular(n)
    dim = isqrt(2n)
    n ≠ TriangularNumber(dim) && throw(ArgumentError("This is not a triangular array!"))
    return dim
end

"""
Return a trianguler array with ``n`` rows set to ``0`` (type ZZ).
"""
ZTriangle(dim::Int) = zeros(ZZ, TriangularNumber(dim))

"""
Return a trianguler array with ``n`` rows set to ``0`` (type QQ).
"""
QTriangle(dim::Int) = zeros(QQ, TriangularNumber(dim))

"""
Return the triangle as a list of rows.
"""
TriangularArray(T::AbstractTriangle) = [row for row in T]

"""
Return the triangle as a list of integers.
"""
TriangleToList(T::AbstractTriangle) = [k for row in T for k in row]

"""
Return row ``n (0 ≤ n)`` of the lower triangular matrix T or the row ``n`` in reversed order if `reversed` is set true.
"""
function Row(T, n, reversed=false)
    s = length(T)
    (s == 0 || n < 0) && return []
    assertTriangular(s)

    t = TriangularNumber(n+1)
    s < t && error("This row is not in the matrix.")
    R = reversed ? range(t, step=-1, stop=t-n) : range(t-n, step=1, stop=t)
    return [T[k] for k in R]
end

"""
Return the row sums of a triangle, if `alternate=true` then the alternating row sums.
"""
function RowSums(T, alternate = false)
    n = length(T)
    n == 0 && return T
    len = assertTriangular(n)
    S = zeros(ZZ, len)
    lo = hi = 0
    step = 0

    while hi < n
        if alternate
            s = sum((-1)^k * T[lo+k+1] for k in 0:hi - lo)
        else
            s = sum(T[k+1] for k in lo:hi)
        end
        step += 1
        S[step] = s
        hi += 1
        lo, hi = hi, hi + step
    end
    S
end

"""
Display a lower triangular matrix.
"""
function ShowAsΔ(T::AbstractTriangle, separator=", ")
    for row in T
        for k in row print(k, separator) end
        println()
    end
end

"""
Display an array as a lower triangular matrix.
"""
function ShowAsΔ(T::AbstractArray, separator=", ")
    if typeof(T) == Array{Array{fmpz,1},1}
        for row in T
            for k in row print(k, separator) end
            println()
        end
    elseif typeof(T) == Array{fmpz,2}
        for row in 1:size(T,1)
            T[row,:]' |> Println
        end
    elseif typeof(T) == Array{fmpz,1}
        dim = assertTriangular(length(T))
        low = high = 1
        for r in 2:dim+1
            [T[j] for j = low:high] |> Println
            low, high = high + 1, high + r
        end
    else
        println(T)
    end
end

"""
Convert a lower triangular array to a square matrix.
"""
function fromΔ(T::AbstractTriangle)
    L = TriangleToList(T)
    dim = assertTriangular(length(L))

    A = fill(ZZ(0), dim, dim)
    k = 1
    for r in 1:dim
        for j = 1:r
            A[r, j] = L[k]
            k += 1
        end
    end
    A
end

"""
Convert a square matrix to a list using only the ``T(n,k)`` with ``0 ≤ k ≤ n``. Returns a ZArray.
"""
function toΔ(M)
    dim = size(M, 1)
    T = ZTriangle(dim)
    k = 1

    if typeof(M[1, 1]) == Nemo.fmpq
        for r in 1:dim, j = 1:r
            T[k] = numerator(M[r, j])
            k += 1
        end
    else
        for r in 1:dim, j = 1:r
            T[k] = M[r, j]
            k += 1
        end
    end
    T
end

"""
Print the triangle in matrix form.
"""
ShowAsMatrix(T) = println(fromΔ(T))

#START-TEST-########################################################

using Test

function test()
    BinTri(n) = [binomial(k, j) for k in 0:n - 1 for j in 0:k]

    X(n, k) = ((k > n || k < 0) && return 0;
        (n == 0 && k == 0) && return 1;
        4 * X(n - 1, k - 1) + (4 * n - 1) * X(n - 1, k); )

    T(n) = [X(j, k) for j in 0:n for k in 0:j]

    @testset "Triangles" begin
        t = T(6)
        a = [3465, 8784, 7136, 2304, 256]
        b = Row(t, 4)
        @test all(a .== b)

        # RowSums
        t = BinTri(8)
        a = [1, 2, 4, 8, 16, 32, 64, 128]
        b = RowSums(t)
        @test all(a .== b)

        a = [1, 0, 0, 0, 0, 0, 0, 0]
        b = RowSums(t, true)  # alternating sum
        @test all(a .== b)
    end
end

function demo()

    println("\n--- Using Triangle\n")

    function X(n, k)
        (k > n || k < 0) && return fmpz(0)
        (n == 0 && k == 0) && return fmpz(1)
        ZZ( 4 * X(n - 1, k - 1) + (4 * n - 1) * X(n - 1, k) )
    end

    TX(dim) = Triangle(dim, X)
    ShowAsΔ(TX(6), "\t")

    println(TriangleToList(TX(6)))
    ShowAsΔ(TriangleToList(TX(6)))

    TL = TriangleToList(TX(6))
    for k in 0:5
        Println(Row(TL, k))
    end

    println()
    ShowAsMatrix(TX(6))

    # The above triangle is much better implemented with a RecTriangle!

    println("\n--- Using RecTriangle\n")

    function r225478(n, k, prevrow::Function)
        (k == 0 && n == 0) && return fmpz(1)
        4 * prevrow(k - 1) + (4 * n - 1) * prevrow(k)
    end

    I225478(dim) = RecTriangle(dim, r225478)
    for row in I225478(6) Println(row) end

    println()
    ShowAsΔ(RecTriangle(6, r225478), "\t")

    println()
    TA225478 = TriangularArray(I225478(6))
    ShowAsΔ(TA225478)

    println()
    println("ΔToListToΔ (test)")
    ShowAsΔ(TriangleToList(I225478(6)))

    println()
    ShowAsMatrix(I225478(6))

    println()
    GC.gc()
    @time TriangularArray(TX(16))
    @time TriangularArray(I225478(16))
end

"""
TriangularArray(TX(16))
    0.063180 seconds (785.95 k allocations: 11.996 MiB)
TriangularArray(I225478(16))
    0.000264 seconds (653 allocations: 15.781 KiB)
"""
function perf()
    #GC.gc()
    #@time TriangularArray(TX(20))
    #@time TriangularArray(I225478(20))
 end

 function main()
     test()
     demo()
     perf()
 end

 main()

 end # module


text:
export SingleLayer

"""
singleLayer

σ(K*s+b)

where K,b are trainable weights

"""
struct SingleLayer
end

mσ(x::AbstractArray{R}) where R<:Real = abs.(x)+log.(R(1) .+ exp.(-R(2)*abs.(x)))
mdσ(x::AbstractArray{R}) where R<:Real = tanh.(x)
md2σ(x::AbstractArray{R}) where R<:Real = one(eltype(x)) .- tanh.(x).^2

"""
evaluate layer for current weights Θ=(K,b)
"""
function (N::SingleLayer)(S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (K,b) = Θ
    return mσ(K*S .+ b)
end

"""
compute matvec J_S N(S,Θ)'*Z
"""
function getJSTmv(N::SingleLayer,Z::AbstractArray{R},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (K,b) = Θ
    return K'*(mdσ(K*S .+ b) .* Z)
end

function getJSTmv(N::SingleLayer,Z::AbstractArray{R,3},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    t1 = reshape(mdσ(K*S .+ b),size(K,1),1,nex)
    t2 = K'*reshape(t1 .* Z, size(K,1),:)
    return reshape(t2,size(K,2),size(Z,2),nex)
end

"""
compute hessian matvec
"""

function getTraceHessAndGrad(N::SingleLayer, w::AbstractArray{R},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    t1 = K*S .+ b
    Jac = reshape(mdσ(t1),size(K,1),1,nex) .* K
    return vec(sum(reshape(md2σ(t1) .* w,size(K,1),:,nex).*(K.^2),dims=(1,2))), Jac
end

function getTraceHessAndGrad(N::SingleLayer, w::AbstractArray{R},Jac::AbstractArray{R,2},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    t1 = K*S .+ b;
    Jac = K * Jac
    trH =  vec(sum(reshape(md2σ(t1) .* w,size(K,1),:,nex).*(Jac).^2,dims=(1,2)))
    Jac = reshape(mdσ(t1),size(K,1),1,nex) .* Jac
    return trH,Jac
end

function getTraceHessAndGrad(N::SingleLayer, w::AbstractArray{R},Jac::AbstractArray{R,3},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    t1 = K*S .+ b;
    Jac = reshape(K* reshape(Jac,size(K,2),:), size(K,1),:,nex)
    trH =  vec(sum(reshape(md2σ(t1) .* w,size(K,1),:,nex).*(Jac).^2,dims=(1,2)))
    Jac = reshape(mdσ(t1),size(K,1),1,nex) .* Jac
    return trH,Jac
end

function getTraceHess(N::SingleLayer, w::AbstractArray{R},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    t1 = K*S .+ b
    return vec(sum(reshape(md2σ(t1) .* w,size(K,1),:,nex).*(K.^2),dims=(1,2)))
end

function getTraceHess(N::SingleLayer, w::AbstractArray{R},Jac::AbstractArray{R,2},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    t1 = K*S .+ b;
    Jac = K * Jac
    trH =  vec(sum(reshape(md2σ(t1) .* w,size(K,1),:,nex).*(Jac).^2,dims=(1,2)))
    return trH
end

function getTraceHess(N::SingleLayer, w::AbstractArray{R},Jac::AbstractArray{R,3},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    t1 = K*S .+ b;
    Jac = reshape(K* reshape(Jac,size(K,2),:), size(K,1),:,nex)
    trH =  vec(sum(reshape(md2σ(t1) .* w,size(K,1),:,nex).*(Jac).^2,dims=(1,2)))
    return trH
end

function getDiagHess(N::SingleLayer, w::AbstractArray{R}, Z::AbstractArray{R},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    return sum((md2σ(K*S .+ b) .* w).*(K*Z).^2,dims=1)
end

function getHessmv(N::SingleLayer, w::AbstractArray{R}, Z::AbstractArray{R},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    return K'*(md2σ(K*S .+ b) .* w .* (K*Z))
end
function getHessmv(N::SingleLayer, w::AbstractArray{R}, Z::AbstractArray{R,3},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    t1 = reshape(md2σ(K*S .+ b).*w,size(K,1),1,nex)
    t2 = t1 .* reshape(K*reshape(Z,size(K,2),:),size(K,1),:,nex)
    t2 = K'* reshape(t2,size(K,1),:)
    return reshape(t2,size(K,2),size(Z,2),nex)
end

"""
compute matvec J_S N(S,Θ)*Z
"""
function getJSmv(N::SingleLayer,Z::AbstractArray{R},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (K,b) = Θ
    return mdσ(K*S .+ b) .* (K * Z)
end

"""
compute matvec J_S N(S,Θ)*Z
"""
function getJSmv(N::SingleLayer,Z::AbstractArray{R,3},S::AbstractArray{R,2},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    # K * Z
    KZ = K*reshape(Z,size(K,2),:)
    return reshape(mdσ(K*S .+ b),size(K,1),1,nex) .* reshape(KZ,size(K,1),size(Z,2),nex)
end

"""
compute matvec J_S(J_S N(S,Θ)'*Z(S))

here we use product rule

J_S N(S,Θ)'*dZ + J_S(N(S,Θ)'*Zfix)
"""
function getJSJSTmv(N::SingleLayer,dz::AbstractVector{R},d2z::AbstractArray{R},s::AbstractVector{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (K,b) = Θ
    t1 = K*s + b
    ndσ =  mdσ(t1)
    return K'* ( Diagonal(md2σ(t1) .* dz) +
                               ndσ .* d2z .* ndσ') *K
end

function getJSJSTmv(N::SingleLayer,dZ::AbstractArray{R},d2Z::AbstractArray{R},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (K,b) = Θ
    t1 = K* S .+ b
    (d,nex) = size(t1)

    H1 = getJSJSTmv(N,dZ,S,Θ)

    t2 = mdσ(t1)
    H2 = reshape(t2,d,1,:) .* d2Z .* reshape(t2,1,d,:)
    s1 = K' * reshape(H2,size(K,1),:)
    s1 = permutedims(reshape(s1,size(K,2),size(K,1),nex),(2,1,3))
    s2 = K'*reshape(s1,size(K,1),:)
    return H1 + permutedims(reshape(s2,size(K,2),size(K,2),nex),(2,1,3))
end

function getJSJSTmv(N::SingleLayer,dZ::AbstractArray{R},S::AbstractArray{R,2},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R<:Real
    (K,b) = Θ
    t1 = K* S .+ b
    (d,nex) = size(t1)

    dZ2  = reshape(dZ .* md2σ(t1),size(dZ,1),1,nex)
    KtdZK = K'*reshape(dZ2.*K,size(K,1),:)
    return reshape(KtdZK,size(K,2),size(K,2),nex)
end

function getGradAndHessian(N::SingleLayer,dZ::AbstractArray{R},S::AbstractArray{R,2},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R<:Real
    # Here, no d2Z is give, so we assume it is zero
    (K,b) = Θ
    t1 = K* S .+ b
    (d,nex) = size(t1)

    t2  = reshape(dZ .* md2σ(t1),size(dZ,1),1,nex)
    KtdZK = K'*reshape(t2.*K,size(K,1),:)
    H =  reshape(KtdZK,size(K,2),size(K,2),nex)

    return K'*(mdσ(t1) .* dZ),H
end

function getGradAndHessian(N::SingleLayer,dZ::AbstractArray{R},d2Z::AbstractArray{R},S::AbstractArray{R},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (K,b) = Θ
    t1 = K * S .+ b
    (d,nex) = size(t1)

    t2  = reshape(dZ .* md2σ(t1),size(dZ,1),1,nex)
    KtdZK = K'*reshape(t2.*K,size(K,1),:)
    H1 =  reshape(KtdZK,size(K,2),size(K,2),nex)

    dσt = mdσ(t1)
    t3 = reshape(dσt,d,1,:) .* d2Z .* reshape(dσt,1,d,:)
    s1 = K' * reshape(t3,size(K,1),:)
    s1 = permutedims(reshape(s1,size(K,2),size(K,1),nex),(2,1,3))
    s2 = K'*reshape(s1,size(K,1),:)
    H2  =  permutedims(reshape(s2,size(K,2),size(K,2),nex),(2,1,3))
    return K'*(dσt .* dZ), H1+H2
end

function getJSJSTmv(N::SingleLayer,d2Z::AbstractArray{R,3},S::AbstractArray{R,2},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}}) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    # K * Z
    t1 = reshape(mdσ(K*S .+ b),size(K,1),1,nex)
    t1 = K'*reshape(t1 .* d2Z, size(K,1),:)
    t1 = reshape(t1,size(K,2),size(d2Z,2),nex)
    permutedims!(t1,t1,(2,1,3))
    t1 = K'*reshape(t1 .* d2Z, size(K,1),:)
    t1 = reshape(t1,size(K,2),size(d2Z,2),nex)
    return t1
end

function getJSJSTmv(N::SingleLayer,d2Z::AbstractArray{R,3},S::AbstractArray{R,2},Θ::Tuple{AbstractArray{R,2},AbstractArray{R,1}},hk::R) where R <: Real
    (d,nex) = size(S)
    (K,b) = Θ
    # K * Z
    t1 = reshape(mdσ(K*S .+ b),size(K,1),1,nex)
    t1 = K'*reshape(t1 .* d2Z, size(K,1),:)
    t1 = d2Z + hk .* reshape(t1,size(K,2),size(d2Z,2),nex)
    permutedims!(t1,t1,(2,1,3))
    t1 = K'*reshape(t1 .* d2Z, size(K,1),:)
    t1 = d2Z + hk .* reshape(t1,size(K,2),size(d2Z,2),nex)
    return t1
end


text:
module ImGuiMakieBackend

using GLMakie
using GLMakie.GLFW
using GLMakie.ModernGL
using LibCImGui
using Libdl

Base.convert(::Type{Cint}, x::GLFW.Key) = Cint(x)

const GLFW_GET_CLIPBOARD_TEXT_FUNCPTR = Ref{Ptr{Cvoid}}(C_NULL)
const GLFW_SET_CLIPBOARD_TEXT_FUNCPTR = Ref{Ptr{Cvoid}}(C_NULL)

include("utils.jl")

function c_get(x::Ptr{NTuple{N,T}}, i) where {N,T}
    unsafe_load(Ptr{T}(x), Integer(i)+1)
end

function c_set!(x::Ptr{NTuple{N,T}}, i, v) where {N,T}
    unsafe_store!(Ptr{T}(x), v, Integer(i)+1)
end

include("context.jl")
include("callbacks.jl")
include("font.jl")
include("device.jl")
include("platform.jl")
include("render.jl")
include("interface.jl")

const IMGUI_MAKIE_CONTEXT = Ref{Context}()

function __init__()
    GLFW_GET_CLIPBOARD_TEXT_FUNCPTR[] = dlsym(dlopen(GLFW.libglfw), :glfwGetClipboardString)
    GLFW_SET_CLIPBOARD_TEXT_FUNCPTR[] = dlsym(dlopen(GLFW.libglfw), :glfwSetClipboardString)

    if !isassigned(IMGUI_MAKIE_CONTEXT)
        IMGUI_MAKIE_CONTEXT[] = create_context()
    end
end

get_context() = IMGUI_MAKIE_CONTEXT[]

function vsynced_renderloop(screen)
    while GLMakie.isopen(screen) && !GLMakie.WINDOW_CONFIG.exit_renderloop[]
        GLMakie.pollevents(screen) # GLFW poll

        # create new cimgui frame
        new_frame(get_context())
        igNewFrame()

        screen.render_tick[] = nothing
        if GLMakie.WINDOW_CONFIG.pause_rendering[]
            sleep(0.1)
        else
            # UI
            get_context().ui_func()
            igShowDemoWindow(Ref(true))
            igShowMetricsWindow(Ref(true))
            @show "vsync"
            # render frame
            GLMakie.make_context_current(screen)
            GLMakie.render_frame(screen)

            # render UI frame
            igRender()
            GLMakie.make_context_current(screen)
            render(get_context())
            if unsafe_load(igGetIO().ConfigFlags) & ImGuiConfigFlags_ViewportsEnable == ImGuiConfigFlags_ViewportsEnable
                backup_current_context = GLFW.GetCurrentContext()
                igUpdatePlatformWindows()
                igRenderPlatformWindowsDefault(C_NULL, pointer_from_objref(ctx))
                GLFW.MakeContextCurrent(backup_current_context)
            end

            # swap buffer
            GLFW.SwapBuffers(GLMakie.to_native(screen))
            yield()
        end
    end
end

function fps_renderloop(screen::GLMakie.Screen, framerate=GLMakie.WINDOW_CONFIG.framerate[])
    time_per_frame = 1.0 / framerate
    while GLMakie.isopen(screen) && !WINDOW_CONFIG.exit_renderloop[]
        t = time_ns()
        GLMakie.pollevents(screen) # GLFW poll

        # create new cimgui frame
        new_frame(get_context())
        igNewFrame()

        screen.render_tick[] = nothing
        if WINDOW_CONFIG.pause_rendering[]
            sleep(0.1)
        else
            # UI
            get_context().ui_func()
            igShowDemoWindow(Ref(true))
            igShowMetricsWindow(Ref(true))
            @show "fps"
            # render frame
            GLMakie.make_context_current(screen)
            GLMakie.render_frame(screen)

            # render UI frame
            igRender()
            GLMakie.make_context_current(screen)
            render(get_context())
            if unsafe_load(igGetIO().ConfigFlags) & ImGuiConfigFlags_ViewportsEnable == ImGuiConfigFlags_ViewportsEnable
                backup_current_context = GLFW.GetCurrentContext()
                igUpdatePlatformWindows()
                igRenderPlatformWindowsDefault(C_NULL, pointer_from_objref(ctx))
                GLFW.MakeContextCurrent(backup_current_context)
            end

            # swap buffer
            GLFW.SwapBuffers(GLMakie.to_native(screen))
            t_elapsed = (time_ns() - t) / 1e9
            diff = time_per_frame - t_elapsed
            if diff > 0.001 # can't sleep less than 0.001
                sleep(diff)
            else # if we don't sleep, we still need to yield explicitely to other tasks
                yield()
            end
        end
    end
end

function renderloop(screen; framerate=WINDOW_CONFIG.framerate[])
    isopen(screen) || error("Screen most be open to run renderloop!")
    try
        if GLMakie.WINDOW_CONFIG.vsync[]
            GLFW.SwapInterval(1)
            vsynced_renderloop(screen)
        else
            GLFW.SwapInterval(0)
            fps_renderloop(screen, framerate)
        end
    catch e
        showerror(stderr, e, catch_backtrace())
        println(stderr)
        rethrow(e)
    finally
        shutdown(get_context())
        destroy!(screen)
    end
end

set_makie_renderloop() = GLMakie.set_window_config!(renderloop=renderloop)
export set_makie_renderloop

end


text:
# Substitution Matrices
# =====================
#
# Substitution matrix types and predefined matrices.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioAlignments.jl/blob/master/LICENSE.md

"""
Supertype of substitution matrix.

The required method:

* `Base.getindex(submat, x, y)`: substitution score/cost from `x` to `y`
"""
abstract type AbstractSubstitutionMatrix{S<:Real} end

"""
Substitution matrix.
"""
struct SubstitutionMatrix{T,S} <: AbstractSubstitutionMatrix{S}
    # square substitution matrix
    data::Matrix{S}
    # score is defined or not
    defined::BitMatrix

    function SubstitutionMatrix{T,S}(data::Matrix{S}, defined::BitMatrix) where {T,S}
        @assert size(data) == size(defined)
        @assert size(data, 1) == size(data, 2) == length(BioSymbols.alphabet(T)) - 1
        return new{T,S}(data, defined)
    end

end

function SubstitutionMatrix(::Type{T},
                            submat::AbstractMatrix{S},
                            defined::AbstractMatrix{Bool},
                            default_match::S,
                            default_mismatch::S) where {T,S}
    alpha = alphabet_without_gap(T)
    n = length(alpha)
    data = Matrix{S}(undef, (n, n))
    for x in alpha, y in alpha
        i = index(x)
        j = index(y)
        if defined[i,j]
            data[i,j] = submat[i,j]
        else
            data[i,j] = x == y ? default_match : default_mismatch
        end
    end
    return SubstitutionMatrix{T,S}(data, copy(defined))
end

function SubstitutionMatrix(::Type{T},
                            submat::AbstractMatrix{S},
                            # assume scores of all substitutions are defined
                            defined=trues(size(submat));
                            default_match=S(0), default_mismatch=S(0)) where {T,S}
    return SubstitutionMatrix(T, submat, defined, default_match, default_mismatch)
end

function SubstitutionMatrix(scores::AbstractDict{Tuple{T,T},S};
                            default_match=S(0), default_mismatch=S(0)) where {T,S}
    n = length(BioSymbols.alphabet(T)) - 1
    submat = Matrix{S}(undef, (n, n))
    defined = falses(n, n)
    for ((x, y), score) in scores
        i = index(x)
        j = index(y)
        submat[i,j] = score
        defined[i,j] = true
    end
    return SubstitutionMatrix(T, submat, defined, default_match, default_mismatch)
end

Base.convert(::Type{Matrix}, submat::SubstitutionMatrix) = copy(submat.data)

@inline function Base.getindex(submat::SubstitutionMatrix{T}, x, y) where T
    i = index(convert(T, x))
    j = index(convert(T, y))
    return submat.data[i,j]
end

function Base.setindex!(submat::SubstitutionMatrix{T}, val, x, y) where T
    i = index(convert(T, x))
    j = index(convert(T, y))
    submat.data[i,j] = val
    submat.defined[i,j] = true
    return submat
end

function Base.copy(submat::SubstitutionMatrix{T,S}) where {T,S}
    return SubstitutionMatrix{T,S}(copy(submat.data), copy(submat.defined))
end

Base.minimum(submat::SubstitutionMatrix) = minimum(submat.data)
Base.maximum(submat::SubstitutionMatrix) = maximum(submat.data)

function Base.show(io::IO, submat::SubstitutionMatrix{T,S}) where {T,S}
    alpha = alphabet_without_gap(T)
    n = length(alpha)
    mat = Matrix{String}(undef, (n, n))
    for (i, x) in enumerate(alpha), (j, y) in enumerate(alpha)
        i′ = index(x)
        j′ = index(y)
        mat[i, j] = string(submat.data[i′, j′])
        if !submat.defined[i′, j′]
            mat[i, j] = underline(mat[i, j])
        end
    end

    # add rows and columns
    rows = map(string, alpha)
    cols = vcat(" ", rows)
    mat = hcat(rows, mat)
    mat = vcat(reshape(cols, 1, length(cols)), mat)

    println(io, summary(submat), ':')
    width = maximum(map(x -> length(x), mat))
    for i in 1:n + 1
        for j in 1:n + 1
            print(io, lpad(mat[i, j], width + ifelse(last(mat[i, j]) == '\U0332', 2, 1)))
        end
        println(io)
    end
    print(io, "(underlined values are default ones)")
end

underline(s) = join([string(c, '\U0332') for c in s])

# Return a vector of all symbols of `T` except the gap symbol.
function alphabet_without_gap(::Type{T}) where T
    return filter!(x -> x != BioSymbols.gap(T), collect(BioSymbols.alphabet(T)))
end

# Return the row/column index of `nt`.
function index(nt::BioSymbols.NucleicAcid)
    return convert(Int, nt)
end

# Return the row/column index of `aa`.
function index(aa::BioSymbols.AminoAcid)
    return convert(Int, aa) + 1
end

"""
Dichotomous substitution matrix.
"""
struct DichotomousSubstitutionMatrix{S} <: AbstractSubstitutionMatrix{S}
    match::S
    mismatch::S
end

function DichotomousSubstitutionMatrix(match::Real, mismatch::Real)
    typ = promote_type(typeof(match), typeof(mismatch))
    return DichotomousSubstitutionMatrix{typ}(match, mismatch)
end

function Base.convert(::Type{SubstitutionMatrix{T,S}},
                      submat::DichotomousSubstitutionMatrix) where {T,S}
    n = length(BioSymbols.alphabet(T)) - 1
    data = Matrix{S}(undef, (n, n))
    fill!(data, submat.mismatch)
    data[diagind(data)] .= submat.match
    return SubstitutionMatrix{T,S}(data, trues(n, n))
end

@inline function Base.getindex(submat::DichotomousSubstitutionMatrix, x, y)
    return ifelse(x == y, submat.match, submat.mismatch)
end

function Base.show(io::IO, submat::DichotomousSubstitutionMatrix)
    match = string(submat.match)
    mismatch = string(submat.mismatch)
    w = max(length(match), length(mismatch))
    println(io, typeof(submat), ':')
    # right aligned
    print(io, "     match = ", lpad(match, w), '\n')
    print(io, "  mismatch = ", lpad(mismatch, w))
end

# Predefined Substitution Matrices
# --------------------------------

function load_submat(::Type{T}, name) where T
    return parse_ncbi_submat(T, joinpath(dirname(@__FILE__), "data", "submat", name))
end

function parse_ncbi_submat(::Type{T}, filepath) where T
    scores = Dict{Tuple{T,T},Int}()
    cols = T[]
    open(filepath) do io
        for line in eachline(io)
            line = chomp(line)
            if startswith(line, "#") || isempty(line)
                continue  # skip comments and empty lines
            end
            if isempty(cols)
                for y in eachmatch(r"[A-Z*]", line)
                    @assert length(y.match) == 1
                    push!(cols, convert(T, first(y.match)))
                end
            else
                x = convert(T, first(line))
                ss = collect(eachmatch(r"-?\d+", line))
                @assert length(ss) == length(cols)
                for (y, s) in zip(cols, ss)
                    scores[(x, y)] = parse(Int, s.match)
                end
            end
        end
    end
    return SubstitutionMatrix(scores, default_match=0, default_mismatch=0)
end

"EDNAFULL (or NUC4.4) substitution matrix"
const EDNAFULL = load_submat(BioSymbols.DNA, "NUC.4.4")

"PAM30 substitution matrix"
const PAM30    = load_submat(BioSymbols.AminoAcid, "PAM30")

"PAM70 substitution matrix"
const PAM70    = load_submat(BioSymbols.AminoAcid, "PAM70")

"PAM250 substitution matrix"
const PAM250   = load_submat(BioSymbols.AminoAcid, "PAM250")

"BLOSUM45 substitution matrix"
const BLOSUM45 = load_submat(BioSymbols.AminoAcid, "BLOSUM45")

"BLOSUM50 substitution matrix"
const BLOSUM50 = load_submat(BioSymbols.AminoAcid, "BLOSUM50")

"BLOSUM62 substitution matrix"
const BLOSUM62 = load_submat(BioSymbols.AminoAcid, "BLOSUM62")

"BLOSUM80 substitution matrix"
const BLOSUM80 = load_submat(BioSymbols.AminoAcid, "BLOSUM80")

"BLOSUM90 substitution matrix"
const BLOSUM90 = load_submat(BioSymbols.AminoAcid, "BLOSUM90")


text:
using JuMP, EAGO, Gurobi, CSV, DataFrames, Ipopt

m = Model(with_optimizer(EAGO.Optimizer, relaxed_optimizer = Gurobi.Optimizer(OutputFlag=0)))

# ----- Variables ----- #
x_Idx = Any[1, 2, 3, 4, 5, 6, 7, 8]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[1], 100.0)
JuMP.set_upper_bound(x[1], 10000.0)
JuMP.set_lower_bound(x[2], 1000.0)
JuMP.set_upper_bound(x[2], 10000.0)
JuMP.set_lower_bound(x[3], 1000.0)
JuMP.set_upper_bound(x[3], 10000.0)
JuMP.set_lower_bound(x[4], 10.0)
JuMP.set_upper_bound(x[4], 1000.0)
JuMP.set_lower_bound(x[5], 10.0)
JuMP.set_upper_bound(x[5], 1000.0)
JuMP.set_lower_bound(x[6], 10.0)
JuMP.set_upper_bound(x[6], 1000.0)
JuMP.set_lower_bound(x[7], 10.0)
JuMP.set_upper_bound(x[7], 1000.0)
JuMP.set_lower_bound(x[8], 10.0)
JuMP.set_upper_bound(x[8], 1000.0)

# ----- Constraints ----- #
@NLconstraint(m, e2, 833.33252*x[4]/x[1]/x[6]+100/x[6]-83333.333/(x[1]*x[6]) <= 1.0)
@NLconstraint(m, e3, 1250*x[5]/x[2]/x[7]+x[4]/x[7]-1250*x[4]/x[2]/x[7] <= 1.0)
@NLconstraint(m, e4, 1250000/(x[3]*x[8])+x[5]/x[8]-2500*x[5]/x[3]/x[8] <= 1.0)
@constraint(m, e5, 0.0025*x[4]+0.0025*x[6] <= 1.0)
@constraint(m, e6, -0.0025*x[4]+0.0025*x[5]+0.0025*x[7] <= 1.0)
@constraint(m, e7, -0.01*x[5]+0.01*x[8] <= 1.0)

# ----- Objective ----- #
@objective(m, Min, x[1] + x[2] + x[3])

JuMP.optimize!(m)

backend_opt = backend(m).optimizer.model.optimizer
last_relaxed_opt = backend_opt.relaxed_optimizer
run_time = backend_opt._run_time
println("run time: $run_time")


text:
# Benchmark the computation of the area of a spherical triangle
#
# Lestringant uses this to calculate twists, but the trig calls are
# cripplingly slow. Here we investigate an allocating version, a
# non-allocating version, and a version accelerated with IntelVectorMath
#
# Conclusion: the IVM non-allocating version wins hands down, and is
# almost 3x faster than its competitor. However, it doesn't permit forward-
# mode automatic differentiation because IVM functions only evaluate on
# Float64 and Float32 types. Overall twist calculation via 2 calls
# to spherical excess with IVM is no faster than twist calculations
# based on parallel transport, and the latter is type-generic.
#
# A possible way forward is to overload IVM.atan!(a, b, c) to operate
# on ForwardDiff.Dual types. That would allow at least fast HVPs,
# and is worth investigating. Note that
#
# atan(y + δy, x + δx) = atan(y, x) + (δy*x - δx*y)/(x^2 + y^2)

using DiscreteElasticRods
DER = DiscreteElasticRods
using BenchmarkTools

# Full formula form Lestringant 2020
function spherical_excess_reference(t1, t2, t3)

    ca = dot3(t2,t3)
    cb = dot3(t3,t1)
    cc = dot3(t1,t2)

    sa = norm3(cross3(t2,t3))
    sb = norm3(cross3(t3,t1))
    sc = norm3(cross3(t1,t2))

    cA = (ca - cb.*cc)./(sb.*sc)
    cB = (cb - cc.*ca)./(sc.*sa)
    cC = (cc - ca.*cb)./(sa.*sb)

    V = dot3(cross3(t1,t2),t3)
    sA = V./(sb.*sc)
    sB = V./(sc.*sa)
    sC = V./(sa.*sb)

    atan.(sA,cA) + atan.(sB,cB) + atan.(sC,cC) - π
end

# Removing common denominators in atan calls
function spherical_excess(t1, t2, t3)

    ca = DER.dot3(t2,t3)
    cb = DER.dot3(t3,t1)
    cc = DER.dot3(t1,t2)
    V = DER.dot3(DER.cross3(t1,t2),t3)

    atan.(V, ca .- cb.*cc) .+
    atan.(V, cb .- cc.*ca) .+
    atan.(V, cc .- ca.*cb) .- π
end

# non-allocating vector triple product
function triple3!(out, a, b, c)
    @inbounds for i in 1:size(out,1)
        out[i] = c[i,1]*(a[i,2]*b[i,3] - a[i,3]*b[i,2]) +
                 c[i,2]*(a[i,3]*b[i,1] - a[i,1]*b[i,3]) +
                 c[i,3]*(a[i,1]*b[i,2] - a[i,2]*b[i,1])
   end
end

# Non-allocating spherical excess
# (ns + 1) x 4 cache
# 90% of this function's runtime is atan calls
function spherical_excess!(γ, t1, t2, t3, cache)

    ca = @view cache[:,1]
    cb = @view cache[:,2]
    cc = @view cache[:,3]
    V = @view cache[:,4]

    triple3!(V, t1, t2, t3)
    DER.dot3!(ca, t2, t3)
    DER.dot3!(cb, t3, t1)
    DER.dot3!(cc, t1, t2)

    γ .= @. atan(V, ca-cb*cc) + atan(V, cb-cc*ca) + atan(V, cc-ca*cb) - π
end

# IVM acceleration
# Note that IVM.atan!(a,b,c) can only operator on Array{Float64,N}, so
# the cache cannot use views or subarrays of any kind
using IntelVectorMath

struct spherical_excess_IVM_cache{T<:Real}
    ca::Vector{T}
    cb::Vector{T}
    cc::Vector{T}
    V::Vector{T}
    tmp1::Vector{T}
    tmp2::Vector{T}
end

function allocate_seIVM(T::Type, N::Int)
    ca = Vector{T}(undef, N)
    cb = Vector{T}(undef, N)
    cc = Vector{T}(undef, N)
    tmp1 = Vector{T}(undef, N)
    tmp2 = Vector{T}(undef, N)
    V = Vector{T}(undef, N)
    spherical_excess_IVM_cache(ca, cb, cc, V, tmp1, tmp2)
end

function spherical_excess_IVM!(γ, t1, t2, t3, cache::spherical_excess_IVM_cache)

    ca, cb, cc = cache.ca, cache.cb, cache.cc
    V, tmp1, tmp2 =  cache.V, cache.tmp1, cache.tmp2

    triple3!(V, t1, t2, t3)
    DER.dot3!(ca, t2, t3)
    DER.dot3!(cb, t3, t1)
    DER.dot3!(cc, t1, t2)

    γ .= -π
    tmp1 .= ca .- cb.*cc
    IVM.atan!(tmp2, V, tmp1)
    γ .+= tmp2

    tmp1 .= cb .- cc.*ca
    IVM.atan!(tmp2, V, tmp1)
    γ .+= tmp2

    tmp1 .= cc .- ca.*cb
    IVM.atan!(tmp2, V, tmp1)
    γ .+= tmp2
end

# Benchmarking
N = 500

# Pure atan comparison
e = rand(N)
d = rand(N)
f = similar(e)
@btime IVM.atan!($f, $e, $d)
@btime $f .= atan.($e, $d)

# Non-allocating triple product
a = rand(N,3)
b = rand(N,3)
c = rand(N,3)
out = zeros(N)
@btime triple3!($out, $a, $b, $c)

# Allocating vs non-allocating vs IVM
cache_sphere = zeros(N,6)
cache_sphere_IVM = allocate_seIVM(Float64, N)
@btime spherical_excess($a, $b, $c)
@btime spherical_excess!($out, $a, $b, $c, $cache_sphere)
@btime spherical_excess_IVM!($out, $a, $b, $c, $cache_sphere_IVM)


text:
# make sure we have the required imports
import Base: iterate, IteratorSize, IteratorEltype, eltype, length

abstract type DependencyGraph end

@inline dependents(dep_graph::DependencyGraph, j) = dependents(dg_iterator(dep_graph), dep_graph, j)

abstract type DGIterStyle end

struct DGView <: DGIterStyle end
struct DGLazy <: DGIterStyle end

# Implementation: range tuples + single partitioned vector

struct DGVector{Iter<:DGIterStyle} <: DependencyGraph
  deps_range::Vector{Tuple{Int,Int}}
  deps_value::Vector{Int}
end

## constructors

function rxnrxn_depgraph(::Iter, model::Network) where Iter <: DGIterStyle
  num_reactions = number_reactions(model)
  reactions = reaction_list(model)
  dg = model.dep_graph

  indexmap = OrderedDict(key => i for (i, key) in enumerate(keys(reactions)))

  deps_range = Vector{Tuple{Int,Int}}(undef, num_reactions)
  deps_value  = Int[]

  idx = 1
  for (key, j) in indexmap
    start = idx
    for dep in dg[key]
      push!(deps_value, indexmap[dep])
      idx += 1
    end
    stop = idx - 1

    deps_range[j] = (start, stop)
  end

  return DGVector{Iter}(deps_range, deps_value)
end

function spcrxn_depgraph(::Iter, model::Network) where Iter <: DGIterStyle
  num_reactions = number_reactions(model)
  num_species = number_species(model)

  species   = species_list(model)
  reactions = reaction_list(model)

  indexmap = OrderedDict(key => i for (i, key) in enumerate(keys(species)))

  deps_range = Vector{Tuple{Int,Int}}(undef, num_species)
  deps_value  = Int[]

  idx = 1
  for (key, i) in indexmap
    start = idx
    for (j, rxn) in enumerate(values(reactions))
      if haskey(rxn.reactants, key)
        push!(deps_value, j)
        idx += 1
      end
    end
    stop = idx - 1
    deps_range[i] = (start, stop)
  end

  return DGVector{Iter}(deps_range, deps_value)
end

function rxnspc_depgraph(::Iter, model::Network) where Iter <: DGIterStyle
  num_reactions = number_reactions(model)
  num_species = number_species(model)

  species   = species_list(model)
  reactions = reaction_list(model)

  indexmap = OrderedDict(key => i for (i, key) in enumerate(keys(species)))

  deps_range = Vector{Tuple{Int,Int}}(undef, num_reactions)
  deps_value  = Int[]

  idx = 1

  for (j, rxn) in enumerate(values(reactions))
    start = idx
    for (key, i) in indexmap
      if haskey(rxn.products, key)
        push!(deps_value, i)
        idx += 1
      end
    end
    stop = idx - 1
    deps_range[j] = (start, stop)
  end

  return DGVector{Iter}(deps_range, deps_value)
end

## iteration using views

@inline function dependents_view(dg::DGVector, j)
  @inbounds (start, stop) = dg.deps_range[j]

  return view(dg.deps_value, start:stop)
end

## iteration using laziness

struct DGVectorIterator
  deps::Vector{Int}
  start::Int
  stop::Int
end

Base.IteratorSize(::DGVectorIterator) = Base.HasLength()
Base.IteratorEltype(::DGVectorIterator) = Base.HasEltype()

Base.length(iter::DGVectorIterator) = iter.stop - iter.start + 1
Base.eltype(iter::DGVectorIterator) = eltype(iter.deps)

Base.iterate(iter::DGVectorIterator) = iterate(iter, iter.start)

function Base.iterate(iter::DGVectorIterator, state)
  return state > iter.stop ? nothing : (iter.deps[state], state + 1)
end

@inline function dependents_lazy(dg::DGVector, j)
  @inbounds (start, stop) = dg.deps_range[j]

  return DGVectorIterator(dg.deps_value, start, stop)
end

## iteration protocol

@inline dg_iterator(dg::DGVector{Iter}) where Iter <: DGIterStyle = Iter()

@inline dependents(::DGView, dep_graph::DGVector, j) = dependents_view(dep_graph, j)
@inline dependents(::DGLazy, dep_graph::DGVector, j) = dependents_lazy(dep_graph, j)


text:
function extract_edit_menu(robj::RenderObject, edit_screen, isvisible)
    data = map(robj.uniforms) do kv
        label = replace(string(kv[1]), "_", " ")*":"
        label => kv[2]
    end
    extract_edit_menu(data, edit_screen, isvisible)
end

makesignal2(v) = Signal(v)
makesignal2(s::Signal) = s
makesignal2(v::GPUArray) = v
keytype(::Dict{K, T}) where {K, T} = K
keytype(::AbstractVector{Pair{K, T}}) where {K, T} = K
keytype(::AbstractVector{Pair{K}}) where {K} = K
keytype(x) = Any

function extract_edit_menu(
        edit_dict, edit_screen, isvisible = Signal(true);
        icon_size = 32, knob_scale = 1.6mm
    )

    screen_w = widths(edit_screen)[1] - round(Int, 10mm)
    pos = 1mm
    widget_text = 4mm
    visses = Pair{String, Context{GLAbstraction.DeviceUnit}}[]
    K = keytype(edit_dict)
    signal_dict = Dict{K, Any}()
    for (k, v) in edit_dict
        label = String(string(k))
        s = makesignal2(v)
        if applicable(widget, s, edit_screen)
            vis, sig = widget(s, edit_screen,
                visible = isvisible,
                text_scale = widget_text,
                area = (screen_w, icon_size),
                knob_scale = knob_scale
            )
            push!(visses, label => vis)
            signal_dict[k] = sig
        end
    end
    visualize(visses), signal_dict
end


text:
using Test
using Bijectors, ForwardDiff, LinearAlgebra
using Random: seed!

seed!(1)

@testset "PlanarLayer" begin
    for i in 1:4
        flow = PlanarLayer(2)
        z = randn(2, 20)
        forward_diff = log(abs(det(ForwardDiff.jacobian(t -> flow(t), z))))
        our_method = sum(forward(flow, z).logabsdetjac)

        @test our_method ≈ forward_diff
        @test inv(flow)(flow(z)) ≈ z rtol=0.2
        @test (inv(flow) ∘ flow)(z) ≈ z rtol=0.2
    end

    w = ones(10, 1)
    u = zeros(10, 1)
    b = ones(1)
    flow = PlanarLayer(w, u, b)
    z = ones(10, 100)
    @test inv(flow)(flow(z)) ≈ z
end

@testset "RadialLayer" begin
    for i in 1:4
        flow = RadialLayer(2)
        z = randn(2, 20)
        forward_diff = log(abs(det(ForwardDiff.jacobian(t -> flow(t), z))))
        our_method = sum(forward(flow, z).logabsdetjac)

        @test our_method ≈ forward_diff
        @test inv(flow)(flow(z)) ≈ z rtol=0.2
        @test (inv(flow) ∘ flow)(z) ≈ z rtol=0.2
    end

    α_ = ones(1)
    β = ones(1)
    z_0 = zeros(10, 1)
    z = ones(10, 100)
    flow = RadialLayer(α_, β, z_0)
    @test inv(flow)(flow(z)) ≈ z
end

@testset "Flows" begin
    d = MvNormal(zeros(2), ones(2))
    b = PlanarLayer(2)
    flow = transformed(d, b)  # <= Radial flow

    y = rand(flow)
    @test logpdf(flow, y) != 0.0

    x = rand(d)
    y = flow.transform(x)
    res = forward(flow.transform, x)
    lp = logpdf_forward(flow, x, res.logabsdetjac)

    @test res.rv ≈ y
    @test logpdf(flow, y) ≈ lp rtol=0.1

    # flow with unconstrained-to-constrained
    d1 = Beta()
    b1 = inv(bijector(d1))
    d2 = InverseGamma()
    b2 = inv(bijector(d2))

    x = rand(d) .+ 10
    y = b(x)

    sb = stack(b1, b1)
    @test all((sb ∘ b)(x) .≤ 1.0)

    sb = stack(b1, b2)
    cb = (sb ∘ b)
    y = cb(x)
    @test (0 ≤ y[1] ≤ 1.0) && (0 < y[2])
end


text:
# Test functions for tracing in different kinds of analytical fields.

include("utility/dipole.jl")

"""
    test_trace_asymptote(issingle=false)

Test streamline tracing by plotting vectors and associated streamlines through a
simple velocity field where Vx=x, Vy=-y. Support for single and double precision
data.
"""
function test_trace_asymptote(issingle=false)

   # Start by creating a velocity vector field.
   if issingle
      xmax, ymax = 200.0f0, 20.0f0
      x = -10.0f0:0.25f0:xmax
      y = -10.0f0:0.25f0:ymax
   else
      xmax, ymax = 200.0, 20.0
      x = -10.0:0.25:xmax
      y = -10.0:0.25:ymax
   end

   xgrid = [i for j in y, i in x]
   ygrid = [j for j in y, i in x]

   if issingle
      vx = xgrid * 1.0f0
      vy = ygrid * -1.0f0

      xstart = 1.0f0
      ystart = 10.0f0
   else
      vx = xgrid * 1.0
      vy = ygrid * -1.0

      xstart = 1.0
      ystart = 10.0
   end

   ds = 0.1
   gridtype = "meshgrid"

   x1, y1 = trace(vx, vy, xstart, ystart, x, y; alg=RK4(), ds, gridtype)
   x2, y2 = trace(vx, vy, xstart, ystart, x, y; alg=Euler(), ds, gridtype)

   # analytical solution const = x*y
   c = xstart * ystart

   if issingle
      if isapprox(x1[end]*y1[end], c, atol=1e-2) && isapprox(x2[end]*y2[end], c, atol=0.13)
         return true
      else
         return false
      end
   else
      if isapprox(x1[end]*y1[end], c, atol=1e-5) && isapprox(x2[end]*y2[end], c, atol=0.12)
         return true
      else
         return false
      end
   end
end

"Trace field lines through a Earth like magnetic dipole field."
function test_trace_dipole()

   # Start by creating a field of unit vectors.
   x = -100.0:5.0:101.0
   y = -100.0:5.0:101.0

   bx, by = b_hat(x, y)

   # Trace through this field.
   xstart = 10.0
   ystart = 25.0
   ds = 0.1
   gridtype = "meshgrid"

   x1, y1 = trace(bx, by, xstart, ystart, x, y; alg=RK4(), ds, gridtype)
   x2, y2 = trace(bx, by, xstart, ystart, x, y; alg=Euler(), ds, gridtype)

   if length(x1) == 260 && isapprox(x1[end]*y1[end], 7462.42688, atol=1e-5) &&
      length(x2) == 262 && isapprox(x2[end]*y2[end], 7601.60590, atol=1e-5)
      return true
   else
      return false
   end
end


text:
using Documenter, Genie

makedocs()

max_stars_repo_name: essenciary/genie-todo-mvc

max_issues_repo_name: essenciary/genie-todo-mvc

max_forks_repo_name: essenciary/genie-todo-mvc


text:
@testset "Jeffrey Problems" begin
    (p, q, r, x, ) = @variables p q r x

    #= ::Package:: =#

    #= ::Title:: =#
    #=David*Jeffrey - Rectifying*Transformations*for*Trig*Integration*(1997)=#

    #= ::Subsection::Closed:: =#
    #=Problem*(1.2)=#

    @test_int [3/(5 - 4*cos(x)), x, 2, x + 2*atan(sin(x)/(2 - cos(x)))]

    #= ::Subsection::Closed:: =#
    #=Problem*(1.4)=#

    @test_int [(cos(x) + 2*sin(x) + 1)/(cos(x)^2 - 2*sin(x)*cos(x) + 2*sin(x) + 3), x, -43, -atan((2*cos(x) - sin(x))/(2 + sin(x)))]

    #= ::Subsection::Closed:: =#
    #=Problem*(1.5)=#

    @test_int [(2 + 5*sin(x) + cos(x))/(4*cos(x) + sin(x)*cos(x) - 2*sin(x) - 2*sin(x)^2), x, -25, -log(1 - 3*cos(x) + sin(x)) + log(3 + cos(x) + sin(x))]

    #= ::Subsection::Closed:: =#
    #=Problem*(3.3)=#

    @test_int [(7*cos(x) + 2*sin(x) + 3)/(3*cos(x)^2 - sin(x)*cos(x) + 4*cos(x) - 5*sin(x) + 1), x, -32, -log(1 + cos(x) - 2*sin(x)) + log(3 + cos(x) + sin(x))]

    #= ::Subsection::Closed:: =#
    #=Problem =#

    @test_int [(5*cos(x)^2 + 4*cos(x) - 1)/(4*cos(x)^3 - 3*cos(x)^2 - 4*cos(x) - 1), x, -2, x - 2*atan(sin(x)/(3 + cos(x))) - 2*atan((3*sin(x) + 7*cos(x)*sin(x))/(1 + 2*cos(x) + 5*cos(x)^2))]

    #= ::Subsection::Closed:: =#
    #=Problem =#

    @test_int [(7*cos(x)^2 + 2*cos(x) - 5)/(4*cos(x)^3 - 9*cos(x)^2 + 2*cos(x) - 1), x, -2, x - 2*atan((2*cos(x)*sin(x))/(1 - cos(x) + 2*cos(x)^2))]

    #= ::Subsection::Closed:: =#
    #=Problem*(3.4)=#

    @test_int [3/(5 + 4*sin(x)), x, 2, x + 2*atan(cos(x)/(2 + sin(x)))]

    #= ::Subsection::Closed:: =#
    #=Problem*(3.6)=#

    @test_int [2/(1 + cos(x)^2), x, 3, sqrt(2)*x - sqrt(2)*atan((cos(x)*sin(x))/(1 + sqrt(2) + cos(x)^2))]

    #= ::Subsection::Closed:: =#
    #=Problem*(3.8)=#

    @test_int [1/(p + q*cos(x) + r*sin(x)), x, 3, (2*atan((r + (p - q)*tan(x/2))/sqrt(p^2 - q^2 - r^2)))/sqrt(p^2 - q^2 - r^2)]
    end

max_stars_repo_path: test/0 Independent test suites/Jeffrey Problems.jl


text:
using HDF5
include("./karger.jl")
using .Karger

# USAGE: julia src/julia/calculate_power_watershed.jl filename
# where filename is just the name (without path or .h5).
# The file has to be an HDF5 file with fields n, edges, weights and seeds. The output
# will be written in a file with the same name in the "power_watershed" directory

graph_path = "results/graphs/" * ARGS[1] * ".h5"
result_path = "results/power_watershed/" * ARGS[1] * ".h5"

n = h5read(graph_path, "n")
edges = h5read(graph_path, "edges")
weights = h5read(graph_path, "weights")
min_weight = minimum(weights)
max_weight = maximum(weights)
# normalize weights to range [0, 1]
weights ./= (max_weight - min_weight)
weights .+= min_weight
# discretize to 8bit
weights .*= 255
weights .= round.(weights)

seeds = h5read(graph_path, "seeds")
g = Graph(n, transpose(edges) .+ 1, weights)
mkpath(dirname(result_path))
cp(graph_path, result_path, force=true)

pots = power_watershed_multi(g, seeds)
segmentation = Vector{Int}(undef, n)
max_pots = zeros(n)
for (key, pot) in pots
    h5write(result_path, "potential/" * string(key), pot)
    for i in 1:n
        if pot[i] > max_pots[i]
            segmentation[i] = key
            max_pots[i] = pot[i]
        end
    end
end
h5write(result_path, "segmentation", segmentation)


text:
include("../funcitons/layer_backend/linear.jl")

mutable struct Linear <: Layer
    W :: Union{Nothing, Variable}
    b :: Union{Nothing, Variable}
    in_dim :: Union{Nothing, Int}
    out_dim :: Int
    initial_method :: Union{String, Function}
    function Linear(out_dim; in_dim=nothing, initial_method="xavier")
        if in_dim isa Nothing
            b = Variable(zeros(1, out_dim), name="b")
            layer = new(nothing, b, in_dim, out_dim, initial_method)
            return layer
        else
            if initial_method == "xavier"
                W = xavier(in_dim, out_dim)
            elseif initial_method == "he"
                W = he(in_dim, out_dim)
            else initial_method isa Function
                W = initial_method(in_dim, out_dim)
            end
            W = Variable(W, name="W")
            b = Variable(zeros(1, out_dim), name="b")
            layer = new(W, b, in_dim, out_dim, initial_method)
            return layer
        end
    end
end

parameters(linear::Linear) = (linear.W, linear.b)

function forward(layer::Linear, x)
    if layer.W isa Nothing
        layer.in_dim = size(x.values)[2]
        if layer.initial_method == "xavier"
            W = xavier(layer.in_dim, layer.out_dim)
        elseif layer.initial_method == "he"
            W = he(layer.in_dim, layer.out_dim)
        else layer.initial_method isa Function
            W = layer.initial_method(layer.in_dim, layer.out_dim)
        end
        layer.W = Variable(W, name="W")
    end
    return linear(x, layer.W, layer.b)
end