text:
using Dierckx
# -------------- Vectors and Paths  -----------------

export Path

struct Path{T, NI}
   x::Vector{T}
   valNI::Type{Val{NI}}
end
Path(x::Vector) = Path(x,Val{length(x)})

# vectorise an object of type Path
function Base.vec{T}(x::Vector{T})
   return cat(1, x...)
end

# convert an object of type Path to type Vector
function Base.convert{T, NI}(::Type{Path{T,NI}}, X::Vector{<: AbstractFloat})
   return [ X[(n-1)*(length(X)÷NI)+1 : n*(length(X)÷NI)] for n=1:NI]
end

export serial, palindrome

# ------------- Path Traversing ------------

"""
energy evaluations along a path in the String/NEB-type methods can be performed
in a 'serial' manner or in a 'palindrome' manner:

* `serial` : traverse path always in the same direction.
* `palindrome` : the order of computing the energy of the images along a path is reversed after each iteration.

### Parameters:
* `direction` : order of traversing the images along the path
"""
@with_kw type serial
   direction = (M, nit) -> 1:M
end

@with_kw type palindrome
   direction = (M, nit) -> M-mod(nit,2)*(M-1):2*mod(nit,2)-1:M-mod(nit+1,2)*(M-1)
end

# ------------- Reparametrisation of Paths ------------

"""
parametrisation functions of paths for String/NEB-type methods
"""
# interpolate and reparametrise path
function parametrise!{T}(dxds::Vector{T}, x::Vector{T}, ds::T; parametrisation=linspace(0.,1.,length(x)))

   param = [0; [sum(ds[1:i]) for i in 1:length(ds)]]
   param /= param[end]; param[end] = 1.

   # cubic splines interpolation
   S = [Spline1D(param, [x[i][j] for i=1:length(x)], w = ones(length(x)), k = 3, bc = "error") for j=1:length(x[1])]

   # reparametrise
   xref = [[Sj(s) for s in parametrisation] for Sj in S ]
   dxdsref = [[derivative(Sj, s) for s in parametrisation] for Sj in S]
   d²xds² = [[derivative(Sj, s, nu=2) for Sj in S] for s in parametrisation]

   # return parametrised variables
   x_ = cat(2, xref...)
   dxds_ = cat(2, dxdsref...)
   [x[i] = x_[i,:] for i=1:length(x)]
   [dxds[i] = dxds_[i,:] for i=1:length(x)]
   return d²xds²
end

# distribute images uniformly with respect to P
function redistribute{T,NI}(X::Vector{Float64}, path_type::Type{Path{T,NI}}, precon, precon_scheme)

   x = convert(path_type, X)
   t = deepcopy(x)

   Np = length(precon);
   function P(i) return precon[mod(i-1,Np)+1, 1]; end
   function P(i, j) return precon[mod(i-1,Np)+1, mod(j-1,Np)+1]; end

   ds = [dist(precon_scheme, P, x, i) for i=1:length(x)-1]
   parametrise!(t, x, ds)

   return vec(x)
end

# -------------- Preconditioning for Paths  -----------------

export localPrecon, globalPrecon

"""
the two preconditioning schemes implemented for String/NEB-type methods

`localPrecon` : the preconditioner acts as a coordinate transformation of the
state space.
`globalPrecon` : the preconditioner is applied to the force directly, each point
along the path is preconditioned independently.

### Parameters:
* `precon` : preconditioner
* `precon_prep!` : update function for preconditioner
* `distance` : distance metric

### Functionality:
* `dist` : distance of neighbouring images along path
* `dot_P` : preconditioned dot product of associated images in two paths x, y
* `point_norm` : local P-norm
* `proj_grad` : projected gradient
* `forcing` : reference vector of preconditioned forces along path
* `elastic_force` : Hooke's law elastic focres along path
* `maxres` : residual error
"""

@with_kw type localPrecon
   precon = I
   precon_prep! = (P, x) -> P
   distance = (P, x1, x2) -> norm(P, x2 - x1)
end

@with_kw type globalPrecon
   precon = I
   precon_prep! = (P, x) -> P
   distance = (P, x1, x2) -> norm(x2 - x1)
end

dist(precon_scheme::localPrecon, P, x, i) = precon_scheme.distance(0.5*(P(i)+P(i+1)), x[i], x[i+1])
dist(precon_scheme::globalPrecon, P, x, i) = precon_scheme.distance(x[i], x[i+1])

dot_P(precon_scheme::localPrecon, x, P, y) = sum([dot(x[i], P(i), y[i]) for i=1:length(x)])
dot_P(precon_scheme::globalPrecon, x, P, y) = sum([dot(x[i], y[i]) for i=1:length(x)])

point_norm(precon_scheme::localPrecon, P, dxds) = [ 1; [norm(P(i), dxds[i])
                                                    for i=2:length(dxds)-1]; 1 ]
point_norm(precon_scheme::globalPrecon, P, dxds) = [norm(dxds[i]) for i=1:length(dxds)]

proj_grad(precon_scheme::localPrecon, P, ∇E, dxds) = -[P(i) \ ∇E[i] - dot(∇E[i],dxds[i])*dxds[i] for i=1:length(dxds)]
proj_grad(precon_scheme::globalPrecon, P, ∇E, dxds) = vec(-[∇E[i] - dot(∇E[i],dxds[i])*dxds[i] for i=1:length(dxds)])

forcing(precon_scheme::localPrecon, P, neg∇E⟂) = return vec(neg∇E⟂)
forcing(precon_scheme::globalPrecon, P, neg∇E⟂) = vec(P) \ neg∇E⟂

function elastic_force(precon_scheme::localPrecon, P, κ, dxds, d²xds²)
    return - [ [zeros(dxds[1])];
               κ*[dot(d²xds²[i], P(i), dxds[i]) * dxds[i] for i=2:length(dxds)-1];
               [zeros(dxds[1])] ]
end
function elastic_force(precon_scheme::globalPrecon, P, κ, dxds, d²xds²)
    return vec(-[ [zeros(dxds[1])];
               κ*[dot(d²xds²[i], dxds[i]) * dxds[i] for i=2:N-1];
               [zeros(dxds[1])] ])
end

maxres(precon_scheme::localPrecon, P, ∇E⟂) =  maximum([norm(P(i)*∇E⟂[i],Inf)
                                                for i = 1:length(∇E⟂)])
maxres(precon_scheme::globalPrecon, P, ∇E⟂) = vecnorm(∇E⟂, Inf)


text:
module GitLinks

    import SimpleLockFiles
    import SimpleLockFiles: acquire_lock

    include("Types/GitLink.jl")

    include("Utils/rand_str.jl")
    include("Utils/toml_utils.jl")
    include("Utils/runcmd.jl")

    include("GitUtils/check_gitdir.jl")
    include("GitUtils/repo_format.jl")
    include("GitUtils/remote_HEAD_hash.jl")
    include("GitUtils/config.jl")
    include("GitUtils/commits_behind_remote.jl")
    include("GitUtils/HEAD_hash.jl")
    include("GitUtils/curr_branch.jl")
    include("GitUtils/curr_remotes.jl")
    include("GitUtils/fetch.jl")
    include("GitUtils/hard_pull.jl")
    include("GitUtils/is_up_to_day.jl")
    include("GitUtils/list_commits.jl")
    include("GitUtils/nuke_remote.jl")
    include("GitUtils/ready_to_push.jl")
    include("GitUtils/soft_push.jl")
    include("GitUtils/url_from_file.jl")

    include("TreeStruct/dir_and_files.jl")
    include("TreeStruct/utils.jl")

    include("Lock_system/lock_file.jl")

    include("DevLand/create_local_upstream.jl")

    include("Server/instantiate.jl")
    include("Server/loop_frec.jl")
    include("Server/is_pull_required.jl")
    include("Server/is_push_required.jl")
    include("Server/stage_sync.jl")
    include("Server/signals.jl")
    include("Server/run_sync_loop.jl")
    include("Server/sync_link.jl")
    include("Server/tokens.jl")

    include("Client/readwdir.jl")
    include("Client/download.jl")
    include("Client/writewdir.jl")
    include("Client/events.jl")
    include("Client/stage.jl")
    include("Client/upload.jl")
    include("Client/ping.jl")
    include("Client/git_status.jl")

    export GitLink, instantiate, stage, readwdir, writewdir
    export repo_dir, state_dir
    export create_local_upstream
    export waitfor_pull, waitfor_stage, waitfor_push
    export if_pull, if_stage, if_push
    export is_push_required, is_pull_required
    export sync_link, upload, download, ping, git_status

end

ext: jl


text:
function qfunc end
function wigner end
function coherentspinstate end
function qfuncsu2 end
function wignersu2 end


text:
module DFControl
const DFC = DFControl
export DFC

const CONFIG_DIR = occursin("cache", first(Base.DEPOT_PATH)) ?
                   abspath(Base.DEPOT_PATH[2], "config", "DFControl") :
                   abspath(Base.DEPOT_PATH[1], "config", "DFControl")
const DEPS_DIR = joinpath(dirname(@__DIR__), "deps")

const PYTHONPATH = Sys.iswindows() ? joinpath(DEPS_DIR, "python2", "python") :
                   joinpath(dirname(@__DIR__), "deps", "python2", "bin", "python")

const CIF2CELLPATH = Sys.iswindows() ?
                     joinpath(DEPS_DIR, "python2", "Scripts", "cif2cell") :
                     joinpath(dirname(@__DIR__), "deps", "python2", "bin", "cif2cell")

config_path(path...) = joinpath(CONFIG_DIR, path...)

using LinearAlgebra
using Reexport
@reexport using StaticArrays

const Point{N,F} = SVector{N,F}
const Point3{F}  = SVector{3,F}
const Vec3{F}    = SVector{3,F}
const Mat3{F}    = SMatrix{3,3,F,9}
const Mat4{F}    = SMatrix{4,4,F,16}
const Vec{N,T}   = SVector{N,T}

Point{N,T}(x::T) where {N,T} = Point{N,T}(x, x, x)

function Base.convert(::Type{Point3{T}}, x::Vector{T}) where {T<:AbstractFloat}
    return Point3{T}(x[1], x[2], x[3])
end

export Point, Point3, Vec3, Mat3, Mat4

using Parameters, StructTypes, Dates
include("types.jl")
include("utils.jl")

# Common functions

include("Structures/Structures.jl")
include("Calculations/Calculations.jl")
include("Jobs/Jobs.jl")
include("FileIO/FileIO.jl")

include("Servers/Servers.jl")
include("Service/Service.jl")
include("Resource/Resource.jl")
include("Client/Client.jl")

include("Display/Display.jl")

@reexport using .Client

function __init__()
    Servers.maybe_create_localhost()
    # Client.maybe_start_server("localhost")
    # if !haskey(ENV, "IS_DAEMON")
    #     init_daemon()
    # else
    #     global_logger(DFControl.daemon_logger())
    # end
    return
end

# if ccall(:jl_generating_output, Cint, ()) == 1
#     Job(joinpath(dirname(dirname(pathof(DFControl))),"test/testassets/reference_job"))
# end
end


text:
################################################################################
## Linear Regression
##   y ~ MvNormal(X * (beta0 .* gamma), 1)
##   gamma ~ DiscreteUniform(0, 1)
################################################################################

using MCPhylo

## Data
n, p = 25, 10
X = randn(n, p)
beta0 = randn(p)
gamma0 = rand(0:1, p)
y = X * (beta0 .* gamma0) + randn(n)

## Log-transformed Posterior(gamma) + Constant
logf = function(gamma::DenseVector)
  logpdf(MvNormal(X * (beta0 .* gamma), 1.0), y)
end

## MCMC Simulation with Binary MCMC Model Composition
t = 10000
sim1 = Chains(t, p, names = map(i -> "gamma[$i]", 1:p))
sim2 = Chains(t, p, names = map(i -> "gamma[$i]", 1:p))
gamma1 = BMC3Variate(zeros(p), logf)
gamma2 = BMC3Variate(zeros(p), logf, k=Vector{Int}[[i] for i in 1:p])
for i in 1:t
  sample!(gamma1)
  sample!(gamma2)
  sim1[i, :, 1] = gamma1
  sim2[i, :, 1] = gamma2
end
describe(sim1)
describe(sim2)

#mark that we got to the end of the test file succesfully
@test true


text:
using StatsBase, Distances, Printf

export DTW, FixedPenaltyDTW, FixPenDTW

struct DTW <: SemiMetric
    ground_dist::Metric
end

function (d::DTW)(
    S1::Vector{T}, S2::Vector{T}
    ) where {T}
    if length(S1) < length(S2)  # This ensures first seq is longest
        d(S2, S1)
    else
        @assert (length(S1)>0) & (length(S2)>0) "both args must be either of type Nothing or of nonzero length."
        d_g = d.ground_dist
        prev_row = pushfirst!(fill(Inf, length(S2)), 0.0);
        curr_row = fill(Inf, length(S2) + 1);

        for i = 1:length(S1)
            # curr_row[1] = prev_row[1]
            for j = 1:(length(S2))
                # @show i, j, prev_row[j], d.ground_dist(S1[i], S2[j])
                cost = d_g(S1[i],S2[j])
                curr_row[j+1] = cost + min(prev_row[j], prev_row[j+1], curr_row[j])
            end
            # @show curr_row
            copy!(prev_row, curr_row)
        end
        return curr_row[end]
    end
end

function (d::DTW)(
    S1::Nothing, S2::Vector{T}
    ) where {T}

    d_g = d.ground_dist
    z = 0.0
    for x in S2
        z += d_g(x, nothing)
    end
    return z + length(S2)

end
(d::DTW)(S1::Vector{T}, S2::Nothing) where {T} = d(S2,S1)
(d::DTW)(S1::Nothing, S2::Nothing) where {T} = 0.0

function print_info(
    d::DTW,
    S1::Vector{T}, S2::Vector{T}
    ) where {T}

    d = d.ground_dist
    # First find the substitution matrix
    C = fill(Inf, length(S1)+1, length(S2)+1)
    C[1,1] = 0.0

    for j in 1:length(S2)
        for i in 1:length(S1)
            cost = d(S1[i], S2[j])
            C[i+1,j+1] = cost + min(C[i+1,j], C[i,j+1], C[i,j])
        end
    end
    # Now retrace steps to determine an optimal matching
    i, j = size(C)
    pairs = Tuple{Int,Int}[]
    pushfirst!(pairs, (i-1,j-1))
    while (i ≠ 2) | (j ≠ 2)
        i_tmp, j_tmp = argmin(view(C, (i-1):i, (j-1):j)).I
        i = i - 2 + i_tmp
        j = j - 2 + j_tmp
        pushfirst!(pairs, (i-1,j-1))
    end
    # @show outputs
    title = "Optimal Coupling"
    println(title)
    println("-"^length(title), "\n")
    # for statement in outputs
    #     println(statement)
    # end
    i_tmp,j_tmp = (0,0)
    for (i,j) in pairs
        if i == i_tmp
            println(" "^length(@sprintf("%s",S1[i])) * " ↘ $(S2[j])")
        elseif j == j_tmp
            println("$(S1[i]) ↗ " * " "^length(@sprintf("%s",S2[j])))
        else
            println("$(S1[i]) → $(S2[j])")
        end
        i_tmp, j_tmp = (i,j)
    end

end

# Penalised
# ---------

struct FixedPenaltyDTW <: Metric
    ground_dist::Metric
    ρ::Real
end

const FixPenDTW = FixedPenaltyDTW

function (d::FixPenDTW)(
    S1::Vector{T}, S2::Vector{T}
    ) where {T}
    if length(S1) < length(S2)  # This ensures first seq is longest
        d(S2, S1)
    else
        d_g = d.ground_dist
        prev_row = pushfirst!(fill(Inf, length(S2)), 0.0);
        curr_row = fill(Inf, length(S2) + 1);

        for i = 1:length(S1)
            # curr_row[1] = prev_row[1]
            for j = 1:(length(S2))
                # @show i, j, prev_row[j], d.ground_dist(S1[i], S2[j])
                cost = d_g(S1[i],S2[j])
                curr_row[j+1] = cost + min(
                    prev_row[j], # New pair (no warping)
                    prev_row[j+1] + d.ρ, # Warping
                    curr_row[j] + d.ρ)  # Warping
            end
            # @show curr_row
            copy!(prev_row, curr_row)
        end
        return curr_row[end]
    end
end

function (d::FixPenDTW)(
    S1::Nothing, S2::Vector{T}
    ) where {T}

    d_g = d.ground_dist
    z = 0.0
    for x in S2
        z += d_g(x, nothing)
    end
    return z + length(S2) * d.ρ

end
(d::FixPenDTW)(S1::Vector{T}, S2::Nothing) where {T} = d(S2,S1)
(d::FixPenDTW)(S1::Nothing, S2::Nothing) where {T} = 0.0#

function print_info(
    d::Union{FixPenDTW},
    S1::Vector{T}, S2::Vector{T}
    ) where {T}

    d_g = d.ground_dist
    # First find the substitution matrix
    C = fill(Inf, length(S1)+1, length(S2)+1)
    C[1,1] = 0.0

    for j in 1:length(S2)
        for i in 1:length(S1)
            cost = d_g(S1[i], S2[j])
            C[i+1,j+1] = cost + min(C[i+1,j] + d.ρ, C[i,j+1] + d.ρ, C[i,j])
        end
    end
    # Now retrace steps to determine an optimal matching
    i, j = size(C)
    pairs = Tuple{Int,Int}[]
    warp_shift_mat = [0 d.ρ; d.ρ 0] # Extra bit here different from DTW
    pushfirst!(pairs, (i-1,j-1))
    while (i ≠ 2) | (j ≠ 2)
        i_tmp, j_tmp = argmin(view(C, (i-1):i, (j-1):j) + warp_shift_mat).I
        i = i - 2 + i_tmp
        j = j - 2 + j_tmp
        pushfirst!(pairs, (i-1,j-1))
    end
    # @show outputs
    title = "Fixed Penalty DTW Print-out with $d_g Ground Distance and ρ=$(d.ρ)"
    println(title)
    println("-"^length(title), "\n")
    println("The cheapest way to do the tranformation...\n")
    println(S1, "---->", S2)
    println("\n...is the following series of edits...\n")
    # for statement in outputs
    #     println(statement)
    # end
    i_tmp,j_tmp = (0,0)
    for (i,j) in pairs
        if i == i_tmp
            println(" "^length(@sprintf("%s",S1[i])) * " ↘ $(S2[j])")
        elseif j == j_tmp
            println("$(S1[i]) ↗ " * " "^length(@sprintf("%s",S2[j])))
        else
            println("$(S1[i]) → $(S2[j])")
        end
        i_tmp, j_tmp = (i,j)
    end
end


text:
using DataAndPlotting
using Test

arr = initialize()

#rm("../temp/array.h5")
#output(arr,"../temp/array.h5")

#arr2 = input("../temp/array.h5")

#@test  isequal(arr2, arr)


text:
## Exercise 9-6
# Write a function called isabecedarian that returns true if the letters in a word appear in alphabetical order (double letters are ok). How many abecedarian words are there?

println("Ans: ")

function isabecedarian(word)
    let
        word = uppercase(word)
        prev = word[1]

        for letter in word
            if letter < prev
                return false
            end
            prev = letter
        end
    end

    return true
end

# println(isabecedarian("abcecd"))
count = 0
for word in eachline("words.txt")
    if isabecedarian(word)
        global count += 1
    end
end

println("There are $count abecedarian words.")

println("End.")


text:
import Base: haskey, getindex, setindex!, delete!

export RedBlackTree

###
#
# Red Black Tree
#
# Supported functionality
#   - haskey
#
###

# TODO Can this be re-used somewhere?
# TODO Put this in tree.jl or try to get this out of DataStructures?
abstract Tree{K,V}

type RedBlackTreeNode{K,V} <: Tree{K,V}
    key::K
    value::V
    left::Union(Nothing, Tree{K,V})
    right::Union(Nothing, Tree{K,V})
end

type RedBlackTree{K,V} <: Tree{K,V}
    root::Union(Nothing, RedBlackTreeNode{K,V})
    red::Bool

    #RedBlackTree() = new(EmptyTree{K,V}())
    # Well, this is kind of weird.. is there a way to do this so that we don't have to this to avoid instantiating instances of EmptyTree?
    # RedBlackTree(key::V, value::K) = new(RedBlackTreeNode{K,V}(key, value, EmptyTree{K,V}(), EmptyTree{K,V}()))
    # Is this right?
    RedBlackTree(key::V, value::K) = new(RedBlackTreeNode{K,V}(key, value, nothing, nothing))
end

_isred(tree::Nothing) = false
_isred(tree::RedBlackTree) = tree.red

haskey(tree::Nothing, key) = false
function haskey(tree::RedBlackTree, key)
    rootkey = tree.root
    if key == rootkey
        return true
    elseif key < rootkey
        return haskey(tree.left, key)
    else
        return haskey(tree.right, key)
    end
end

export RedBlackTree


text:
export GameTree, GameNode, GameEdge
export play, move, reward, isfinal

mutable struct GameAction
    id::Int
    P::Float32
    stateid::Int
    N::Int
    Q::Float32
end

GameAction(id::Int, P::Float32) = GameAction(id, P, 0, 0, 0.0)

function backup!(a::GameAction, value::Float32)
    a.Q = (a.N * a.Q + value) / (a.N + 1)
    a.N += 1
end

mutable struct GameState
    data
    p
    v
    reward
    actions::Vector{GameAction}
end

GameState(data, p, v) = GameState(data, p, v, nothing, GameAction[])

isfinal(state::GameState) = isempty(state.actions)

function getpi(state::GameState)
    N = sum(a -> a.N, state.actions)
    @assert N > 0
    pi = zeros(state.p.data)
    for a in state.actions
        pi[a.id] = a.N / N
    end
    pi
end

mutable struct GameTree
    dict::Dict
    states::Vector{GameState}
    model
end

function GameTree(initdata::T, model) where T
    tree = GameTree(Dict{T,Int}(), GameState[], model)
    register!(tree, initdata)
    tree
end

function train(initdata, model)
    opt = SGD()
    for iter = 1:100
        println("iter: $iter")
        nplays = 100
        prog = Progress(nplays)
        states = GameState[]
        for i = 1:nplays
            path = play(initdata, model)
            append!(states, path)
            ProgressMeter.next!(prog)
        end
        opt.rate = 0.1 / length(states)
        continue

        p = concat(2, map(s -> s.p, states))
        pi = cat(2, map(s -> getpi(s), states)...)
        v = concat(2, map(s -> s.v, states))
        r = map(s -> Float32(s.reward), states)
        r = reshape(r, 1, length(r))
        l = softmax_crossentropy(Var(pi), p)
        l += mse(v, Var(r))

        loss = sum(l.data) / length(states)
        println("Loss:\t$loss")

        params = gradient!(l)
        foreach(opt, params)
    end

    println("play start")
    states = play(initdata, model)
    for s in states
        println(string(s.data))
        println(s.p.data)
        println(s.v.data)
        println()
    end
end

function play(initdata, model)
    tree = GameTree(initdata, model)
    states = GameState[]
    state = tree.states[1]
    while !isfinal(state)
        push!(states, state)
        for i = 1:100
            search!(tree, state)
        end
        state = move(tree, state)
    end

    r = -state.reward
    #@assert r == 1 || r == 0
    for i = length(states):-1:1
        states[i].reward = r
        r = -r
    end
    states
end

function move(tree::GameTree, state::GameState)
    maxa = state.actions[1]
    for a in state.actions
        if a.N > maxa.N
            maxa = a
        end
    end
    tree.states[maxa.stateid]
end

function select(state::GameState)
    maxa = state.actions[1]
    length(state.actions) == 1 && return maxa
    maxU = typemin(Float64)
    N = sum(a -> a.N, state.actions)
    cpuct = 3
    for a in state.actions
        U = a.Q + cpuct * a.P * sqrt(N) / (a.N + 1)
        if U > maxU
            maxU = U
            maxa = a
        end
    end
    maxa
end

function register!(tree::GameTree, data)
    id = get!(tree.dict, data, length(tree.dict)+1)
    if id > length(tree.states)
        r = reward(data)
        if r <= -10.0
            actions = getactions(data)
            p, v = tree.model(data)
            state = GameState(data, p, v)
            Ps = softmax(p.data[actions])
            for i = 1:length(actions)
                P = Ps[i]
                push!(state.actions, GameAction(actions[i],P))
            end
        else
            actions = Int[]
            p = nothing
            v = Var(Float32[r])
            state = GameState(data, p, v)
        end
        state.reward = r
        push!(tree.states, state)
        @assert length(tree.states) == length(tree.dict)
    end
    id
end

function search!(tree::GameTree, initstate::GameState)
    state = initstate
    actions = GameAction[]
    while !isfinal(state)
        action = select(state)
        push!(actions, action)
        if action.stateid == 0
            data = move(state.data, action.id)
            action.stateid = register!(tree, data)
            state = tree.states[action.stateid]
            break
        else
            state = tree.states[action.stateid]
        end
    end

    v = state.v.data[1]
    for i = length(actions):-1:1
        backup!(actions[i], v)
        v = -v
    end
end

export TicTacToe

struct TicTacToeGame
end

const WHITE = Int8(1)
const BLACK = Int8(2)
const BLANK = Int8(3)

mutable struct TicTacToe
    player::Int
    board::Vector{Int8}
    count::Int
end
TicTacToe() = TicTacToe(WHITE, fill(BLANK,9), 0)

Base.isequal(x::TicTacToe, y::TicTacToe) = isequal(x.board, y.board)
Base.hash(x::TicTacToe) = hash(x.board)

function Base.string(ttt::TicTacToe)
    s = map(ttt.board) do i
        i == WHITE && return "○"
        i == BLACK && return "×"
        i == BLANK && return "_"
        throw("")
    end
    p = ttt.player == WHITE ? "○" : "×"
    "player: $p\n$(s[1]) $(s[2]) $(s[3])\n$(s[4]) $(s[5]) $(s[6])\n$(s[7]) $(s[8]) $(s[9])"
end

function getactions(ttt::TicTacToe)
    actions = Int[]
    for i = 1:length(ttt.board)
        ttt.board[i] == BLANK && push!(actions,i)
    end
    actions
end

function reward(ttt::TicTacToe)
    pats = [(1,2,3), (4,5,6), (7,8,9), (1,4,7), (2,5,8), (3,6,9), (1,5,9), (3,5,7)]
    for p in pats
        ttt.board[p[1]] == BLANK && continue
        if ttt.board[p[1]] == ttt.board[p[2]] == ttt.board[p[3]]
            @assert ttt.board[p[1]] != ttt.player
            return -1.0
        end
    end
    ttt.count < 9 && return -100.0
    0.0
end

function move(ttt::TicTacToe, act::Int)
    @assert ttt.board[act] == BLANK
    board = copy(ttt.board)
    board[act] = ttt.player
    player = ttt.player == WHITE ? BLACK : WHITE
    TicTacToe(player, board, ttt.count+1)
end

mutable struct NN
    g
end

function NN()
    T = Float32
    embeds = zerograd(Uniform(-0.001,0.001)(T,10,3))
    h = lookup(Node(embeds), Node(name="x"))
    h = Linear(T,90,90)(h)
    h = tanh(h)
    p = Linear(T,90,9)(h)
    v = Linear(T,90,1)(h)
    v = tanh(v)
    NN(Graph(p,v))
end

function (nn::NN)(ttt::TicTacToe)
    x = ttt.board
    p = ttt.player
    if p == BLACK
        x = copy(x)
        for i = 1:length(x)
            x[i] == BLANK && continue
            x[i] = x[i] == WHITE ? BLACK : WHITE
        end
    end
    nn.g(Var(x))
end


text:
using GraphPCA
using Base.Test

# X matrix
X1 = 10*randn(20,4)
X2 = X1 + randn(20,4)
X = [X1 X2]

#Construct the Laplacian matrix
L = eye(8)
for i in 1:4
  j = i + 4
  L[i,j] = -1
  L[j,i] = -1
end

λ = eigmax(Symmetric(X'*X))
η = eigmax(Symmetric(L))
eet = ones(8,8)
n = 8

#Pure PCA
XtX = Symmetric(I - X'*X/λ)
Dpca,Vpca = eig(XtX, 1:4)
Qpca = Vpca
Upca = X*Qpca
Xpca = Upca*Qpca'

#Pure Laplacian Embedding
Ln = Symmetric(L/η)
Dn, Vn = eig(Ln, 1:4)
Qlap = Vn
Ulap = X*Qlap
Xlap = Ulap*Qlap'

#Laplacian Embedding with Offset
# Lnp = Symmetric(L/η + eet/n)
# Dnp,Vnp = eig(Lnp, 1:4)
# Qlapp = Vnp
# Ulapp = X*Qlapp
# Xlapp = Ulapp*Qlapp'

GLPCA0 = GLPCA(X, 0, L, 4)
GLPCA1 = GLPCA(X, 1, L, 4)
@test GLPCA1.U*GLPCA1.Q' ≈ Xlap
@test GLPCA0.U*GLPCA0.Q' ≈ Xpca

max_stars_repo_name: mihirparadkar/GraphPCA

max_issues_repo_name: mihirparadkar/GraphPCA

max_forks_repo_name: mihirparadkar/GraphPCA


text:
function sinkhornTransport(a,b,K,U,lambda,stoppingCriterion="marginalDifference",p_norm=Inf,tolerance=.5e-2,maxIter= 5000,VERBOSE=0)
    # Compute the N Sinkhorn divergence for all the pairs
    # D = [d(a_1,b_1),...,d(a_N,b_N)]
    #--------------------------------------------------
    # Inputs:
    #--------------------------------------------------
    # a : is either
    #   - a d1 x 1 column vector in the probability simpolex (non-negative summing to 1). This is the [1-vs-N mode]
    #   - a d1 x N matrix, where each column vector is in the probability simplex.This is the [N x 1-vs-1 mode]
    # b : is a d2 x N matrix of N vectors in the probability simplex
    #
    # K is a d1 x d2 matrix, equal to exp(-lambda M), where M is the d1 x d2 matrix of pairwise distances between
    # bins described in a and bins in the b_1,...b_N histograms.
    # In the most simple case d_1=d_2 and M is simply a distance matrix (zero on the diagonal and such that m_ij < m_ik + m_kj
    #
    # U = K.*M is a d1 x d2 matrix, pre-stored to speed up the computation of the distances.
    #--------------------------------------------------
    # Optional Inputs:
    #--------------------------------------------------
    # Outputs:
    #--------------------------------------------------
    # D : vector of N dual-sinkhorn divergences, or upper bounds to the EMD.
    #
    # L : vector of N lower bounds to the original OT problem, a.k.a EMD. This is computed by using the dual variables of
    # the smoothed problem, which, when modified adequately, are feasible for the original (non-smoothed) OT dual problem
    #
    # u : d1 x N matrix of left scalings
    # v : d2 x N matrix of right scalings
    # The smoothed optimal transport between (a_i,b_i) can be recovered as
    # T_i = diag(u(:,i)) * K * diag(v(:,i));
    #
    # or, equivalently and substantially faster:
    # TO-DO: the julia syntax is not working now
    # T_i = broadcast(*,v(:,i)',(broadcast(*,u(:,i),K)))
    # Check this for details: https://docs.julialang.org/en/v0.6.1/manual/arrays/
    #
    # Relevant paper:
    # M. Cuturi,
    # Sinkhorn Distances : Lightspeed Computation of Optimal Transport,
    # Advances in Neural Information Processing Systems (NIPS) 26, 2013
    #
    # Optional Input Check:
    #--------------------------------------------------
    # First make sure that a and b are two dimensional array
    # TO-DO:
    #
    ## Checking the type of computation: 1-vs-N points or many pairs
    a_r1, a_c = size(a)
    b_r, b_c = size(b)

    if a_c == 1
        ONE_VS_N = true # We are computing [D(a,b_1), ... , D(a,b_N)]
    elseif a_c == b_c
        ONE_VS_N = false # We are computing [D(a_1,b_1), ... , D(a_N,b_N)]
    else
        error("The first parameter a is either a column vector in the probability simplex, or N column vectors in the probability simplex where N is size(b,2)")
    end
    ## Checking dimensionality:
    if b_c > b_r
        BIGN = true
    else
        BIGN = false
    end
    ## Small changes in the 1-vs-N case to go a bit faster.
    if ONE_VS_N # if computing 1-vs-N make sure all components of a are >0. Otherwise we can get rid of some   lines of K to go faster.
        I = find(a) # Changed from Matlab
        someZeroValues = false
        if length(I)<a_r1 # need to update some vectors and matrices if a does not have full support
            someZeroValues=true
            K = K[I,:]
            U = U[I,:]
            a = a[I]
        end
        ainvK= broadcast(/,K,a) # precomputation of this matrix saves a d1 x N Schur product at each iteration.
    end
    # display(ainvK[1:5,1:5])

    a_r = length(a)
    ## Fixed point counter
    compt=0;
    ## Initialization of
    u =  (1/a_r).*ones(a_r,b_c)./a_r #broadcast(/,ones(a_r,b_c),a_r)#
    v = zeros(b_r,b_c)
    if cmp(stoppingCriterion, "distanceRelativeDecrease")==0
        Dold = ones(1,b_c); # initialization of vector of distances.
    end
    while compt < maxIter
        if ONE_VS_N # 1-vs-N mode
            if BIGN
                u = 1./(ainvK*broadcast(/,b,(K'*u)))
            else
                u=1./(ainvK*broadcast(/,b,(u'*K)'))
            end
        else
            if BIGN
                u=broadcast(/,a,(K*broadcast(/,b,(u'*K)')))
            else
                u=broadcast(/,a,(K*broadcast(/,b,(K'*u))))
            end
        end
        compt=compt+1
        # check the stopping criterion every 20 fixed point iterations
        # or, if that's the case, before the final iteration to store the most
        # recent value for the matrix of right scaling factors v.
        if mod(compt,20)==1 || compt==maxIter
        # split computations to recover right and left scalings.
            if BIGN
                v=broadcast(/,b,(K'*u)) # main iteration of Sinkhorn's algorithm
            else
                v=broadcast(/,b,(u'*K)')
            end
            if ONE_VS_N # 1-vs-N mode
                u=1./(ainvK*v)
            else
                u=broadcast(/,a,(K*v))
            end
            # check stopping criterion
            if cmp(stoppingCriterion, "distanceRelativeDecrease")==0
                D=sum(broadcast(*,u,(U*v)))
                Criterion=norm(broadcast(/,D,Dold)-1,p_norm)
                if Criterion<tolerance || isnan(Criterion)
                    break
                end
                 Dold=D
            elseif cmp(stoppingCriterion,"marginalDifference")==0
                Criterion=norm(sum(abs.(broadcast(*,v,(K'*u))-b)),p_norm)
                if Criterion<tolerance || isnan(Criterion)
                    break
                end
            else
                error("Stopping Criterion not recognized")
            end# norm of all || . ||_1 differences between the marginal of the current solution with the actual marginals.
            compt=compt+1
            if VERBOSE>0
                display("Iteration : $compt, Criterion: $Criterion")
            end
            # if any(isnan(Criterion)), # stop all computation if a computation of one of the pairs goes wrong.
            #     error("NaN values have appeared during the fixed point iteration. This problem appears because of insufficient machine precision when processing computations with a regularization value of lambda that is too high. Try again with a reduced regularization parameter lambda or with a thresholded metric matrix M.");
            # end
        end
    end
    if cmp(stoppingCriterion,"marginalDifference")==0 # if we have been watching marginal differences, we need to compute the vector of distances.
        D=sum(u.*(U*v))
    end

    alpha = log.(u)
    beta = log.(v)
    [beta[t]==-Inf ? beta[t]=0 : beta[t]=beta[t] for t in eachindex(beta)]
    if ONE_VS_N
        L= (a'* alpha + sum(b.*beta))/lambda
    else
        [alpha[t]==-Inf ? alpha[t]=0 : alpha[t]=alpha[t] for t in eachindex(alpha)]
        L= (sum(broadcast(*,a,alpha)) + sum(broadcast(*,b,beta)))/lambda
    end
    # display(L)
    if ONE_VS_N && someZeroValues
        uu=u
        u=zeros(a_r1,b_c)
        u[I,:]=uu
    end
    return D, L, u, v
end


text:
include("../oracles/shortest_path_oracle.jl")
include("../oracles/portfolio_oracle.jl")
include("../solver/util.jl")
include("../solver/sgd.jl")
include("../solver/reformulation.jl")
include("../solver/random_forests_po.jl")
include("../solver/validation_set.jl")
include("../experiments/replication_functions.jl")

using Random, Distributions, Test, CSV

function test_sp_replication()
    Random.seed!(8392)

    envOracle = setup_gurobi_env(method_type = :default, use_time_limit = true)
    envReform = setup_gurobi_env(method_type = :method3, use_time_limit = true)

    final_results = shortest_path_replication(5,
        100, 25, 10000,
        5, 8, 0.5;
        num_lambda = 10, lambda_max = 1.0, lambda_min_ratio = 10.0^(-8), regularization = :lasso,
        gurobiEnvOracle = envOracle, gurobiEnvReform = envReform,
        different_validation_losses = false)

    println(final_results)
end

function test_port_replication()
    Random.seed!(780)

    envOracle = setup_gurobi_env(method_type = :default, use_time_limit = true)
    envReform = setup_gurobi_env(method_type = :method3, use_time_limit = true)

    final_results = portfolio_replication(200, 4, 1000, 250, 10000, 5,
        8, 1;
        num_lambda = 10, lambda_max = 100.0, lambda_min_ratio = 10.0^(-8),
        gurobiEnvOracle = envOracle, gurobiEnvReform = envReform, different_validation_losses = false,
        data_type = :poly_kernel)

    println(final_results)
end

max_stars_repo_name: paulgrigas/SmartPredictThenOptimize

max_issues_repo_name: paulgrigas/SmartPredictThenOptimize

max_forks_repo_name: paulgrigas/SmartPredictThenOptimize


text:
using Flux: onehotbatch
using WAV
using BSON

# This is a custom fork of MFCC for the time being, since the project hasn't
# accepted my pull request that updated the depreacted `iceil` to `ceil`
# https://github.com/maetshju/MFCC.jl
using MFCC

# Define constants that will be used
const TRAINING_DATA_DIR = "TIMIT/TRAIN"
const TEST_DATA_DIR = "TIMIT/TEST"

const TRAINING_OUT_DIR = "train"
const TEST_OUT_DIR = "test"

# Make dictionary to map from phones to class numbers
const PHONES = split("h#	q	eh	dx	iy	r	ey	ix	tcl	sh	ow	z	s	hh	aw	m	t	er	l	w	aa	hv	ae	dcl	y	axr	d	kcl	k	ux	ng	gcl	g	ao	epi	ih	p	ay	v	n	f	jh	ax	en	oy	dh	pcl	ah	bcl	el	zh	uw	pau	b	uh	th	ax-h	em	ch	nx	eng")
translations = Dict(phone=>i for (i, phone) in enumerate(PHONES))
translations["sil"] = translations["h#"]
const PHONE_TRANSLATIONS = translations

# The frame length and interval are guesses because the paper does not
# specify them. Most literature seems to use these values, however.
FRAME_LENGTH = 0.025 # ms
FRAME_INTERVAL = 0.010 # ms

"""
    makeFeatures(wavFname, phnFname)

Extracts Mel filterbanks and associated labels from `wavFname` and `phnFaname`.
"""
function makeFeatures(phnFname, wavFname)
    samps, sr = wavread(wavFname)
    samps = vec(samps)

    feats = mfcc(samps, sr, minfreq=64.0, maxfreq=8000.0)[1]
    frames = powspec(samps, sr; wintime=FRAME_LENGTH, steptime=FRAME_INTERVAL)

    # I am assuming we want log energy because Mel scale is a log scale
    # If not, we can remove the log

    # Transpose so that we get each row as an observation of variables, as
    # opposed to each column as an observation of variables; not that one is
    # inherently better than the other, but most work I've encountered seems to
    # prefer rows as observations, and that's how I think about the problems too
    energies = log.(sum(frames', dims=2))

    feats[:,1] = energies[:,1]

    local lines
    open(phnFname, "r") do f
        lines = readlines(f)
    end

    boundaries = Vector()
    labels = Vector()

    # first field in the file is the beginning sample number, which isn't
    # needed for calculating where the labels are
    for line in lines
        _, boundary, label = split(line)
        boundary = parse(Int64, boundary)
        push!(boundaries, boundary)
        push!(labels, label)
    end

    labelInfo = collect(zip(boundaries, labels))
    labelInfoIdx = 1
    boundary, label = labelInfo[labelInfoIdx]
    nSegments = length(labelInfo)

    frameLengthSamples = FRAME_LENGTH * sr
    frameIntervalSamples = FRAME_INTERVAL * sr
    halfFrameLength = FRAME_LENGTH / 2

    # Begin generating sequence labels by looping through the acoustic
    # sample numbers

    labelSequence = Vector() # Holds the sequence of labels

    idxsToDelete = Vector() # To store indices for frames labeled as 'q'
    for i=1:size(feats, 1)
        win_end = frameLengthSamples + (i-1)*frameIntervalSamples

        # Move on to next label if current frame of samples is more than half
        # way into next labeled section and there are still more labels to
        # iterate through
        if labelInfoIdx < nSegments && win_end - boundary > halfFrameLength

            labelInfoIdx += 1
            boundary, label = labelInfo[labelInfoIdx]
        end

        if label == "q"
            push!(idxsToDelete, i)
            continue
        end

        push!(labelSequence, label)
    end

    # Remove the frames that were labeld as 'q'
    feats = feats[[i for i in 1:size(feats,1) if !(i in Set(idxsToDelete
    ))],:]

    feats_deltas = deltas(feats)
    #feats_delta_deltas = deltas(feats_deltas)
    features = hcat(feats, feats_deltas)
    #features = hcat(feats, feats_deltas, feats_delta_deltas)
    return (features, labelSequence)
end

"""
    createData(data_dir, out_dir)

Extracts data from files in `data_dir` and saves results in `out_dir`.
"""
function createData(data_dir, out_dir)
    for (root, dirs, files) in walkdir(data_dir)

        # Exclude the files that are part of the speaker accent readings
        files = [x for x in files if ! occursin("SA", x)]

        phnFnames = [x for x in files if occursin("PHN", x) && ! occursin("BUCK", x)]
        wavFnames = [x for x in files if occursin("WAV", x)]

        one_dir_up = basename(root)
        println(root)

        for (wavFname, phnFname) in zip(wavFnames, phnFnames)
            phn_path = joinpath(root, phnFname)
            wav_path = joinpath(root, wavFname)

            x, y = makeFeatures(phn_path, wav_path)

	    # Convert labels to class numbers
            y = [PHONE_TRANSLATIONS[x] for x in y]

            # Generate class nums; there are 61 total classe. However,
	    # because we're using connectionist temporal classification
	    # loss, we need an extra class, so we go up to 62.
            class_nums = [n for n in 1:62]
            y = onehotbatch(y, class_nums)'

            base, _ = splitext(phnFname)
            dat_name = one_dir_up * base * ".bson"
            dat_path = joinpath(out_dir, dat_name)
            BSON.@save dat_path x y
            # save(dat_path, "x", x, "y", y)
        end
    end
end

createData(TRAINING_DATA_DIR, TRAINING_OUT_DIR)
createData(TEST_DATA_DIR, TEST_OUT_DIR)

max_stars_repo_name: maetshju/flux-ctc-grad

max_issues_repo_name: maetshju/flux-ctc-grad

max_forks_repo_name: maetshju/flux-ctc-grad


text:
# Generates bind.jl
# To run execute `julia --project=. bind_generator.jl` in this directory

using Clang.Generators
using HDF5_jll

include_dir = normpath(HDF5_jll.artifact_dir, "include")
#clang_dir = joinpath(include_dir, "clang-c")

# wrapper generator options
options = load_options(joinpath(@__DIR__, "bind_generator.toml"))

# add compiler flags, e.g. "-DXXXXXXXXX"
args = get_default_args()
push!(args, "-I$include_dir -DH5_BUILT_AS_DYNAMIC_LIB")

headers = [normpath(include_dir, "hdf5.h")]
#headers = [joinpath(clang_dir, header) for header in readdir(clang_dir) if endswith(header, ".h")]
# there is also an experimental `detect_headers` function for auto-detecting top-level headers in the directory
# headers = detect_headers(clang_dir, args)

# create context
ctx = create_context(headers, args, options)

# run generator
#build!(ctx)
build!(ctx, BUILDSTAGE_NO_PRINTING)

function rewrite!(e::Expr)
    if e.head == :function
        body = e.args[2]
        for s in body.args
            if s.head == :call && s.args[1] == :ccall
                b = ccall_to_bind(s)
                e.head = b.head
                e.args = b.args
            end
        end
    end
end

function rewrite!(dag::ExprDAG)
    for node in get_nodes(dag)
        for expr in get_exprs(node)
            rewrite!(expr)
        end
    end
end

function ccall_to_bind(e::Expr)
    ## Extraction the information from the ccall
    func = e.args[2].args[1].value
    retval = e.args[3]
    types = e.args[4].args
    args = e.args[5:end]

    ## Reformat the information
    try
        m = match(r"([A-Z0-9]+)(.*)", string(func))
        s = lowercase(m.captures[1]) * "_" * m.captures[2]
        func = Symbol(s)
    catch err
        @error err backtrace()
    end
    arg_type_pairs = [:($a::$t)  for (a,t) in zip(args, types)]

    ## Build the macro expression
    template = :(@bind f()::ret "Error unspecified")
    # Remove line annotation
    template.args[2] = Expr(:nothing)
    template.args[3].args[2] = retval
    template.args[3].args[1].args[1] = func
    append!(template.args[3].args[1].args, arg_type_pairs)
    template.args[4] = "Error in " * string(func) * " (not annotated)"
    template
end

rewrite!(ctx.dag)

build!(ctx, BUILDSTAGE_PRINTING_ONLY)


text:
#########
# Setup #
#########

using CSV
using DataFrames
using SynapticFilter

i = parse(Int, ENV["LSB_JOBINDEX"]) # set this to a fixed number when not on the cluster

num_η = 10
η = 2.0 ^ (-(i % num_η)) # learning rate as a function of the job index

ρ         = 40    # input firing rate
τm        = 0.025 # membrane time constant
ε0        = 1     # EPSP amplitude           <------ !!!
τ         = 200   # OU time constant
σs        = 1     # OU variance
g0        = 20    # baseline output firing rate
β0        = 0.1   # gain for numblocks = 1 when scaling beta with dimension
β         = 0.1   # fixed beta
blocksize = 8     # size of input block
τblock    = 1.0   # duration of block activation

# values for the number of blocks to iterate over
numblock_rng = [1, 2, 3, 4, 6, 8, 11, 16, 23, 32, 45, 64, 91, 128, 181, 256]

epoch      = ceil(Int, τ)                   # duration of one epoch
num_epochs = parse.(Int, ARGS)[1]           # number of epochs to simulate
num_burnin = parse.(Int, ARGS)[2]           # number of epochs to burn in
ts_per_sec = 1000                           # number of timesteps per seconds
                                            # (reciprocal of timestep)
num_timesteps = num_epochs*epoch*ts_per_sec # total number of timesteps
burnin        = num_burnin*epoch*ts_per_sec # number of timesteps for burn in
timestep      = 1/ts_per_sec                # time step

results = DataFrame(
    filter       = String[],
    eta          = Float64[],
    sparse_input = Bool[],
    beta         = Float64[],
    dim          = Int[],
    mse          = Float64[]
)
l = ReentrantLock()

##################
# Saving results #
##################

function saveresults(df)
    #"/cluster/home/ssurace/perf_sims/csv/nonsparse_scaling_vs_sparse_nonscaling"
    CSV.write("csv/nonsparse_scaling_vs_sparse_nonscaling_grad"  * lpad(i, 4, "0") * ".csv", df)
end#saveresults

###########################
# Performance simulations #
###########################

smodel = OUModel(τ, σs)

@time begin # measure total time for running performance sims

# performance sims for GradientRule, loop over number of blocks in parallel
Threads.@threads for numblocks in numblock_rng
    dim = numblocks * blocksize
    β1 = β0 / sqrt(numblocks)

    imodel = PoissonExpModel(ρ, τm, ε0, dim)
    omodel = ExpGainModel(g0, β1)
    filter = GradientRule(dim, smodel, omodel, η)

    sim = Simulation(imodel, smodel, omodel, filter)

    mse = ComputeMSE(
        sim; num_timesteps = num_timesteps, timestep = timestep, burnin = burnin
    )
    println("GradientRule for numblocks = ", numblocks, ", nonsparse input:   ", mse)

    lock(l) do
        push!(
            results,
            Dict(
                :filter       => "GradientRule",
                :eta          => η,
                :sparse_input => false,
                :beta         => β1,
                :dim          => dim,
                :mse          => mse
            )
        )
    end
    saveresults(results)

    imodel = BlockPoissonExpModel(ρ, τm, ε0, numblocks, blocksize, τblock)
    omodel = ExpGainModel(g0, β)

    sim = Simulation(imodel, smodel, omodel, filter)

    mse = ComputeMSE(
        sim; num_timesteps = num_timesteps, timestep = timestep, burnin = burnin
    )
    println("GradientRule for numblocks = ", numblocks, ", sparse input: ", mse)

    lock(l) do
        push!(
            results,
            Dict(
                :filter       => "GradientRule",
                :eta          => η,
                :sparse_input => true,
                :beta         => β,
                :dim          => dim,
                :mse          => mse
            )
        )
    end
    saveresults(results)
end

println("
TOTAL SIMULATION TIME: ")
end# @time


text:
# Copyright 2015 TWO SIGMA OPEN SOURCE, LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#        http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module Beaker

using JSON
using Requests

core_url = string("http://beaker:", ENV["beaker_core_password"], "@127.0.0.1:", ENV["beaker_core_port"])

function setsession(id)
  global sessionid = id
end

function get(var)
  reply = Requests.json(Requests.get(string(core_url, "/rest/namespace/get"); query = Dict("name" => var, "session" => sessionid)))
  if !reply["defined"]
    error(string("name '", var, "' is not defined in notebook namespace"))
  end
  return reply["value"]
end

function set4(var, val, unset, sync)
  args = Dict("name" => var, "session" => sessionid, "sync" => sync)
  if !unset
    args["value"] = JSON.json(val)
  end
  reply = readall(Requests.post(string(core_url, "/rest/namespace/set"); data = args))
  if reply != "ok"
    error(reply)
  end
end

function set(var, val)
  return set4(var, val, false, true)
end

end

max_stars_repo_name: tbenst/beaker-webppl

max_issues_repo_name: tbenst/beaker-webppl

max_forks_repo_name: tbenst/beaker-webppl


text:
include("Exceptions.jl")
include("SerialDevices/SerialDevices.jl")
include("SCPI/SCPIInstruments.jl")
include("DictToStruct.jl")
include("BidirectionalChannel.jl")
include("ScannerCoordinates.jl")

function wait(::Nothing)
 # NOP
end

# I only add this here until https://github.com/JuliaLang/julia/pull/42272 is decided.
Base.convert(::Type{IPAddr}, str::AbstractString) = parse(IPAddr, str)


text:
import Base.print

"""
Terminal
"""
struct Terminal
    buffers::Vector{Buffer}
    current::Ref{UInt8}
    cursor_hidden::Ref{Bool}
    terminal_size::Ref{Rect}
    keyboard_buffer::Vector{Char}
    stdout_channel::Channel{String}
    stdin_channel::Channel{Char}
    kind::String
    wait::Float64
    ispaused::Ref{Bool}
    function Terminal(stdout_io = stdout, stdin_io = stdin)
        width, height = terminal_size()
        rect = Rect(1, 1, width, height)
        buffers = [Buffer(rect), Buffer(rect)]
        current = 1
        cursor_hidden = false
        stdout_channel = Channel{String}(Inf)
        stdin_channel = Channel{Char}(Inf)
        ispaused = Ref{Bool}(false)
        stdout_task = @async while true
            if ispaused[] == true
                sleep(1)
                continue
            end
            print(stdout_io, take!(stdout_channel))
        end
        stdin_task = @async while true
            if ispaused[] == true
                sleep(1)
                continue
            end
            c = Char(read(stdin_io, 1)[])
            put!(stdin_channel, c)
        end
        t = new(buffers, current, cursor_hidden, rect, Char[], stdout_channel, stdin_channel, get(ENV, "TERM", ""), 1 / 1000, ispaused)
        TERMINAL[] = t
        return t
    end
end

"""
Get key press
"""
function get_event(t)
    return if isready(t.stdin_channel)
        take!(t.stdin_channel)
    else
        nothing
    end
end

function update_channel(t, args...)
    for arg in args
        update_channel(t, string(arg))
    end
end
update_channel(t, arg) = update_channel(t, string(arg))
update_channel(t, arg::String) = put!(t.stdout_channel, arg)

const TERMINAL = Ref{Terminal}()

"""
Flush terminal contents
"""
function flush(t::Terminal, diff = true)
    previous_buffer = t.buffers[END - t.current[]]
    current_buffer = t.buffers[t.current[]]
    if diff
        draw(t, previous_buffer, current_buffer)
    else
        draw(t, current_buffer)
    end
    while isready(t.stdout_channel)
        sleep(1e-6)
    end
    update(t)
end

"""
Move cursor
"""
move_cursor(t::Terminal, row, col) = update_channel(t, Terminals.CSI, INDICES[row], ';', INDICES[col], 'H')
"""
Move cursor up
"""
move_cursor_up(t::Terminal, row = 1) = update_channel(t, Terminals.CSI, INDICES[row], 'A')
"""
Move cursor down
"""
move_cursor_down(t::Terminal, row = 1) = update_channel(t, Terminals.CSI, INDICES[row], 'B')
"""
Move cursor right
"""
move_cursor_right(t::Terminal, col = 1) = update_channel(t, Terminals.CSI, INDICES[col], 'C')
"""
Move cursor left
"""
move_cursor_left(t::Terminal, col = 1) = update_channel(t, Terminals.CSI, INDICES[col], 'D')
"""
Move cursor home
"""
move_cursor_home(t::Terminal) = update_channel(t, HOMEPOSITION)

"""
Clear screen
"""
clear_screen(t::Terminal) = update_channel(t, CLEARSCREEN)
"""
Clear screen from cursor up onwards
"""
clear_screen_from_cursor_up(t::Terminal) = update_channel(t, CLEARBELOWCURSOR)
"""
Clear screen from cursor down onwards
"""
clear_screen_from_cursor_down(t::Terminal) = update_channel(t, CLEARABOVECURSOR)

"""
Clear line
"""
clear_line(t::Terminal) = update_channel(t, CLEARLINE)
"""
Clear line from cursor right onwards
"""
clear_line_from_cursor_right(t::Terminal) = update_channel(t, CLEARAFTERCURSOR)
"""
Clear line from cursor left onwards
"""
clear_line_from_cursor_left(t::Terminal) = update_channel(t, CLEARBEFORECURSOR)

"""
Hide cursor
"""
hide_cursor(t::Terminal) = update_channel(t, HIDECURSOR)
"""
Show cursor
"""
show_cursor(t::Terminal) = update_channel(t, SHOWCURSOR)

"""
Save cursor
"""
save_cursor(t::Terminal) = update_channel(t, SAVECURSOR)
"""
Restore cursor
"""
restore_cursor(t::Terminal) = update_channel(t, RESTORECURSOR)

"""
Change cursor to blinking block
"""
change_cursor_to_blinking_block(t::Terminal) = update_channel(t, CURSORBLINKINGBLOCK)
"""
Change cursor to steady block
"""
change_cursor_to_steady_block(t::Terminal) = update_channel(t, CURSORSTEADYBLOCK)
"""
Change cursor to blinking underline
"""
change_cursor_to_blinking_underline(t::Terminal) = update_channel(t, CURSORBLINKINGUNDERLINE)
"""
Change cursor to steady underline
"""
change_cursor_to_steady_underline(t::Terminal) = update_channel(t, CURSORSTEADYUNDERLINE)
"""
Change cursor to blinking ibeam
"""
change_cursor_to_blinking_ibeam(t::Terminal) = update_channel(t, CURSORBLINKINGIBEAM)
"""
Change cursor to steady ibeam
"""
change_cursor_to_steady_ibeam(t::Terminal) = update_channel(t, CURSORSTEADYIBEAM)

"""
Reset terminal settings
"""
reset(t::Terminal) = update_channel(t, RESET)

"""
Alternate screen TUI mode
"""
tui_mode(t::Terminal) = update_channel(t, TUIMODE)
"""
Default mode
"""
default_mode(t::Terminal) = update_channel(t, DEFAULTMODE)

"""
Get current buffer
"""
current_buffer(t::Terminal)::Buffer = t.buffers[t.current[]]
current_buffer()::Buffer = current_buffer(TERMINAL[])

"""
Reset terminal
"""
reset(t::Terminal, buffer::Int) = reset(t.buffers[buffer])

"""
Update terminal
"""
function update(t::Terminal)
    reset(t.buffers[END - t.current[]])
    t.current[] = END - t.current[]
    w, h = terminal_size()
    if t.terminal_size[].width != w || t.terminal_size[].height != h
        resize(t, w, h)
        clear_screen(t)
        move_cursor_home(t)
        # TODO: we need to redraw here.
    end
end

"""
Draw terminal
"""
function draw(t::Terminal, buffer::Buffer)
    save_cursor(t)
    move_cursor_home(t)
    iob = IOBuffer()
    for cell in permutedims(buffer.content)[:]
        print(iob, cell.style, cell.char, inv(cell.style))
    end
    update_channel(t, String(take!(iob)))
    restore_cursor(t)
    while isready(t.stdout_channel)
        sleep(1e-6)
    end
    update(t)
end

"""
Draw terminal
"""
function draw(t::Terminal, buffer1::Buffer, buffer2::Buffer)
    save_cursor(t)
    b1 = buffer1.content[:]
    b2 = buffer2.content[:]
    move_cursor_home(t)
    R, C = size(buffer2.content)
    for r = 1:R, c = 1:C
        if buffer1.content[r, c] != buffer2.content[r, c]
            move_cursor(t, r, c)
            cell = buffer2.content[r, c]
            update_channel(t, cell.style, cell.char, inv(cell.style))
        end
    end
    restore_cursor(t)
end

"""
Resize terminal
"""
function resize(t::Terminal, w::Int, h::Int)
    rect = Rect(1, 1, w, h)
    t.terminal_size[] = rect
    t.buffers[1] = Buffer(rect)
    t.buffers[2] = Buffer(rect)
end

"""
Locate cursor
"""
function locate_cursor(t::Terminal)
    ucs = Char[]
    while length(t.keyboard_buffer) > 0
        push!(ucs, popfirst!(t.keyboard_buffer))
    end
    location = UInt8[]
    with_raw_mode() do
        print(stdout, LOCATECURSOR)
        Base.flush(stdout)
        total_read = 0
        now = time()
        c = 0x00
        bell = UInt8(BELL)
        channel = Channel(1)
        t = @async begin
            while c != UInt8('R')
                c = read(stdin, 1)[]
                push!(location, c)
                total_read += 1
            end
        end
        while c != UInt8('R') && time() - now < CONTROL_SEQUENCE_TIMEOUT
            sleep(1e-3)
        end
        if !istaskdone(t)
            t.exception = InterruptException()
        end
    end
    row, col = split(String(Char.(location)[3:end-1]), ";")
    return parse(Int, row), parse(Int, col)
end


text:
using DataFrames, MLDataUtils, KernelFunctions, MVDA, StableRNGs
using LinearAlgebra, Statistics
using Polyester

# using MKL
BLAS.set_num_threads(8)

run = function(dir, example, algorithm, data, sparse2dense::Bool=false; at::Real=0.8, nfolds::Int=5, kwargs...)
    # Create MVDAProblem instance w/ RBFKernel and construct grids.
    @info "Creating MVDAProblem for $(example)"
    dist = Float64[]
    for j in axes(data[2], 1)
        @views yⱼ, xⱼ = data[1][j], data[2][j, :]
        @batch per=core threadlocal=Float64[]::Vector{Float64} for i in axes(data[2], 1)
            @views yᵢ, xᵢ = data[1][i], data[2][i, :]
            if yᵢ != yⱼ
                push!(threadlocal, norm(xᵢ-xⱼ))
            end
        end
        foreach(subset -> append!(dist, subset), threadlocal)
    end
    σ = median(dist)
    problem = MVDAProblem(data..., intercept=true, kernel=σ*RBFKernel())
    n, p, c = MVDA.probdims(problem)
    n_train = round(Int, n * at * (nfolds-1)/nfolds)
    n_validate = round(Int, n * at * 1/nfolds)
    n_test = round(Int, n * (1-at))
    fill!(problem.coeff.all, 1/(n_train+1))
    fill!(problem.coeff_prev.all, 1/(n_train+1))
    ϵ_grid = [MVDA.maximal_deadzone(problem)]

    if example == "HAR"
        s_grid = sort!(collect(range(0.0, 1.0, length=21)), rev=sparse2dense)
    else
        s_grid = sort!([1-k/n_train for k in n_train:-1:0], rev=sparse2dense)
    end
    grids = (ϵ_grid, s_grid)

    # Collect data for cross-validation replicates.
    @info "CV split: $(n_train) Train / $(n_validate) Validate / $(n_test) Test"
    subsets = (n_train, n_validate, n_test)
    rng = StableRNG(1903)
    replicates = MVDA.cv_estimation(algorithm, problem, grids;
        at=at,                  # propagate CV / Test split
        nfolds=nfolds,          # propagate number of folds
        rng=rng,                # random number generator for reproducibility
        nouter=10^2,            # outer iterations
        ninner=10^6,            # inner iterations
        gtol=1e-3,              # tolerance on gradient for convergence of inner problem
        dtol=1e-3,              # tolerance on distance for convergence of outer problem
        rtol=0.0,               # use strict distance criteria
        lambda=1e-3,            # regularization level
        tol=1e-6,               # convergence tolerance in initialization
        maxiter=10^4,           # maximum iterations in initialization
        nesterov_threshold=100, # delay on Nesterov acceleration
        show_progress=true,     # display progress over replicates
        kwargs...               # propagate other keywords
    )

    # Write results to disk.
    @info "Writing to disk"
    traversal = sparse2dense ? "S2D" : "D2S"
    filename = joinpath(dir, "$(example)-L-path=$(traversal).dat")
    kernel = "RBF"
    problem_info = (example, subsets, p, c, sparse2dense, kernel)
    MVDA.save_cv_results(filename, problem_info, replicates, grids)

    return nothing
end

# Gaussian Clouds
ex1 = function()
    n_cv, n_test = 250, 10^3
    nsamples = n_cv + n_test
    nclasses = 3
    data = MVDA.simulate_gaussian_clouds(nsamples, nclasses; sigma=0.25, rng=StableRNG(1903))
    (MMSVD(), data, "clouds", 5, n_cv / nsamples)
end

# Nested Circles
ex2 = function()
    n_cv, n_test = 250, 10^3
    nsamples = n_cv + n_test
    nclasses = 3
    data = MVDA.simulate_nested_circles(nsamples, nclasses; p=8//10, rng=StableRNG(1903))
    (MMSVD(), data, "circles", 5, n_cv / nsamples)
end

# Waveform
ex3 = function()
    n_cv, n_test = 375, 10^3
    nsamples = n_cv + n_test
    nfeatures = 21
    data = MVDA.simulate_waveform(nsamples, nfeatures; rng=StableRNG(1903))
    (MMSVD(), data, "waveform", 5, n_cv / nsamples)
end

# Zoo
ex4 = function()
    df = MVDA.dataset("zoo")
    data = (Vector(df[!,1]), Matrix{Float64}(df[!,2:end]))
    (MMSVD(), data, "zoo", 3, 0.9)
end

# Vowel
ex5 = function()
    df = MVDA.dataset("vowel")
    data = (Vector(df[!,1]), Matrix{Float64}(df[!,2:end]))
    (MMSVD(), data, "vowel", 5, 528 / 990)
end

# HAR
# ex6 = function()
#     df = MVDA.dataset("HAR")
#     data = (Vector(df[!,1]), Matrix{Float64}(df[!,2:end]))
#     (SD(), data, "HAR", 5, 7352 / 10299)
# end

examples = (ex1, ex2, ex3, ex4, ex5,)
dir = ARGS[1]
@info "Output directory: $(dir)"

for f in examples
    algorithm, data, example, nfolds, split = f()
    run(dir, example, algorithm, data, false;
        at=split,           # CV set / Test set split
        nfolds=nfolds,      # number of folds
        nreplicates=50,     # number of CV replicates
    )
end


text:
## cal𝗚
function cal𝗚!(ap::ReproducingKernel)
    𝓖 = ap.𝓖
    𝗚 = ap.𝗠[:∇̃]
    n = get𝑛𝒒(ap)
    fill!(𝗚,0.0)
    for ξ in 𝓖
        w = get𝑤(ap,ξ)
        𝒒 = get𝒒(ap,ξ)
        for I in 1:n
            for J in I:n
                𝗚[I,J] += w*𝒒[I]*𝒒[J]
            end
        end
    end
    cholesky!(𝗚)
    U⁻¹ = inverse!(𝗚)
    𝗚⁻¹ = UUᵀ!(U⁻¹)
    return 𝗚⁻¹
end

function cal𝗚!(ap::ReproducingKernel{𝝃,:Quadratic1D,𝑠,𝜙,:Seg2}) where {𝝃<:AbstractNode,𝑠,𝜙}
    𝗚⁻¹ = ap.𝗠[:∇̃]
    fill!(𝗚⁻¹,0.0)
    𝐿 = get𝐿(ap)
    𝗚⁻¹[1] =  4.0/𝐿
    𝗚⁻¹[2] = -6.0/𝐿
    𝗚⁻¹[3] = 12.0/𝐿
    return 𝗚⁻¹
end

function cal𝗚!(ap::ReproducingKernel{𝝃,:Cubic1D,𝑠,𝜙,:Seg2}) where {𝝃<:AbstractNode,𝑠,𝜙}
    𝗚⁻¹ = ap.𝗠[:∇̃]
    fill!(𝗚⁻¹,0.0)
    𝐿 = get𝐿(ap)
    𝗚⁻¹[1] =    9.0/𝐿
    𝗚⁻¹[2] =  -36.0/𝐿
    𝗚⁻¹[3] =  192.0/𝐿
    𝗚⁻¹[4] =   30.0/𝐿
    𝗚⁻¹[5] = -180.0/𝐿
    𝗚⁻¹[6] =  180.0/𝐿
    return 𝗚⁻¹
end

function cal𝗚!(ap::ReproducingKernel{𝝃,:Linear2D,𝑠,𝜙,:Tri3}) where {𝝃<:AbstractNode,𝑠,𝜙}
    𝗚⁻¹ = ap.𝗠[:∇̃]
    fill!(𝗚⁻¹,0.0)
    𝐴 = get𝐴(ap)
    𝗚⁻¹[1] = 1.0/𝐴
    return 𝗚⁻¹
end

function cal𝗚!(ap::ReproducingKernel{𝝃,:Quadratic2D,𝑠,𝜙,:Tri3}) where {𝝃<:AbstractNode,𝑠,𝜙}
    𝗚⁻¹ = ap.𝗠[:∇̃]
    fill!(𝗚⁻¹,0.0)
    𝐴 = get𝐴(ap)
    𝗚⁻¹[1] =   9.0/𝐴
    𝗚⁻¹[2] = -12.0/𝐴
    𝗚⁻¹[3] =  24.0/𝐴
    𝗚⁻¹[4] = -12.0/𝐴
    𝗚⁻¹[5] =  12.0/𝐴
    𝗚⁻¹[6] =  24.0/𝐴
    return 𝗚⁻¹
end

function cal𝗚!(ap::ReproducingKernel{𝝃,:Cubic2D,𝑠,𝜙,:Tri3}) where {𝝃<:AbstractNode,𝑠,𝜙}
    𝗚⁻¹ = ap.𝗠[:∇̃]
    fill!(𝗚⁻¹,0.0)
    𝐴 = get𝐴(ap)
    𝗚⁻¹[1] =   36.0/𝐴
    𝗚⁻¹[2] = -120.0/𝐴
    𝗚⁻¹[3] =  600.0/𝐴
    𝗚⁻¹[4] = -120.0/𝐴
    𝗚⁻¹[5] =  300.0/𝐴
    𝗚⁻¹[6] =  600.0/𝐴
    𝗚⁻¹[7] =   90.0/𝐴
    𝗚⁻¹[8] = -540.0/𝐴
    𝗚⁻¹[9] = -180.0/𝐴
    𝗚⁻¹[10] =  540.0/𝐴
    𝗚⁻¹[11] =  180.0/𝐴
    𝗚⁻¹[12] = -720.0/𝐴
    𝗚⁻¹[13] = -720.0/𝐴
    𝗚⁻¹[14] =  540.0/𝐴
    𝗚⁻¹[15] = 1440.0/𝐴
    𝗚⁻¹[16] =   90.0/𝐴
    𝗚⁻¹[17] = -180.0/𝐴
    𝗚⁻¹[18] = -540.0/𝐴
    𝗚⁻¹[19] =   90.0/𝐴
    𝗚⁻¹[20] =  540.0/𝐴
    𝗚⁻¹[21] =  540.0/𝐴
    return 𝗚⁻¹
end

function cal𝗚₂!(ap::ReproducingKernel{𝝃,:Quadratic2D,𝑠,𝜙,:Tri3}) where {𝝃<:AbstractNode,𝑠,𝜙}
    𝗚⁻¹ = ap.𝗠[:∇̃²]
    fill!(𝗚⁻¹,0.0)
    𝐴 = get𝐴(ap)
    𝗚⁻¹[1] = 1.0/𝐴
    return 𝗚⁻¹
end
function cal𝗚₂!(ap::ReproducingKernel{𝝃,:Cubic2D,𝑠,𝜙,:Tri3}) where {𝝃<:AbstractNode,𝑠,𝜙}
    𝗚⁻¹ = ap.𝗠[:∇̃²]
    fill!(𝗚⁻¹,0.0)
    𝐴 = get𝐴(ap)
    𝗚⁻¹[1] =   9.0/𝐴
    𝗚⁻¹[2] = -12.0/𝐴
    𝗚⁻¹[3] =  24.0/𝐴
    𝗚⁻¹[4] = -12.0/𝐴
    𝗚⁻¹[5] =  12.0/𝐴
    𝗚⁻¹[6] =  24.0/𝐴
    return 𝗚⁻¹
end

function cal𝗚₂!(ap::ReproducingKernel{𝝃,:Quartic2D,𝑠,𝜙,:Tri3}) where {𝝃<:AbstractNode,𝑠,𝜙}
    𝗚⁻¹ = ap.𝗠[:∇̃²]
    fill!(𝗚⁻¹,0.0)
    𝐴 = get𝐴(ap)
    𝗚⁻¹[1] =   36.0/𝐴
    𝗚⁻¹[2] = -120.0/𝐴
    𝗚⁻¹[3] =  600.0/𝐴
    𝗚⁻¹[4] = -120.0/𝐴
    𝗚⁻¹[5] =  300.0/𝐴
    𝗚⁻¹[6] =  600.0/𝐴
    𝗚⁻¹[7] =   90.0/𝐴
    𝗚⁻¹[8] = -540.0/𝐴
    𝗚⁻¹[9] = -180.0/𝐴
    𝗚⁻¹[10] =  540.0/𝐴
    𝗚⁻¹[11] =  180.0/𝐴
    𝗚⁻¹[12] = -720.0/𝐴
    𝗚⁻¹[13] = -720.0/𝐴
    𝗚⁻¹[14] =  540.0/𝐴
    𝗚⁻¹[15] = 1440.0/𝐴
    𝗚⁻¹[16] =   90.0/𝐴
    𝗚⁻¹[17] = -180.0/𝐴
    𝗚⁻¹[18] = -540.0/𝐴
    𝗚⁻¹[19] =   90.0/𝐴
    𝗚⁻¹[20] =  540.0/𝐴
    𝗚⁻¹[21] =  540.0/𝐴
    return 𝗚⁻¹
end

## get shape for SNode
function get∇∇̃²𝝭(ap::ReproducingKernel,ξ::SNode)
    𝝭 = ap.𝝭[:∂1]
    ∂∂²𝝭∂x²∂x = ap.𝝭[:∂x²∂x]
    ∂∂²𝝭∂x²∂y = ap.𝝭[:∂x²∂y]
    ∂∂²𝝭∂x∂y∂x = ap.𝝭[:∂x∂y∂x]
    ∂∂²𝝭∂x∂y∂y = ap.𝝭[:∂x∂y∂y]
    ∂∂²𝝭∂y²∂x = ap.𝝭[:∂y²∂x]
    ∂∂²𝝭∂y²∂y = ap.𝝭[:∂y²∂y]
    i = ξ.id
    index = ξ.index
    for j in 1:length(ap.𝓒)
        𝝭[j] = ξ.𝝭[:∂1][index[i]+j]
        ∂∂²𝝭∂x²∂x[j] = ξ.𝝭[:∂x²∂x][index[i]+j]
        ∂∂²𝝭∂x²∂y[j] = ξ.𝝭[:∂x²∂y][index[i]+j]
        ∂∂²𝝭∂x∂y∂x[j] = ξ.𝝭[:∂x∂y∂x][index[i]+j]
        ∂∂²𝝭∂x∂y∂y[j] = ξ.𝝭[:∂x∂y∂y][index[i]+j]
        ∂∂²𝝭∂y²∂x[j] = ξ.𝝭[:∂y²∂x][index[i]+j]
        ∂∂²𝝭∂y²∂y[j] = ξ.𝝭[:∂y²∂y][index[i]+j]
    end
    return 𝝭, ∂∂²𝝭∂x²∂x, ∂∂²𝝭∂x²∂y, ∂∂²𝝭∂x∂y∂x, ∂∂²𝝭∂x∂y∂y, ∂∂²𝝭∂y²∂x, ∂∂²𝝭∂y²∂y
end

function get∇̄𝝭(ap::ReproducingKernel,ξ::SNode)
    𝝭 = ap.𝝭[:∂1]
    ∂𝝭∂x = ap.𝝭[:∂x]
    ∂𝝭∂y = ap.𝝭[:∂y]
    ∂𝝭∂z = ap.𝝭[:∂z]
    i = ξ.id
    index = ξ.index
    for j in 1:length(ap.𝓒)
        𝝭[j] = ξ.𝝭[:∂1][index[i]+j]
        ∂𝝭∂x[j] = ξ.𝝭[:∂̄x][index[i]+j]
        ∂𝝭∂y[j] = haskey(ξ.𝝭,:∂̄y) ? ξ.𝝭[:∂̄y][index[i]+j] : 0.0
        ∂𝝭∂z[j] = haskey(ξ.𝝭,:∂̄z) ? ξ.𝝭[:∂̄z][index[i]+j] : 0.0
    end
    return 𝝭, ∂𝝭∂x, ∂𝝭∂y, ∂𝝭∂z
end
## RK gradient smoothing
function set∇̃𝝭!(aps::Vector{T}) where T<:ReproducingKernel
    for ap in aps
        set∇̃𝝭!(ap)
    end
end
set∇̃𝝭!(ap::T) where T<:ReproducingKernel{SNode} = set∇̃𝝭!(ap,ap)

function set∇̃²𝝭!(aps::Vector{T}) where T<:ReproducingKernel
    for ap in aps
        set∇̃²𝝭!(ap)
    end
end
set∇̃²𝝭!(ap::T) where T<:ReproducingKernel{SNode} = set∇̃²𝝭!(ap,ap)

function set∇̃𝝭!(gps::Vector{T},aps::Vector{S}) where {T<:ReproducingKernel,S<:ReproducingKernel}
    if length(gps) ≠ length(aps)
        error("Miss match element numbers")
    else
        for i in 1:length(gps)
            set∇̃𝝭!(gps[i],aps[i])
        end
    end
end

function set∇̃²𝝭!(gps::Vector{T},aps::Vector{S}) where {T<:ReproducingKernel,S<:ReproducingKernel}
    if length(gps) ≠ length(aps)
        error("Miss match element numbers")
    else
        for i in 1:length(gps)
            set∇̃²𝝭!(gps[i],aps[i])
        end
    end
end

function set∇∇̃²𝝭!(gps::Vector{T},aps::Vector{S}) where {T<:ReproducingKernel,S<:ReproducingKernel}
    if length(gps) ≠ length(aps)
        error("Miss match element numbers")
    else
        for i in 1:length(gps)
            set∇∇̃²𝝭!(gps[i],aps[i])
        end
    end
end

function set∇̃𝝭!(as::Vector{T},bs::Vector{S},cs::Vector{R}) where {T<:ReproducingKernel,S<:ReproducingKernel,R<:ReproducingKernel}
    if length(as) ≠ length(bs) || length(bs) ≠ length(cs)
        error("Miss match element numbers")
    else
        for i in 1:length(as)
            set∇̃𝝭!(as[i],bs[i],cs[i])
        end
    end
end
function set∇̃𝝭!(gp::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Seg2},ap::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Seg2}) where {𝒑,𝑠,𝜙}
    n₁ =  1.0
    n₂ = -1.0
    𝓒 = gp.𝓒
    𝓖 = gp.𝓖
    for ξ̂ in 𝓖
        𝒒̂ = get𝒒(gp,ξ̂)
        𝗚⁻¹ = cal𝗚!(gp)
        𝒒̂ᵀ𝗚⁻¹ = 𝒒̂*𝗚⁻¹
        ∂𝝭∂x = gp.𝝭[:∂x]
        fill!(∂𝝭∂x,0.0)
        for ξ in ap.𝓖
            w = ξ.w/2
            wᵇ = ξ.wᵇ
            nᵇ₁ = 0.0
            nᵇ₁ += ξ.ξ ==  1.0 ? n₁ : 0.0
            nᵇ₁ += ξ.ξ == -1.0 ? n₂ : 0.0
            𝝭 = get𝝭(ap,ξ)
            𝒒, ∂𝒒∂ξ = get∇𝒒(gp,ξ)
            W₁ = 𝒒̂ᵀ𝗚⁻¹*𝒒*nᵇ₁*wᵇ + 𝒒̂ᵀ𝗚⁻¹*∂𝒒∂ξ*n₁*w
            for i in 1:length(𝓒)
                ∂𝝭∂x[i] += 𝝭[i]*W₁
            end
        end
        for i in 1:length(𝓒)
            ξ̂.𝝭[:∂x][ξ̂.index[ξ̂.id]+i] = ∂𝝭∂x[i]
        end
    end
end

function set∇̃𝝭!(gp::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Tri3},ap::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Tri3}) where {𝒑,𝑠,𝜙}
    x₁ = gp.𝓒[1].x;y₁ = gp.𝓒[1].y
    x₂ = gp.𝓒[2].x;y₂ = gp.𝓒[2].y
    x₃ = gp.𝓒[3].x;y₃ = gp.𝓒[3].y
    n₁₁ = y₃-y₂;n₂₁ = y₁-y₃;n₃₁ = y₂-y₁
    n₁₂ = x₂-x₃;n₂₂ = x₃-x₁;n₃₂ = x₁-x₂
    𝓒 = gp.𝓒
    𝓖 = gp.𝓖
    for ξ̂ in 𝓖
        𝒒̂ = get𝒑₁(gp,ξ̂)
        𝗚⁻¹ = cal𝗚!(gp)
        𝒒̂ᵀ𝗚⁻¹ = 𝒒̂*𝗚⁻¹
        ∂𝝭∂x = gp.𝝭[:∂x]
        ∂𝝭∂y = gp.𝝭[:∂y]
        fill!(∂𝝭∂x,0.0)
        fill!(∂𝝭∂y,0.0)
        for ξ in ap.𝓖
            w = ξ.w
            wᵇ = ξ.wᵇ
            𝝭 = get𝝭(ap,ξ)
            𝒒, ∂𝒒∂ξ, ∂𝒒∂η = get∇𝒑₁(ap,ξ)
            𝒒̂ᵀ𝗚⁻¹𝒒 =  𝒒̂ᵀ𝗚⁻¹*𝒒
            𝒒̂ᵀ𝗚⁻¹∂𝒒∂ξ = 𝒒̂ᵀ𝗚⁻¹*∂𝒒∂ξ
            𝒒̂ᵀ𝗚⁻¹∂𝒒∂η = 𝒒̂ᵀ𝗚⁻¹*∂𝒒∂η
            nᵇ₁ = 0.0;nᵇ₂ = 0.0
            ξ.ξ == 0.0 ? (nᵇ₁ += n₁₁;nᵇ₂ += n₁₂) : nothing
            ξ.η == 0.0 ? (nᵇ₁ += n₂₁;nᵇ₂ += n₂₂) : nothing
            ξ.ξ+ξ.η ≈ 1.0 ? (nᵇ₁ += n₃₁;nᵇ₂ += n₃₂) : nothing
            b₁ = 𝒒̂ᵀ𝗚⁻¹∂𝒒∂ξ*n₁₁ + 𝒒̂ᵀ𝗚⁻¹∂𝒒∂η*n₂₁
            b₂ = 𝒒̂ᵀ𝗚⁻¹∂𝒒∂ξ*n₁₂ + 𝒒̂ᵀ𝗚⁻¹∂𝒒∂η*n₂₂
            W₁ = 𝒒̂ᵀ𝗚⁻¹𝒒*nᵇ₁*wᵇ + b₁*w/2
            W₂ = 𝒒̂ᵀ𝗚⁻¹𝒒*nᵇ₂*wᵇ + b₂*w/2
            for i in 1:length(𝓒)
                ∂𝝭∂x[i] += 𝝭[i]*W₁
                ∂𝝭∂y[i] += 𝝭[i]*W₂
            end
        end
        for i in 1:length(𝓒)
            ξ̂.𝝭[:∂x][ξ̂.index[ξ̂.id]+i] = ∂𝝭∂x[i]
            ξ̂.𝝭[:∂y][ξ̂.index[ξ̂.id]+i] = ∂𝝭∂y[i]
        end
    end
end

function set∇̃²𝝭!(gp::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Tri3},ap::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Tri3}) where {𝒑,𝑠,𝜙}
    x₁ = gp.𝓒[1].x;y₁ = gp.𝓒[1].y
    x₂ = gp.𝓒[2].x;y₂ = gp.𝓒[2].y
    x₃ = gp.𝓒[3].x;y₃ = gp.𝓒[3].y
    𝐴 = get𝐴(gp)
    n₁₁ = y₃-y₂;n₂₁ = y₁-y₃;n₃₁ = y₂-y₁
    n₁₂ = x₂-x₃;n₂₂ = x₃-x₁;n₃₂ = x₁-x₂
    s₁₁ = -n₁₂;s₂₁ = -n₂₂;s₃₁ = -n₃₂
    s₁₂ =  n₁₁;s₂₂ =  n₂₁;s₃₂ =  n₃₁
    𝓒 = gp.𝓒
    𝓖 = gp.𝓖
    for ξ̂ in 𝓖
        𝒒̂ = get𝒑₂(gp,ξ̂)
        𝗚⁻¹ = cal𝗚₂!(gp)
        𝒒̂ᵀ𝗚⁻¹ = 𝒒̂*𝗚⁻¹
        ∂²𝝭∂x² = gp.𝝭[:∂x²]
        ∂²𝝭∂x∂y = gp.𝝭[:∂x∂y]
        ∂²𝝭∂y² = gp.𝝭[:∂y²]
        fill!(∂²𝝭∂x²,0.0)
        fill!(∂²𝝭∂x∂y,0.0)
        fill!(∂²𝝭∂y²,0.0)
        for ξ in ap.𝓖
            w = ξ.w
            wᵇ = ξ.wᵇ
            𝝭,∂𝝭∂x,∂𝝭∂y = get∇𝝭(ap,ξ)
            𝒒, ∂𝒒∂ξ, ∂𝒒∂η, ∂²𝒒∂ξ², ∂²𝒒∂ξ∂η, ∂²𝒒∂η² = get∇²𝒑₂(ap,ξ)
            𝒒̂ᵀ𝗚⁻¹𝒒 =  𝒒̂ᵀ𝗚⁻¹*𝒒
            𝒒̂ᵀ𝗚⁻¹∂𝒒∂ξ = 𝒒̂ᵀ𝗚⁻¹*∂𝒒∂ξ
            𝒒̂ᵀ𝗚⁻¹∂𝒒∂η = 𝒒̂ᵀ𝗚⁻¹*∂𝒒∂η
            𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ² = 𝒒̂ᵀ𝗚⁻¹*∂²𝒒∂ξ²
            𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ∂η = 𝒒̂ᵀ𝗚⁻¹*∂²𝒒∂ξ∂η
            𝒒̂ᵀ𝗚⁻¹∂²𝒒∂η² = 𝒒̂ᵀ𝗚⁻¹*∂²𝒒∂η²

            q₁₁ = 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ²*n₁₁*n₁₁ + 2*𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ∂η*n₁₁*n₂₁ + 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂η²*n₂₁*n₂₁
            q₁₂ = 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ²*n₁₁*n₁₂ + 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ∂η*(n₁₁*n₂₂+n₁₂*n₂₁) + 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂η²*n₂₁*n₂₂
            q₂₂ = 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ²*n₁₂*n₁₂ + 2*𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ∂η*n₁₂*n₂₂ + 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂η²*n₂₂*n₂₂

            q₁ = 𝒒̂ᵀ𝗚⁻¹∂𝒒∂ξ*n₁₁ + 𝒒̂ᵀ𝗚⁻¹∂𝒒∂η*n₂₁
            q₂ = 𝒒̂ᵀ𝗚⁻¹∂𝒒∂ξ*n₁₂ + 𝒒̂ᵀ𝗚⁻¹∂𝒒∂η*n₂₂

            q₁n₁ = 0.0;q₂n₂ = 0.0;q₁n₂ = 0.0;q₂n₁ = 0.0
            mn₁₁n₁ = 0.0;mn₁₁n₂ = 0.0;mn₁₂n₁ = 0.0;mn₁₂n₂ = 0.0;mn₂₂n₁ = 0.0;mn₂₂n₂ = 0.0
            ms₁₁ = 0.0;ms₁₂ = 0.0;ms₂₂ = 0.0
            Δms₁₁ = 0.0;Δms₁₂ = 0.0;Δms₂₂ = 0.0
            𝐿₁² = n₁₁^2+n₁₂^2
            𝐿₂² = n₂₁^2+n₂₂^2
            𝐿₃² = n₃₁^2+n₃₂^2
            if ξ.ξ == 0.0
                q₁n₁ += q₁*n₁₁
                q₁n₂ += q₁*n₁₂
                q₂n₁ += q₂*n₁₁
                q₂n₂ += q₂*n₁₂
                mn₁₁n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₁₁*n₁₁*n₁₁/𝐿₁²
                mn₁₁n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₁₁*n₁₁*n₁₂/𝐿₁²
                mn₁₂n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₁₁*n₁₂*n₁₁/𝐿₁²
                mn₁₂n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₁₁*n₁₂*n₁₂/𝐿₁²
                mn₂₂n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₁₂*n₁₂*n₁₁/𝐿₁²
                mn₂₂n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₁₂*n₁₂*n₁₂/𝐿₁²
                ms₁₁ += (q₁*s₁₁+q₂*s₁₂)*n₁₁*s₁₁/𝐿₁²
                ms₁₂ += (q₁*s₁₁+q₂*s₁₂)*0.5*(n₁₁*s₁₂+n₁₂*s₁₁)/𝐿₁²
                ms₂₂ += (q₁*s₁₁+q₂*s₁₂)*n₁₂*s₁₂/𝐿₁²
                if ξ.η == 0.0
                    Δms₁₁ =  𝒒̂ᵀ𝗚⁻¹𝒒*(n₁₁*s₁₁/𝐿₁²-n₂₁*s₂₁/𝐿₂²)
                    Δms₂₂ =  𝒒̂ᵀ𝗚⁻¹𝒒*(n₁₂*s₁₂/𝐿₁²-n₂₂*s₂₂/𝐿₂²)
                    Δms₁₂ =  𝒒̂ᵀ𝗚⁻¹𝒒*0.5*((n₁₁*s₁₂+n₁₂*s₁₁)/𝐿₁²-(n₂₁*s₂₂+n₂₂*s₂₁)/𝐿₂²)
                end
            end
            if  ξ.η == 0.0
                q₁n₁ += q₁*n₂₁
                q₁n₂ += q₁*n₂₂
                q₂n₁ += q₂*n₂₁
                q₂n₂ += q₂*n₂₂
                mn₁₁n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₂₁*n₂₁*n₂₁/𝐿₂²
                mn₁₁n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₂₁*n₂₁*n₂₂/𝐿₂²
                mn₁₂n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₂₁*n₂₂*n₂₁/𝐿₂²
                mn₁₂n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₂₁*n₂₂*n₂₂/𝐿₂²
                mn₂₂n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₂₂*n₂₂*n₂₁/𝐿₂²
                mn₂₂n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₂₂*n₂₂*n₂₂/𝐿₂²
                ms₁₁ += (q₁*s₂₁+q₂*s₂₂)*n₂₁*s₂₁/𝐿₂²
                ms₁₂ += (q₁*s₂₁+q₂*s₂₂)*0.5*(n₂₁*s₂₂+n₂₂*s₂₁)/𝐿₂²
                ms₂₂ += (q₁*s₂₁+q₂*s₂₂)*n₂₂*s₂₂/𝐿₂²
                if ξ.ξ+ξ.η ≈ 1.0
                    Δms₁₁ =  𝒒̂ᵀ𝗚⁻¹𝒒*(n₂₁*s₂₁/𝐿₂²-n₃₁*s₃₁/𝐿₃²)
                    Δms₂₂ =  𝒒̂ᵀ𝗚⁻¹𝒒*(n₂₂*s₂₂/𝐿₂²-n₃₂*s₃₂/𝐿₃²)
                    Δms₁₂ =  𝒒̂ᵀ𝗚⁻¹𝒒*0.5*((n₂₁*s₂₂+n₂₂*s₂₁)/𝐿₂²-(n₃₁*s₃₂+n₃₂*s₃₁)/𝐿₃²)
                end
            end
            if ξ.ξ+ξ.η ≈ 1.0
                q₁n₁ += q₁*n₃₁
                q₁n₂ += q₁*n₃₂
                q₂n₁ += q₂*n₃₁
                q₂n₂ += q₂*n₃₂
                mn₁₁n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₃₁*n₃₁*n₃₁/𝐿₃²
                mn₁₁n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₃₁*n₃₁*n₃₂/𝐿₃²
                mn₁₂n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₃₁*n₃₂*n₃₁/𝐿₃²
                mn₁₂n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₃₁*n₃₂*n₃₂/𝐿₃²
                mn₂₂n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₃₂*n₃₂*n₃₁/𝐿₃²
                mn₂₂n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₃₂*n₃₂*n₃₂/𝐿₃²
                ms₁₁ += (q₁*s₃₁+q₂*s₃₂)*n₃₁*s₃₁/𝐿₃²
                ms₁₂ += (q₁*s₃₁+q₂*s₃₂)*0.5*(n₃₁*s₃₂+n₃₂*s₃₁)/𝐿₃²
                ms₂₂ += (q₁*s₃₁+q₂*s₃₂)*n₃₂*s₃₂/𝐿₃²
                if ξ.ξ == 0.0
                    Δms₁₁ =  𝒒̂ᵀ𝗚⁻¹𝒒*(n₃₁*s₃₁/𝐿₃²-n₁₁*s₁₁/𝐿₁²)
                    Δms₂₂ =  𝒒̂ᵀ𝗚⁻¹𝒒*(n₃₂*s₃₂/𝐿₃²-n₁₂*s₁₂/𝐿₁²)
                    Δms₁₂ =  𝒒̂ᵀ𝗚⁻¹𝒒*0.5*((n₃₁*s₃₂+n₃₂*s₃₁)/𝐿₃²-(n₁₁*s₁₂+n₁₂*s₁₁)/𝐿₁²)
                end
            end

            W₁₁₁ = mn₁₁n₁*wᵇ
            W₁₁₂ = mn₁₁n₂*wᵇ
            W₁₂₁ = mn₁₂n₁*wᵇ
            W₁₂₂ = mn₁₂n₂*wᵇ
            W₂₂₁ = mn₂₂n₁*wᵇ
            W₂₂₂ = mn₂₂n₂*wᵇ
            W₁₁ = (q₁₁*w + 2*(q₁n₁+ms₁₁)*wᵇ)/4/𝐴 + Δms₁₁
            W₁₂ = (q₁₂*w + (q₁n₂+q₂n₁+2*ms₁₂)*wᵇ)/4/𝐴 + Δms₁₂
            W₂₂ = (q₂₂*w + 2*(q₂n₂+ms₂₂)*wᵇ)/4/𝐴 + Δms₂₂
            for i in 1:length(𝓒)
                ∂²𝝭∂x²[i] += 𝝭[i]*W₁₁ + ∂𝝭∂x[i]*W₁₁₁ + ∂𝝭∂y[i]*W₁₁₂
                ∂²𝝭∂x∂y[i] += 𝝭[i]*W₁₂ + ∂𝝭∂x[i]*W₁₂₁ + ∂𝝭∂y[i]*W₁₂₂
                ∂²𝝭∂y²[i] += 𝝭[i]*W₂₂ + ∂𝝭∂x[i]*W₂₂₁ + ∂𝝭∂y[i]*W₂₂₂
            end
        end
        for i in 1:length(𝓒)
            ξ̂.𝝭[:∂x²][ξ̂.index[ξ̂.id]+i] = ∂²𝝭∂x²[i]
            ξ̂.𝝭[:∂x∂y][ξ̂.index[ξ̂.id]+i] = ∂²𝝭∂x∂y[i]
            ξ̂.𝝭[:∂y²][ξ̂.index[ξ̂.id]+i] = ∂²𝝭∂y²[i]
        end
    end
end

function set∇∇̃²𝝭!(gp::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Seg2},ap::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Tri3}) where {𝒑,𝑠,𝜙}
    x₁ = ap.𝓒[1].x;y₁ = ap.𝓒[1].y
    x₂ = ap.𝓒[2].x;y₂ = ap.𝓒[2].y
    x₃ = ap.𝓒[3].x;y₃ = ap.𝓒[3].y
    𝐴 = get𝐴(ap)
    n₁₁ = y₃-y₂;n₂₁ = y₁-y₃;n₃₁ = y₂-y₁
    n₁₂ = x₂-x₃;n₂₂ = x₃-x₁;n₃₂ = x₁-x₂
    s₁₁ = -n₁₂;s₂₁ = -n₂₂;s₃₁ = -n₃₂
    s₁₂ =  n₁₁;s₂₂ =  n₂₁;s₃₂ =  n₃₁
    𝐿₁² = n₁₁^2+n₁₂^2
    𝐿₂² = n₂₁^2+n₂₂^2
    𝐿₃² = n₃₁^2+n₃₂^2
    𝓒 = gp.𝓒
    𝓖 = gp.𝓖
    for ξ̂ in 𝓖
        𝒒̂,∂𝒒̂∂ξ,∂𝒒̂∂η = get∇𝒑₂(ap,ξ̂)
        𝗚⁻¹ = cal𝗚₂!(ap)
        𝒒̂ᵀ𝗚⁻¹ = 𝒒̂*𝗚⁻¹
        ∂𝒒̂ᵀ∂ξ𝗚⁻¹ = ∂𝒒̂∂ξ*𝗚⁻¹
        ∂𝒒̂ᵀ∂η𝗚⁻¹ = ∂𝒒̂∂η*𝗚⁻¹
        ∂²𝝭∂x² = gp.𝝭[:∂x²]
        ∂²𝝭∂x∂y = gp.𝝭[:∂x∂y]
        ∂²𝝭∂y² = gp.𝝭[:∂y²]
        ∂∂²𝝭∂x²∂x = gp.𝝭[:∂x²∂x]
        ∂∂²𝝭∂x²∂y = gp.𝝭[:∂x²∂y]
        ∂∂²𝝭∂x∂y∂x = gp.𝝭[:∂x∂y∂x]
        ∂∂²𝝭∂x∂y∂y = gp.𝝭[:∂x∂y∂y]
        ∂∂²𝝭∂y²∂x = gp.𝝭[:∂y²∂x]
        ∂∂²𝝭∂y²∂y = gp.𝝭[:∂y²∂y]
        fill!(∂²𝝭∂x²,0.0)
        fill!(∂²𝝭∂x∂y,0.0)
        fill!(∂²𝝭∂y²,0.0)
        fill!(∂∂²𝝭∂x²∂x,0.0)
        fill!(∂∂²𝝭∂x²∂y,0.0)
        fill!(∂∂²𝝭∂x∂y∂x,0.0)
        fill!(∂∂²𝝭∂x∂y∂y,0.0)
        fill!(∂∂²𝝭∂y²∂x,0.0)
        fill!(∂∂²𝝭∂y²∂y,0.0)
        for ξ in ap.𝓖
            w = ξ.w
            wᵇ = ξ.wᵇ
            𝝭,∂𝝭∂x,∂𝝭∂y = get∇𝝭(ap,ξ)
            𝒒, ∂𝒒∂ξ, ∂𝒒∂η, ∂²𝒒∂ξ², ∂²𝒒∂ξ∂η, ∂²𝒒∂η² = get∇²𝒑₂(ap,ξ)

            𝒒̂ᵀ𝗚⁻¹𝒒 =  𝒒̂ᵀ𝗚⁻¹*𝒒
            ∂𝒒̂ᵀ∂ξ𝗚⁻¹𝒒 =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹*𝒒
            ∂𝒒̂ᵀ∂η𝗚⁻¹𝒒 =  ∂𝒒̂ᵀ∂η𝗚⁻¹*𝒒
            ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒 =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹𝒒*n₁₁+∂𝒒̂ᵀ∂η𝗚⁻¹𝒒*n₂₁
            ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒 =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹𝒒*n₁₂+∂𝒒̂ᵀ∂η𝗚⁻¹𝒒*n₂₂

            𝒒̂ᵀ𝗚⁻¹∂𝒒∂ξ = 𝒒̂ᵀ𝗚⁻¹*∂𝒒∂ξ
            ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂ξ = ∂𝒒̂ᵀ∂ξ𝗚⁻¹*∂𝒒∂ξ
            ∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂ξ = ∂𝒒̂ᵀ∂η𝗚⁻¹*∂𝒒∂ξ
            ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂ξ =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂ξ*n₁₁+∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂ξ*n₂₁
            ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂ξ =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂ξ*n₁₂+∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂ξ*n₂₂

            𝒒̂ᵀ𝗚⁻¹∂𝒒∂η = 𝒒̂ᵀ𝗚⁻¹*∂𝒒∂η
            ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂η = ∂𝒒̂ᵀ∂ξ𝗚⁻¹*∂𝒒∂η
            ∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂η = ∂𝒒̂ᵀ∂η𝗚⁻¹*∂𝒒∂η
            ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂η =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂η*n₁₁+∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂η*n₂₁
            ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂η =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂η*n₁₂+∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂η*n₂₂

            𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ² = 𝒒̂ᵀ𝗚⁻¹*∂²𝒒∂ξ²
            ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂²𝒒∂ξ² = ∂𝒒̂ᵀ∂ξ𝗚⁻¹*∂²𝒒∂ξ²
            ∂𝒒̂ᵀ∂η𝗚⁻¹∂²𝒒∂ξ² = ∂𝒒̂ᵀ∂η𝗚⁻¹*∂²𝒒∂ξ²
            ∂𝒒̂ᵀ∂x𝗚⁻¹∂²𝒒∂ξ² =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂²𝒒∂ξ²*n₁₁+∂𝒒̂ᵀ∂η𝗚⁻¹∂²𝒒∂ξ²*n₂₁
            ∂𝒒̂ᵀ∂y𝗚⁻¹∂²𝒒∂ξ² =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂²𝒒∂ξ²*n₁₂+∂𝒒̂ᵀ∂η𝗚⁻¹∂²𝒒∂ξ²*n₂₂

            𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ∂η = 𝒒̂ᵀ𝗚⁻¹*∂²𝒒∂ξ∂η
            ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂²𝒒∂ξ∂η = ∂𝒒̂ᵀ∂ξ𝗚⁻¹*∂²𝒒∂ξ∂η
            ∂𝒒̂ᵀ∂η𝗚⁻¹∂²𝒒∂ξ∂η = ∂𝒒̂ᵀ∂η𝗚⁻¹*∂²𝒒∂ξ∂η
            ∂𝒒̂ᵀ∂x𝗚⁻¹∂²𝒒∂ξ∂η =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂²𝒒∂ξ∂η*n₁₁+∂𝒒̂ᵀ∂η𝗚⁻¹∂²𝒒∂ξ∂η*n₂₁
            ∂𝒒̂ᵀ∂y𝗚⁻¹∂²𝒒∂ξ∂η =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂²𝒒∂ξ∂η*n₁₂+∂𝒒̂ᵀ∂η𝗚⁻¹∂²𝒒∂ξ∂η*n₂₂

            𝒒̂ᵀ𝗚⁻¹∂²𝒒∂η² = 𝒒̂ᵀ𝗚⁻¹*∂²𝒒∂η²
            ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂²𝒒∂η² = ∂𝒒̂ᵀ∂ξ𝗚⁻¹*∂²𝒒∂η²
            ∂𝒒̂ᵀ∂η𝗚⁻¹∂²𝒒∂η² = ∂𝒒̂ᵀ∂η𝗚⁻¹*∂²𝒒∂η²
            ∂𝒒̂ᵀ∂x𝗚⁻¹∂²𝒒∂η² =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂²𝒒∂η²*n₁₁+∂𝒒̂ᵀ∂η𝗚⁻¹∂²𝒒∂η²*n₂₁
            ∂𝒒̂ᵀ∂y𝗚⁻¹∂²𝒒∂η² =  ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂²𝒒∂η²*n₁₂+∂𝒒̂ᵀ∂η𝗚⁻¹∂²𝒒∂η²*n₂₂

            q₁₁ = 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ²*n₁₁*n₁₁ + 2*𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ∂η*n₁₁*n₂₁ + 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂η²*n₂₁*n₂₁
            ∂q₁₁∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹∂²𝒒∂ξ²*n₁₁*n₁₁ + 2*∂𝒒̂ᵀ∂x𝗚⁻¹∂²𝒒∂ξ∂η*n₁₁*n₂₁ + ∂𝒒̂ᵀ∂x𝗚⁻¹∂²𝒒∂η²*n₂₁*n₂₁
            ∂q₁₁∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹∂²𝒒∂ξ²*n₁₁*n₁₁ + 2*∂𝒒̂ᵀ∂y𝗚⁻¹∂²𝒒∂ξ∂η*n₁₁*n₂₁ + ∂𝒒̂ᵀ∂y𝗚⁻¹∂²𝒒∂η²*n₂₁*n₂₁

            q₁₂ = 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ²*n₁₁*n₁₂ + 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ∂η*(n₁₁*n₂₂+n₁₂*n₂₁) + 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂η²*n₂₁*n₂₂
            ∂q₁₂∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹∂²𝒒∂ξ²*n₁₁*n₁₂ + ∂𝒒̂ᵀ∂x𝗚⁻¹∂²𝒒∂ξ∂η*(n₁₁*n₂₂+n₁₂*n₂₁) + ∂𝒒̂ᵀ∂x𝗚⁻¹∂²𝒒∂η²*n₂₁*n₂₂
            ∂q₁₂∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹∂²𝒒∂ξ²*n₁₁*n₁₂ + ∂𝒒̂ᵀ∂y𝗚⁻¹∂²𝒒∂ξ∂η*(n₁₁*n₂₂+n₁₂*n₂₁) + ∂𝒒̂ᵀ∂y𝗚⁻¹∂²𝒒∂η²*n₂₁*n₂₂

            q₂₂ = 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ²*n₁₂*n₁₂ + 2*𝒒̂ᵀ𝗚⁻¹∂²𝒒∂ξ∂η*n₁₂*n₂₂ + 𝒒̂ᵀ𝗚⁻¹∂²𝒒∂η²*n₂₂*n₂₂
            ∂q₂₂∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹∂²𝒒∂ξ²*n₁₂*n₁₂ + 2*∂𝒒̂ᵀ∂x𝗚⁻¹∂²𝒒∂ξ∂η*n₁₂*n₂₂ + ∂𝒒̂ᵀ∂x𝗚⁻¹∂²𝒒∂η²*n₂₂*n₂₂
            ∂q₂₂∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹∂²𝒒∂ξ²*n₁₂*n₁₂ + 2*∂𝒒̂ᵀ∂y𝗚⁻¹∂²𝒒∂ξ∂η*n₁₂*n₂₂ + ∂𝒒̂ᵀ∂y𝗚⁻¹∂²𝒒∂η²*n₂₂*n₂₂

            q₁ = 𝒒̂ᵀ𝗚⁻¹∂𝒒∂ξ*n₁₁ + 𝒒̂ᵀ𝗚⁻¹∂𝒒∂η*n₂₁
            ∂q₁∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂ξ*n₁₁ + ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂η*n₂₁
            ∂q₁∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂ξ*n₁₁ + ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂η*n₂₁
            q₂ = 𝒒̂ᵀ𝗚⁻¹∂𝒒∂ξ*n₁₂ + 𝒒̂ᵀ𝗚⁻¹∂𝒒∂η*n₂₂
            ∂q₂∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂ξ*n₁₂ + ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂η*n₂₂
            ∂q₂∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂ξ*n₁₂ + ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂η*n₂₂

            q₁n₁ = 0.0;q₂n₂ = 0.0;q₁n₂ = 0.0;q₂n₁ = 0.0
            ∂q₁∂xn₁ = 0.0;∂q₂∂xn₂ = 0.0;∂q₁∂xn₂ = 0.0;∂q₂∂xn₁ = 0.0
            ∂q₁∂yn₁ = 0.0;∂q₂∂yn₂ = 0.0;∂q₁∂yn₂ = 0.0;∂q₂∂yn₁ = 0.0
            mn₁₁n₁ = 0.0;mn₁₁n₂ = 0.0;mn₁₂n₁ = 0.0;mn₁₂n₂ = 0.0;mn₂₂n₁ = 0.0;mn₂₂n₂ = 0.0
            ∂mn₁₁∂xn₁ = 0.0;∂mn₁₁∂xn₂ = 0.0;∂mn₁₂∂xn₁ = 0.0;∂mn₁₂∂xn₂ = 0.0;∂mn₂₂∂xn₁ = 0.0;∂mn₂₂∂xn₂ = 0.0
            ∂mn₁₁∂yn₁ = 0.0;∂mn₁₁∂yn₂ = 0.0;∂mn₁₂∂yn₁ = 0.0;∂mn₁₂∂yn₂ = 0.0;∂mn₂₂∂yn₁ = 0.0;∂mn₂₂∂yn₂ = 0.0
            ms₁₁ = 0.0;ms₁₂ = 0.0;ms₂₂ = 0.0
            ∂ms₁₁∂x = 0.0;∂ms₁₂∂x = 0.0;∂ms₂₂∂x = 0.0
            ∂ms₁₁∂y = 0.0;∂ms₁₂∂y = 0.0;∂ms₂₂∂y = 0.0
            Δms₁₁ = 0.0;Δms₁₂ = 0.0;Δms₂₂ = 0.0
            Δ∂ms₁₁∂x = 0.0;Δ∂ms₁₂∂x = 0.0;Δ∂ms₂₂∂x = 0.0
            Δ∂ms₁₁∂y = 0.0;Δ∂ms₁₂∂y = 0.0;Δ∂ms₂₂∂y = 0.0
            if ξ.ξ == 0.0
                q₁n₁ += q₁*n₁₁
                ∂q₁∂xn₁ += ∂q₁∂x*n₁₁
                ∂q₁∂yn₁ += ∂q₁∂y*n₁₁
                q₁n₂ += q₁*n₁₂
                ∂q₁∂xn₂ += ∂q₁∂x*n₁₂
                ∂q₁∂yn₂ += ∂q₁∂y*n₁₂
                q₂n₁ += q₂*n₁₁
                ∂q₂∂xn₁ += ∂q₂∂x*n₁₁
                ∂q₂∂yn₁ += ∂q₂∂y*n₁₁
                q₂n₂ += q₂*n₁₂
                ∂q₂∂xn₂ += ∂q₂∂x*n₁₂
                ∂q₂∂yn₂ += ∂q₂∂y*n₁₂
                mn₁₁n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₁₁*n₁₁*n₁₁/𝐿₁²
                ∂mn₁₁∂xn₁ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₁₁*n₁₁*n₁₁/𝐿₁²
                ∂mn₁₁∂yn₁ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₁₁*n₁₁*n₁₁/𝐿₁²
                mn₁₁n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₁₁*n₁₁*n₁₂/𝐿₁²
                ∂mn₁₁∂xn₂ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₁₁*n₁₁*n₁₂/𝐿₁²
                ∂mn₁₁∂yn₂ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₁₁*n₁₁*n₁₂/𝐿₁²
                mn₁₂n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₁₁*n₁₂*n₁₁/𝐿₁²
                ∂mn₁₂∂xn₁ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₁₁*n₁₂*n₁₁/𝐿₁²
                ∂mn₁₂∂yn₁ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₁₁*n₁₂*n₁₁/𝐿₁²
                mn₁₂n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₁₁*n₁₂*n₁₂/𝐿₁²
                ∂mn₁₂∂xn₂ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₁₁*n₁₂*n₁₂/𝐿₁²
                ∂mn₁₂∂yn₂ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₁₁*n₁₂*n₁₂/𝐿₁²
                mn₂₂n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₁₂*n₁₂*n₁₁/𝐿₁²
                ∂mn₂₂∂xn₁ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₁₂*n₁₂*n₁₁/𝐿₁²
                ∂mn₂₂∂yn₁ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₁₂*n₁₂*n₁₁/𝐿₁²
                mn₂₂n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₁₂*n₁₂*n₁₂/𝐿₁²
                ∂mn₂₂∂xn₂ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₁₂*n₁₂*n₁₂/𝐿₁²
                ∂mn₂₂∂yn₂ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₁₂*n₁₂*n₁₂/𝐿₁²
                ms₁₁ += (q₁*s₁₁+q₂*s₁₂)*n₁₁*s₁₁/𝐿₁²
                ∂ms₁₁∂x += (∂q₁∂x*s₁₁+∂q₂∂x*s₁₂)*n₁₁*s₁₁/𝐿₁²
                ∂ms₁₁∂y += (∂q₁∂y*s₁₁+∂q₂∂y*s₁₂)*n₁₁*s₁₁/𝐿₁²
                ms₁₂ += (q₁*s₁₁+q₂*s₁₂)*0.5*(n₁₁*s₁₂+n₁₂*s₁₁)/𝐿₁²
                ∂ms₁₂∂x += (∂q₁∂x*s₁₁+∂q₂∂x*s₁₂)*0.5*(n₁₁*s₁₂+n₁₂*s₁₁)/𝐿₁²
                ∂ms₁₂∂y += (∂q₁∂y*s₁₁+∂q₂∂y*s₁₂)*0.5*(n₁₁*s₁₂+n₁₂*s₁₁)/𝐿₁²
                ms₂₂ += (q₁*s₁₁+q₂*s₁₂)*n₁₂*s₁₂/𝐿₁²
                ∂ms₂₂∂x += (∂q₁∂x*s₁₁+∂q₂∂x*s₁₂)*n₁₂*s₁₂/𝐿₁²
                ∂ms₂₂∂y += (∂q₁∂y*s₁₁+∂q₂∂y*s₁₂)*n₁₂*s₁₂/𝐿₁²
                if ξ.η == 0.0
                    Δms₁₁ =  𝒒̂ᵀ𝗚⁻¹𝒒*(n₁₁*s₁₁/𝐿₁²-n₂₁*s₂₁/𝐿₂²)
                    Δ∂ms₁₁∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*(n₁₁*s₁₁/𝐿₁²-n₂₁*s₂₁/𝐿₂²)
                    Δ∂ms₁₁∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*(n₁₁*s₁₁/𝐿₁²-n₂₁*s₂₁/𝐿₂²)
                    Δms₂₂ =  𝒒̂ᵀ𝗚⁻¹𝒒*(n₁₂*s₁₂/𝐿₁²-n₂₂*s₂₂/𝐿₂²)
                    Δ∂ms₂₂∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*(n₁₂*s₁₂/𝐿₁²-n₂₂*s₂₂/𝐿₂²)
                    Δ∂ms₂₂∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*(n₁₂*s₁₂/𝐿₁²-n₂₂*s₂₂/𝐿₂²)
                    Δms₁₂ =  𝒒̂ᵀ𝗚⁻¹𝒒*0.5*((n₁₁*s₁₂+n₁₂*s₁₁)/𝐿₁²-(n₂₁*s₂₂+n₂₂*s₂₁)/𝐿₂²)
                    Δ∂ms₁₂∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*0.5*((n₁₁*s₁₂+n₁₂*s₁₁)/𝐿₁²-(n₂₁*s₂₂+n₂₂*s₂₁)/𝐿₂²)
                    Δ∂ms₁₂∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*0.5*((n₁₁*s₁₂+n₁₂*s₁₁)/𝐿₁²-(n₂₁*s₂₂+n₂₂*s₂₁)/𝐿₂²)
                end
            end
            if  ξ.η == 0.0
                q₁n₁ += q₁*n₂₁
                ∂q₁∂xn₁ += ∂q₁∂x*n₂₁
                ∂q₁∂yn₁ += ∂q₁∂y*n₂₁
                q₁n₂ += q₁*n₂₂
                ∂q₁∂xn₂ += ∂q₁∂x*n₂₂
                ∂q₁∂yn₂ += ∂q₁∂y*n₂₂
                q₂n₁ += q₂*n₂₁
                ∂q₂∂xn₁ += ∂q₂∂x*n₂₁
                ∂q₂∂yn₁ += ∂q₂∂y*n₂₁
                q₂n₂ += q₂*n₂₂
                ∂q₂∂xn₂ += ∂q₂∂x*n₂₂
                ∂q₂∂yn₂ += ∂q₂∂y*n₂₂
                mn₁₁n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₂₁*n₂₁*n₂₁/𝐿₂²
                ∂mn₁₁∂xn₁ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₂₁*n₂₁*n₂₁/𝐿₂²
                ∂mn₁₁∂yn₁ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₂₁*n₂₁*n₂₁/𝐿₂²
                mn₁₁n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₂₁*n₂₁*n₂₂/𝐿₂²
                ∂mn₁₁∂xn₂ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₂₁*n₂₁*n₂₂/𝐿₂²
                ∂mn₁₁∂yn₂ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₂₁*n₂₁*n₂₂/𝐿₂²
                mn₁₂n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₂₁*n₂₂*n₂₁/𝐿₂²
                ∂mn₁₂∂xn₁ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₂₁*n₂₂*n₂₁/𝐿₂²
                ∂mn₁₂∂yn₁ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₂₁*n₂₂*n₂₁/𝐿₂²
                mn₁₂n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₂₁*n₂₂*n₂₂/𝐿₂²
                ∂mn₁₂∂xn₂ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₂₁*n₂₂*n₂₂/𝐿₂²
                ∂mn₁₂∂yn₂ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₂₁*n₂₂*n₂₂/𝐿₂²
                mn₂₂n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₂₂*n₂₂*n₂₁/𝐿₂²
                ∂mn₂₂∂xn₁ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₂₂*n₂₂*n₂₁/𝐿₂²
                ∂mn₂₂∂yn₁ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₂₂*n₂₂*n₂₁/𝐿₂²
                mn₂₂n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₂₂*n₂₂*n₂₂/𝐿₂²
                ∂mn₂₂∂xn₂ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₂₂*n₂₂*n₂₂/𝐿₂²
                ∂mn₂₂∂yn₂ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₂₂*n₂₂*n₂₂/𝐿₂²
                ms₁₁ += (q₁*s₂₁+q₂*s₂₂)*n₂₁*s₂₁/𝐿₂²
                ∂ms₁₁∂x += (∂q₁∂x*s₂₁+∂q₂∂x*s₂₂)*n₂₁*s₂₁/𝐿₂²
                ∂ms₁₁∂y += (∂q₁∂y*s₂₁+∂q₂∂y*s₂₂)*n₂₁*s₂₁/𝐿₂²
                ms₁₂ += (q₁*s₂₁+q₂*s₂₂)*0.5*(n₂₁*s₂₂+n₂₂*s₂₁)/𝐿₂²
                ∂ms₁₂∂x += (∂q₁∂x*s₂₁+∂q₂∂x*s₂₂)*0.5*(n₂₁*s₂₂+n₂₂*s₂₁)/𝐿₂²
                ∂ms₁₂∂y += (∂q₁∂y*s₂₁+∂q₂∂y*s₂₂)*0.5*(n₂₁*s₂₂+n₂₂*s₂₁)/𝐿₂²
                ms₂₂ += (q₁*s₂₁+q₂*s₂₂)*n₂₂*s₂₂/𝐿₂²
                ∂ms₂₂∂x += (∂q₁∂x*s₂₁+∂q₂∂x*s₂₂)*n₂₂*s₂₂/𝐿₂²
                ∂ms₂₂∂y += (∂q₁∂y*s₂₁+∂q₂∂y*s₂₂)*n₂₂*s₂₂/𝐿₂²
                if ξ.ξ+ξ.η ≈ 1.0
                    Δms₁₁ =  𝒒̂ᵀ𝗚⁻¹𝒒*(n₂₁*s₂₁/𝐿₂²-n₃₁*s₃₁/𝐿₃²)
                    Δ∂ms₁₁∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*(n₂₁*s₂₁/𝐿₂²-n₃₁*s₃₁/𝐿₃²)
                    Δ∂ms₁₁∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*(n₂₁*s₂₁/𝐿₂²-n₃₁*s₃₁/𝐿₃²)
                    Δms₂₂ =  𝒒̂ᵀ𝗚⁻¹𝒒*(n₂₂*s₂₂/𝐿₂²-n₃₂*s₃₂/𝐿₃²)
                    Δ∂ms₂₂∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*(n₂₂*s₂₂/𝐿₂²-n₃₂*s₃₂/𝐿₃²)
                    Δ∂ms₂₂∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*(n₂₂*s₂₂/𝐿₂²-n₃₂*s₃₂/𝐿₃²)
                    Δms₁₂ =  𝒒̂ᵀ𝗚⁻¹𝒒*0.5*((n₂₁*s₂₂+n₂₂*s₂₁)/𝐿₂²-(n₃₁*s₃₂+n₃₂*s₃₁)/𝐿₃²)
                    Δ∂ms₁₂∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*0.5*((n₂₁*s₂₂+n₂₂*s₂₁)/𝐿₂²-(n₃₁*s₃₂+n₃₂*s₃₁)/𝐿₃²)
                    Δ∂ms₁₂∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*0.5*((n₂₁*s₂₂+n₂₂*s₂₁)/𝐿₂²-(n₃₁*s₃₂+n₃₂*s₃₁)/𝐿₃²)
                end
            end
            if ξ.ξ+ξ.η ≈ 1.0
                q₁n₁ += q₁*n₃₁
                ∂q₁∂xn₁ += ∂q₁∂x*n₃₁
                ∂q₁∂yn₁ += ∂q₁∂y*n₃₁
                q₁n₂ += q₁*n₃₂
                ∂q₁∂xn₂ += ∂q₁∂x*n₃₂
                ∂q₁∂yn₂ += ∂q₁∂y*n₃₂
                q₂n₁ += q₂*n₃₁
                ∂q₂∂xn₁ += ∂q₂∂x*n₃₁
                ∂q₂∂yn₁ += ∂q₂∂y*n₃₁
                q₂n₂ += q₂*n₃₂
                ∂q₂∂xn₂ += ∂q₂∂x*n₃₂
                ∂q₂∂yn₂ += ∂q₂∂y*n₃₂
                mn₁₁n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₃₁*n₃₁*n₃₁/𝐿₃²
                ∂mn₁₁∂xn₁ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₃₁*n₃₁*n₃₁/𝐿₃²
                ∂mn₁₁∂yn₁ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₃₁*n₃₁*n₃₁/𝐿₃²
                mn₁₁n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₃₁*n₃₁*n₃₂/𝐿₃²
                ∂mn₁₁∂xn₂ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₃₁*n₃₁*n₃₂/𝐿₃²
                ∂mn₁₁∂yn₂ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₃₁*n₃₁*n₃₂/𝐿₃²
                mn₁₂n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₃₁*n₃₂*n₃₁/𝐿₃²
                ∂mn₁₂∂xn₁ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₃₁*n₃₂*n₃₁/𝐿₃²
                ∂mn₁₂∂yn₁ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₃₁*n₃₂*n₃₁/𝐿₃²
                mn₁₂n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₃₁*n₃₂*n₃₂/𝐿₃²
                ∂mn₁₂∂xn₂ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₃₁*n₃₂*n₃₂/𝐿₃²
                ∂mn₁₂∂yn₂ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₃₁*n₃₂*n₃₂/𝐿₃²
                mn₂₂n₁ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₃₂*n₃₂*n₃₁/𝐿₃²
                ∂mn₂₂∂xn₁ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₃₂*n₃₂*n₃₁/𝐿₃²
                ∂mn₂₂∂yn₁ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₃₂*n₃₂*n₃₁/𝐿₃²
                mn₂₂n₂ += 𝒒̂ᵀ𝗚⁻¹𝒒*n₃₂*n₃₂*n₃₂/𝐿₃²
                ∂mn₂₂∂xn₂ += ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*n₃₂*n₃₂*n₃₂/𝐿₃²
                ∂mn₂₂∂yn₂ += ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*n₃₂*n₃₂*n₃₂/𝐿₃²
                ms₁₁ += (q₁*s₃₁+q₂*s₃₂)*n₃₁*s₃₁/𝐿₃²
                ∂ms₁₁∂x += (∂q₁∂x*s₃₁+∂q₂∂x*s₃₂)*n₃₁*s₃₁/𝐿₃²
                ∂ms₁₁∂y += (∂q₁∂y*s₃₁+∂q₂∂y*s₃₂)*n₃₁*s₃₁/𝐿₃²
                ms₁₂ += (q₁*s₃₁+q₂*s₃₂)*0.5*(n₃₁*s₃₂+n₃₂*s₃₁)/𝐿₃²
                ∂ms₁₂∂x += (∂q₁∂x*s₃₁+∂q₂∂x*s₃₂)*0.5*(n₃₁*s₃₂+n₃₂*s₃₁)/𝐿₃²
                ∂ms₁₂∂y += (∂q₁∂y*s₃₁+∂q₂∂y*s₃₂)*0.5*(n₃₁*s₃₂+n₃₂*s₃₁)/𝐿₃²
                ms₂₂ += (q₁*s₃₁+q₂*s₃₂)*n₃₂*s₃₂/𝐿₃²
                ∂ms₂₂∂x += (∂q₁∂x*s₃₁+∂q₂∂x*s₃₂)*n₃₂*s₃₂/𝐿₃²
                ∂ms₂₂∂y += (∂q₁∂y*s₃₁+∂q₂∂y*s₃₂)*n₃₂*s₃₂/𝐿₃²
                if ξ.ξ == 0.0
                    Δms₁₁ =  𝒒̂ᵀ𝗚⁻¹𝒒*(n₃₁*s₃₁/𝐿₃²-n₁₁*s₁₁/𝐿₃²)
                    Δ∂ms₁₁∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*(n₃₁*s₃₁/𝐿₃²-n₁₁*s₁₁/𝐿₃²)
                    Δ∂ms₁₁∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*(n₃₁*s₃₁/𝐿₃²-n₁₁*s₁₁/𝐿₃²)
                    Δms₂₂ =  𝒒̂ᵀ𝗚⁻¹𝒒*(n₃₂*s₃₂/𝐿₃²-n₁₂*s₁₂/𝐿₁²)
                    Δ∂ms₂₂∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*(n₃₂*s₃₂/𝐿₃²-n₁₂*s₁₂/𝐿₁²)
                    Δ∂ms₂₂∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*(n₃₂*s₃₂/𝐿₃²-n₁₂*s₁₂/𝐿₁²)
                    Δms₁₂ =  𝒒̂ᵀ𝗚⁻¹𝒒*0.5*((n₃₁*s₃₂+n₃₂*s₃₁)/𝐿₃²-(n₁₁*s₁₂+n₁₂*s₁₁)/𝐿₁²)
                    Δ∂ms₁₂∂x = ∂𝒒̂ᵀ∂x𝗚⁻¹𝒒*0.5*((n₃₁*s₃₂+n₃₂*s₃₁)/𝐿₃²-(n₁₁*s₁₂+n₁₂*s₁₁)/𝐿₁²)
                    Δ∂ms₁₂∂y = ∂𝒒̂ᵀ∂y𝗚⁻¹𝒒*0.5*((n₃₁*s₃₂+n₃₂*s₃₁)/𝐿₃²-(n₁₁*s₁₂+n₁₂*s₁₁)/𝐿₁²)
                end
            end

            W₁₁₁ = mn₁₁n₁*wᵇ
            ∂W₁₁₁∂x = ∂mn₁₁∂xn₁*wᵇ
            ∂W₁₁₁∂y = ∂mn₁₁∂yn₁*wᵇ
            W₁₁₂ = mn₁₁n₂*wᵇ
            ∂W₁₁₂∂x = ∂mn₁₁∂xn₂*wᵇ
            ∂W₁₁₂∂y = ∂mn₁₁∂yn₂*wᵇ
            W₁₂₁ = mn₁₂n₁*wᵇ
            ∂W₁₂₁∂x = ∂mn₁₂∂xn₁*wᵇ
            ∂W₁₂₁∂y = ∂mn₁₂∂yn₁*wᵇ
            W₁₂₂ = mn₁₂n₂*wᵇ
            ∂W₁₂₂∂x = ∂mn₁₂∂xn₂*wᵇ
            ∂W₁₂₂∂y = ∂mn₁₂∂yn₂*wᵇ
            W₂₂₁ = mn₂₂n₁*wᵇ
            ∂W₂₂₁∂x = ∂mn₂₂∂xn₁*wᵇ
            ∂W₂₂₁∂y = ∂mn₂₂∂yn₁*wᵇ
            W₂₂₂ = mn₂₂n₂*wᵇ
            ∂W₂₂₂∂x = ∂mn₂₂∂xn₂*wᵇ
            ∂W₂₂₂∂y = ∂mn₂₂∂yn₂*wᵇ
            W₁₁ = (q₁₁*w + 2*(q₁n₁+ms₁₁)*wᵇ)/4/𝐴 + Δms₁₁
            ∂W₁₁∂x = (∂q₁₁∂x*w + 2*(∂q₁∂xn₁+∂ms₁₁∂x)*wᵇ)/4/𝐴 + Δ∂ms₁₁∂x
            ∂W₁₁∂y = (∂q₁₁∂y*w + 2*(∂q₁∂yn₁+∂ms₁₁∂y)*wᵇ)/4/𝐴 + Δ∂ms₁₁∂y
            W₁₂ = (q₁₂*w + (q₁n₂+q₂n₁+2*ms₁₂)*wᵇ)/4/𝐴 + Δms₁₂
            ∂W₁₂∂x = (∂q₁₂∂x*w + (∂q₁∂xn₂+∂q₂∂xn₁+2*∂ms₁₂∂x)*wᵇ)/4/𝐴 + Δ∂ms₁₂∂x
            ∂W₁₂∂y = (∂q₁₂∂y*w + (∂q₁∂yn₂+∂q₂∂yn₁+2*∂ms₁₂∂y)*wᵇ)/4/𝐴 + Δ∂ms₁₂∂y
            W₂₂ = (q₂₂*w + 2*(q₂n₂+ms₂₂)*wᵇ)/4/𝐴 + Δms₂₂
            ∂W₂₂∂x = (∂q₂₂∂x*w + 2*(∂q₂∂xn₂+∂ms₂₂∂x)*wᵇ)/4/𝐴 + Δ∂ms₂₂∂x
            ∂W₂₂∂y = (∂q₂₂∂y*w + 2*(∂q₂∂yn₂+∂ms₂₂∂y)*wᵇ)/4/𝐴 + Δ∂ms₂₂∂y
            for i in 1:length(𝓒)
                ∂²𝝭∂x²[i] += 𝝭[i]*W₁₁ + ∂𝝭∂x[i]*W₁₁₁ + ∂𝝭∂y[i]*W₁₁₂
                ∂∂²𝝭∂x²∂x[i] += 𝝭[i]*∂W₁₁∂x + ∂𝝭∂x[i]*∂W₁₁₁∂x + ∂𝝭∂y[i]*∂W₁₁₂∂x
                ∂∂²𝝭∂x²∂y[i] += 𝝭[i]*∂W₁₁∂y + ∂𝝭∂x[i]*∂W₁₁₁∂y + ∂𝝭∂y[i]*∂W₁₁₂∂y
                ∂²𝝭∂x∂y[i] += 𝝭[i]*W₁₂ + ∂𝝭∂x[i]*W₁₂₁ + ∂𝝭∂y[i]*W₁₂₂
                ∂∂²𝝭∂x∂y∂x[i] += 𝝭[i]*∂W₁₂∂x + ∂𝝭∂x[i]*∂W₁₂₁∂x + ∂𝝭∂y[i]*∂W₁₂₂∂x
                ∂∂²𝝭∂x∂y∂y[i] += 𝝭[i]*∂W₁₂∂y + ∂𝝭∂x[i]*∂W₁₂₁∂y + ∂𝝭∂y[i]*∂W₁₂₂∂y
                ∂²𝝭∂y²[i] += 𝝭[i]*W₂₂ + ∂𝝭∂x[i]*W₂₂₁ + ∂𝝭∂y[i]*W₂₂₂
                ∂∂²𝝭∂y²∂x[i] += 𝝭[i]*∂W₂₂∂x + ∂𝝭∂x[i]*∂W₂₂₁∂x + ∂𝝭∂y[i]*∂W₂₂₂∂x
                ∂∂²𝝭∂y²∂y[i] += 𝝭[i]*∂W₂₂∂y + ∂𝝭∂x[i]*∂W₂₂₁∂y + ∂𝝭∂y[i]*∂W₂₂₂∂y
            end
        end
        for i in 1:length(𝓒)
            ξ̂.𝝭[:∂x²][ξ̂.index[ξ̂.id]+i] = ∂²𝝭∂x²[i]
            ξ̂.𝝭[:∂x²∂x][ξ̂.index[ξ̂.id]+i] = ∂∂²𝝭∂x²∂x[i]
            ξ̂.𝝭[:∂x²∂y][ξ̂.index[ξ̂.id]+i] = ∂∂²𝝭∂x²∂y[i]
            ξ̂.𝝭[:∂x∂y][ξ̂.index[ξ̂.id]+i] = ∂²𝝭∂x∂y[i]
            ξ̂.𝝭[:∂x∂y∂x][ξ̂.index[ξ̂.id]+i] = ∂∂²𝝭∂x∂y∂x[i]
            ξ̂.𝝭[:∂x∂y∂y][ξ̂.index[ξ̂.id]+i] = ∂∂²𝝭∂x∂y∂y[i]
            ξ̂.𝝭[:∂y²][ξ̂.index[ξ̂.id]+i] = ∂²𝝭∂y²[i]
            ξ̂.𝝭[:∂y²∂x][ξ̂.index[ξ̂.id]+i] = ∂∂²𝝭∂y²∂x[i]
            ξ̂.𝝭[:∂y²∂y][ξ̂.index[ξ̂.id]+i] = ∂∂²𝝭∂y²∂y[i]
        end
    end
end

function set∇̄𝝭!(aps::Vector{T}) where T<:AbstractElement
    set𝝭!(aps)
    for ap in aps
        set∇̄𝝭!(ap)
    end
end

function set∇̄𝝭!(ap::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Seg2}) where {𝒑,𝑠,𝜙}
    𝓒 = ap.𝓒
    𝓖 = ap.𝓖
    for ξ̂ in 𝓖
        𝒒̂ = get𝒑₁(ap,ξ̂)
        𝗚⁻¹ = cal𝗚!(ap)
        𝒒̂ᵀ𝗚⁻¹ = 𝒒̂*𝗚⁻¹
        ∂𝝭∂x = ap.𝝭[:∂x]
        fill!(∂𝝭∂x,0.0)
        for ξ in ap.𝓖
            w = ξ.w
            n = ξ.n₁
            𝝭 = get𝝭(ap,ξ)
            𝒒 = get𝒑₁(ap,ξ)
            W₁ = 𝒒̂ᵀ𝗚⁻¹*𝒒*n*w
            for i in 1:length(𝓒)
                ∂𝝭∂x[i] += 𝝭[i]*W₁
            end
        end
        for i in 1:length(𝓒)
            ξ̂.𝝭[:∂̄x][ξ̂.index[ξ̂.id]+i] = ∂𝝭∂x[i]
        end
    end
end

function set∇̄𝝭!(ap::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Tri3}) where {𝒑,𝑠,𝜙}
    𝓒 = ap.𝓒
    𝓖 = ap.𝓖
    for ξ̂ in 𝓖
        𝒒̂ = get𝒑₁(ap,ξ̂)
        𝗚⁻¹ = cal𝗚!(ap)
        𝒒̂ᵀ𝗚⁻¹ = 𝒒̂*𝗚⁻¹
        ∂𝝭∂x = ap.𝝭[:∂x]
        ∂𝝭∂y = ap.𝝭[:∂y]
        fill!(∂𝝭∂x,0.0)
        fill!(∂𝝭∂y,0.0)
        for ξ in ap.𝓖
            w = ξ.w
            n₁ = ξ.n₁
            n₂ = ξ.n₂
            𝝭 = get𝝭(ap,ξ)
            𝒒 = get𝒑₁(ap,ξ)
            𝒒̂ᵀ𝗚⁻¹𝒒 = 𝒒̂ᵀ𝗚⁻¹*𝒒
            W₁ = 𝒒̂ᵀ𝗚⁻¹𝒒*n₁*w
            W₂ = 𝒒̂ᵀ𝗚⁻¹𝒒*n₂*w
            for i in 1:length(𝓒)
                ∂𝝭∂x[i] += 𝝭[i]*W₁
                ∂𝝭∂y[i] += 𝝭[i]*W₂
            end
        end
        for i in 1:length(𝓒)
            ξ̂.𝝭[:∂̄x][ξ̂.index[ξ̂.id]+i] = ∂𝝭∂x[i]
            ξ̂.𝝭[:∂̄y][ξ̂.index[ξ̂.id]+i] = ∂𝝭∂y[i]
        end
    end
end

function set∇̄²𝝭!(ap::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Tri3}) where {𝒑,𝑠,𝜙}
    𝓒 = ap.𝓒
    𝓖 = ap.𝓖
    for ξ̂ in 𝓖
        𝒒̂ = get𝒑₂(ap,ξ̂)
        𝗚⁻¹ = cal𝗚₂!(ap)
        𝒒̂ᵀ𝗚⁻¹ = 𝒒̂*𝗚⁻¹
        M₁₁ = ap.𝝭[:M₁₁]
        M₁₂ = ap.𝝭[:M₁₂]
        M₂₂ = ap.𝝭[:M₂₂]
        fill!(M₁₁,0.0)
        fill!(M₁₂,0.0)
        fill!(M₂₂,0.0)
        for ξ in 𝓖
            w = ξ.w
            n₁ = ξ.n₁
            n₂ = ξ.n₂
            𝐿² = ξ.𝐿^2
            _,∂𝝭∂x,∂𝝭∂y = get∇𝝭(ap,ξ)
            𝒒 = get𝒑₂(ap,ξ)
            𝒒̂ᵀ𝗚⁻¹𝒒 = 𝒒̂ᵀ𝗚⁻¹*𝒒
            W₁₁₁ = 𝒒̂ᵀ𝗚⁻¹𝒒*n₁*n₁*n₁*w/𝐿²
            W₁₁₂ = 𝒒̂ᵀ𝗚⁻¹𝒒*n₁*n₁*n₂*w/𝐿²
            W₂₂₁ = 𝒒̂ᵀ𝗚⁻¹𝒒*n₂*n₂*n₁*w/𝐿²
            W₂₂₂ = 𝒒̂ᵀ𝗚⁻¹𝒒*n₂*n₂*n₂*w/𝐿²
            W₁₂₁ = 𝒒̂ᵀ𝗚⁻¹𝒒*n₁*n₂*n₁*w/𝐿²
            W₁₂₂ = 𝒒̂ᵀ𝗚⁻¹𝒒*n₁*n₂*n₂*w/𝐿²
            for i in 1:length(𝓒)
                M₁₁[i] += ∂𝝭∂x[i]*W₁₁₁ + ∂𝝭∂y[i]*W₁₁₂
                M₁₂[i] += ∂𝝭∂x[i]*W₁₂₁ + ∂𝝭∂y[i]*W₁₂₂
                M₂₂[i] += ∂𝝭∂x[i]*W₂₂₁ + ∂𝝭∂y[i]*W₂₂₂
            end
        end
        n₁ = ξ̂.n₁
        n₂ = ξ̂.n₂
        for i in 1:length(𝓒)
            ξ̂.𝝭[:Mₙₙ][ξ̂.index[ξ̂.id]+i] = M₁₁[i]*n₁*n₁ + 2*M₁₂[i]*n₁*n₂ + M₂₂[i]*n₂*n₂
        end
    end
end

function setV̄ₙ!(ap::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Tri3}) where {𝒑,𝑠,𝜙}
    𝓒 = ap.𝓒
    𝓖 = ap.𝓖
    x₁ = 𝓒[1].x;y₁ = 𝓒[1].y
    x₂ = 𝓒[2].x;y₂ = 𝓒[2].y
    x₃ = 𝓒[3].x;y₃ = 𝓒[3].y
    𝐴 = get𝐴(gp)
    n₁₁ = y₃-y₂;n₂₁ = y₁-y₃;n₃₁ = y₂-y₁
    n₁₂ = x₂-x₃;n₂₂ = x₃-x₁;n₃₂ = x₁-x₂

    for ξ̂ in 𝓖
        _,∂𝒒̂∂ξ,∂𝒒̂∂η = get∇𝒑₂(ap,ξ̂)
        𝗚⁻¹ = cal𝗚₂!(ap)
        ∂𝒒̂ᵀ∂ξ𝗚⁻¹ = ∂𝒒̂∂ξ*𝗚⁻¹
        ∂𝒒̂ᵀ∂η𝗚⁻¹ = ∂𝒒̂∂η*𝗚⁻¹
        V₁₁₁ = ap.𝝭[:V₁₁₁]
        V₁₁₂ = ap.𝝭[:V₁₁₂]
        V₁₂₁ = ap.𝝭[:V₁₂₁]
        V₁₂₂ = ap.𝝭[:V₁₂₂]
        V₂₁₁ = ap.𝝭[:V₂₁₁]
        V₂₁₂ = ap.𝝭[:V₂₁₂]
        V₂₂₁ = ap.𝝭[:V₂₂₁]
        V₂₂₂ = ap.𝝭[:V₂₂₂]
        fill!(V₁₁₁,0.0)
        fill!(V₁₁₂,0.0)
        fill!(V₁₂₁,0.0)
        fill!(V₁₂₂,0.0)
        fill!(V₂₁₁,0.0)
        fill!(V₂₁₂,0.0)
        fill!(V₂₂₁,0.0)
        fill!(V₂₂₂,0.0)
        for ξ in ap.𝓖
            w = ξ.w
            n₁ = ξ.n₁
            n₂ = ξ.n₂
            s₁ = -ξ.n₂
            s₂ = ξ.n₁
            𝐿² = ξ.𝐿^2
            𝝭 = get𝝭(ap,ξ)
            _,∂𝒒∂ξ,∂𝒒∂η = get∇𝒑₂(ap,ξ)
            ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂ξ = ∂𝒒̂ᵀ∂ξ𝗚⁻¹*∂𝒒∂ξ
            ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂η = ∂𝒒̂ᵀ∂ξ𝗚⁻¹*∂𝒒∂η
            ∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂ξ = ∂𝒒̂ᵀ∂η𝗚⁻¹*∂𝒒∂ξ
            ∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂η = ∂𝒒̂ᵀ∂η𝗚⁻¹*∂𝒒∂η
            ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂x = ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂ξ*n₁₁*n₁₁ + (∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂η + ∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂ξ)*n₁₁*n₂₁ + ∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂η*n₂₁*n₂₁
            ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂y = ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂ξ*n₁₁*n₁₂ + ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂η*n₁₁*n₂₂ + ∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂ξ*n₁₂*n₂₁ + ∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂η*n₂₁*n₂₂
            ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂x = ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂ξ*n₁₁*n₁₂ + ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂η*n₁₂*n₂₁ + ∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂ξ*n₁₁*n₂₂ + ∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂η*n₂₁*n₂₂
            ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂y = ∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂ξ*n₁₂*n₁₂ + (∂𝒒̂ᵀ∂ξ𝗚⁻¹∂𝒒∂η + ∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂ξ)*n₁₂*n₂₂ + ∂𝒒̂ᵀ∂η𝗚⁻¹∂𝒒∂η*n₂₂*n₂₂
            W₁₁₁ = ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂x*n₁ + (∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂x*s₁ + ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂y*s₂)*s₁*n₁*w/𝐿²
            W₁₁₂ = ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂x*n₁ + (∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂x*s₁ + ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂y*s₂)*s₁*n₁*w/𝐿²
            W₁₂₁ = ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂x*n₂ + (∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂x*s₁ + ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂y*s₂)*s₁*n₂*w/𝐿²
            W₁₂₂ = ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂x*n₂ + (∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂x*s₁ + ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂y*s₂)*s₁*n₂*w/𝐿²
            W₂₁₁ = ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂y*n₁ + (∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂x*s₁ + ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂y*s₂)*s₂*n₁*w/𝐿²
            W₂₁₂ = ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂y*n₁ + (∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂x*s₁ + ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂y*s₂)*s₂*n₁*w/𝐿²
            W₂₂₁ = ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂y*n₂ + (∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂x*s₁ + ∂𝒒̂ᵀ∂x𝗚⁻¹∂𝒒∂y*s₂)*s₂*n₂*w/𝐿²
            W₂₂₂ = ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂y*n₂ + (∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂x*s₁ + ∂𝒒̂ᵀ∂y𝗚⁻¹∂𝒒∂y*s₂)*s₂*n₂*w/𝐿²
            for i in 1:length(𝓒)
                V₁₁₁[i] += 𝝭[i]*W₁₁₁
                V₁₁₂[i] += 𝝭[i]*W₁₁₂
                V₁₂₁[i] += 𝝭[i]*W₁₂₁
                V₁₂₂[i] += 𝝭[i]*W₁₂₂
                V₂₁₁[i] += 𝝭[i]*W₂₁₁
                V₂₁₂[i] += 𝝭[i]*W₂₁₂
                V₂₂₁[i] += 𝝭[i]*W₂₂₁
                V₂₂₂[i] += 𝝭[i]*W₂₂₂
            end
        end
        n₁ = ξ̂.n₁
        n₂ = ξ̂.n₂
        s₁ = -ξ̂.n₂
        s₂ = ξ̂.n₁
        for i in 1:length(𝓒)
            ξ̂.𝝭[:Vₙ][ξ̂.index[ξ̂.id]+i] = V₁₁₁[i]*n₁ + V₁₂₂[i]*n₁ + V₂₁₁[i]*n₂ + V₂₂₂[i]*n₂ + (V₁₁₁*s₁+V₁₁₂*s₂)*s₁*n₁ + (V₁₂₁*s₁+V₁₂₂*s₂)*s₁*n₂ + (V₂₁₁*s₁+V₂₁₂*s₂)*s₂*n₁ + (V₂₂₁*s₁+V₂₂₂*s₂)*s₂*n₂
        end
    end
end


text:
using SafeTestsets

@safetestset "Test Uniform Mesh Construction" begin
    include("test_uniform_mesh.jl")
end

@safetestset "Test Mesh Attributes" begin
    include("test_mesh_attributes.jl")
end


text:
"""

## sparin!

Function sparin! performs Cholesky factorisation on a symmetric
skyline global matrix. The vector kv is updated.

###Arguments
```julia
sparin!(kv, kdiag)
```

### Arguments
```julia
* kv::Vector{Float64}       : Global stiffness matrix (Updated)
* kdiag::Vector{Int}      : Diagonal elemnt index vector
```
"""
function sparin!(kv::Vector{Float64}, kdiag::Vector{Int})
  local x::Float64
  n = size(kdiag, 1)
  kv[1] = sqrt(kv[1])
  for i in 2:n
    ki = kdiag[i] - i
    l = kdiag[i-1] - ki + 1
    for j in l:i
      x = kv[ki+j]
      kj = kdiag[j] - j
      if j !== 1
        ll = kdiag[j-1] - kj + 1
        ll = max(l, ll)
        if ll !== j
          m = j - 1
          for k in ll:m
            x = x - kv[ki+k]*kv[kj+k]
          end
        end
      end
      kv[ki+j] = x/kv[kj+j]
    end
    kv[ki+i] = sqrt(x)
  end
end


text:
function get_stream(connection::HTTPConnection, stream_identifier::UInt32)
    @assert stream_identifier != 0x0

    for stream in connection.streams
        if stream.stream_identifier == stream_identifier
            return stream
        end
    end

    stream = HTTPStream(stream_identifier, IDLE,
                        connection.settings.initial_window_size, nothing)

    push!(connection.streams, stream)
    return stream
end

function stream_states(connection::HTTPConnection)
    result = Array{Tuple{UInt32, STREAM_STATE}, 1}()

    for i = 1:length(connection.streams)
        push!(result, (connection.streams[i].stream_identifier, connection.streams[i].state))
    end

    return result
end

function get_dependency_parent(connection::HTTPConnection, stream_identifier::UInt32)
    stream = get_stream(connection, stream_identifier)

    if isnothing(stream.priority)
        return nothing
    else
        return get_stream(stream.priority.dependent_stream_identifier)
    end
end

function get_dependency_children(connection::HTTPConnection, stream_identifier::UInt32)
    result = Array{Stream}()

    for i = 1:length(connection.streams)
        stream = connection.streams[i]

        if !isnothing(stream.priority) &&
            stream.priority.dependent_stream_identifier == stream_identifier
            push!(result, stream)
        end
    end

    return result
end

function handle_priority!(connection::HTTPConnection, stream_identifier::UInt32,
                          exclusive::Bool, dependent_stream_identifier::UInt32, weight::UInt8)
    stream = get_stream(connection, stream_identifier)

    if exclusive
        children = get_dependency_children(connection, stream_identifier)

        for i = 1:length(children)
            children[i].priority.dependent_stream_identifier = stream_identifier
        end
    end

    stream.priority = Priority(dependent_stream_identifier, weight)
end

function concurrent_streams_count(connection::HTTPConnection)
    n = 0

    for stream in connection.streams
        if stream.state != IDLE && stream.state != CLOSED
            n += 1
        end
    end

    return n
end