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

"""
    Identity()

The identity transform `x -> x`.
"""
struct Identity <: Stateless end

isrevertible(::Type{Identity}) = true

apply(::Identity, table) = table, nothing

revert(::Identity, newtable, cache) = newtable


text:
using DataFrames

# dataframes:
#
function test_read_df()
    df = dataframe_read("../data/iris/iris150.csv")
    return nrow(df) == 150
end

function test_df_loader()

        trn = DataFrame(x1=randn(16), x2=randn(16),
                        x3=randn(16), x4=randn(16),
                        x5=randn(16), x6=randn(16),
                        x7=randn(16), x8=randn(16),
                        y=["blue", "red", "green", "green",
                           "blue", "red", "green", "green",
                           "blue", "red", "green", "green",
                           "blue", "red", "green", "green"])

        mb1 = dataframe_minibatches(trn, size=4, teaching="y", ignore="x1")

        trn.iy = [1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4]
        @show mb2 = dataframe_minibatches(trn, size=4, teaching="iy", ignore=["x1", "y"])
        @show first(mb2)[2]

        trn.is = [:a,:a,:a,:a,:b,:b,:b,:b,:c,:c,:c,:c,:d,:d,:d,:d]
        mb3 = dataframe_minibatches(trn, size=4, teaching="is", ignore=["iy", "y"])

        return first(mb1)[2] == UInt8[0x01  0x03  0x02  0x02] &&
               first(mb2)[2] == UInt8[0x01 0x01 0x01 0x01]    &&
               mb3 == nothing
end

function test_df_split()
    df = dataframe_read("../data/iris/iris150.csv")
    t,v = dataframe_split(df, fr=0.5, shuffle=true,
          teaching="species", balanced=true)
    return nrow(v) == 75
end

function test_df_class_ids()
    df = dataframe_read("../data/iris/iris150.csv")
    c = mk_class_ids(df.species)
    return c[2] == ["setosa", "versicolor", "virginica"]
end

function test_df_minibatch()
    df = dataframe_read("../data/iris/iris150.csv")
    mb = dataframe_minibatches(df, size=10, teaching="species")
    return size(first(mb)[1]) == (4,10)
end

function test_df_errors()
    df1 = dataframe_read("../data/iris/iris150.ods")
    df2 = dataframe_read("../data/iris/iris150.dat")

    return df1 == nothing &&
           df2 == nothing
end

# test NLP utils:
#
function test_tokenizer()
    tok = WordTokenizer(["I love Julia",
                         "Peter loves Python",
                         "We all marvel Geoff"])
    l = tok(["I love Julia", "Peter loves Python", "We all marvel Geoff"],
            add_ctls=true)

    sentence = tok([8, 12, 5])
    return tok("Julia") == 8 && tok(8) == "Julia" &&
           l[3] == [1, 6, 14, 13, 11, 2] &&
           sentence == "Julia loves Peter"
end

function test_seq_mb()
    tok = WordTokenizer(["I love Julia",
                         "Peter loves Python",
                         "We all marvel Geoff"])

    t = [tok("I love Julia", split_words=true),
         tok("Peter loves Python", split_words=true),
         tok("Peter loves Julia and Scala", split_words=true)]

    mb = seq_minibatch(t, 2, seq_len=4)
    return size(first(mb)) == (4,2)
end

function test_seq_mb_xy()
    tok = WordTokenizer(["I love Julia",
                         "Peter loves Python",
                         "We all marvel Geoff"])

    t = [tok("I love Julia", split_words=true),
         tok("Peter loves Python", split_words=true),
         tok("Peter loves Julia and Scala", split_words=true)]

    y = [1 1 2]
    mb = seq_minibatch(t, y, 2, seq_len=4)
    return first(mb)[2] == [1 1]
end

# test tatoeba on tiny dataset Galician:
#
function test_tatoeba()
    en,glg = get_tatoeba_corpus("glg")
    return en isa AbstractArray
end

function test_seq2seq_mb()
     en,glg = get_tatoeba_corpus("glg")
     ven = WordTokenizer(en)
     vglg = WordTokenizer(glg)
     mb = seq2seq_minibatch(ven.(en, split_words=true),
            vglg.(glg, split_words=true),
            32, seq_len=20,
            pad_x=ven("<pad>"), pad_y=vglg("<pad>"))

    mb1 = first(mb)[1]
    return size(mb1) == (20,32)
end

ext: jl


text:
module MLJOpenML

const OpenML = MLJOpenML
export OpenML

include("openml.jl")

end # module


text:
function make_linkstate_table(tlist,c,L)
    for linktype in tlist
        from=linktype.from
        to=linktype.to
        if isa(from,Array)
            for x in from
                if x<1 || x>c
                    msg="One of the LinkTypes passed to LargeNet contains refers to nodes of state $x,"
                    msg*=" but $x is not a valid node state in this network."
                    throw(ArgumentError(msg))
                end
            end
        else
            if from<0 || from>c
                msg="One of the LinkTypes passed to LargeNet contains refers to nodes of state $from,"
                msg*=" but $from is not a valid node state in this network."
                throw(ArgumentError(msg))
            end
        end
        if isa(to,Array)
            for x in to
                if x<1 || x>c
                    msg="One of the LinkTypes passed to LargeNet contains refers to nodes of state $x,"
                    msg*=" but $x is not a valid node state in this network."
                    throw(ArgumentError(msg))
                end
            end
        else
            if to<0 || to>c
                msg="One of the LinkTypes passed to LargeNet contains refers to nodes of state $to,"
                msg*=" but $to is not a valid node state in this network."
                throw(ArgumentError(msg))
            end
        end
    end

    m=Array{Int,2}(undef,c,c)
    for i=1:c
        for j=1:c
            m[i,j]=L-1
            for rule=1:L-2
                dir=tlist[rule].dir
                ft=tlist[rule].from
                tt=tlist[rule].to
                if dir!=2
                    continue
                end
                if !((i in ft && j in tt) || (j in ft && i in tt) || (ft==0 && j in tt) || (ft==0 && i in tt) || (i in ft && tt==0) ||  (ft==j && tt==0) || (ft==0 && tt==0) )
                    continue
                end
                if (m[i,j]!=L-1)
                    throw(ArgumentError("LinkTyoe number $rule, passed to the FastNet constructure, seems to say that links of the form from state $i to state $j are of type $rule, but they were previously defined as LinkType $(m[i,j])."))
                end
                m[i,j]=rule
            end
        end
    end
    for i=1:c
        for j=1:c
            for rule=1:L-2
                dir=tlist[rule].dir
                ft=tlist[rule].from
                tt=tlist[rule].to
                if dir!=1
                    continue
                end
                if !(j in tt) && tt!=0
                    continue
                end
                if !(i in ft) && ft!=0
                    continue
                end
                if (m[i,j]!=L-1)
                    throw(ArgumentError("LinkTyoe number $rule, passed to the FastNet constructure, seems to say that links of the form from state $i to state $j are of type $rule, but they were previously defined as LinkType $(m[i,j])."))
                end
                m[i,j]=rule
            end
        end
    end
    m
end

function verify_linktype!(from,to,dir)
    if dir!=1 && dir!=2
        throw(ArgumentError("Third argument of LinkType must be 1 (unidirectional) or 2 (bidirection)."))
    end
    if isa(from,Tuple)
        from=[from...]
    end
    if isa(to,Tuple)
        to=[to...]
    end

    if from=="*"
        from=0
    end
    if to=="*"
        to=0
    end

    if isa(from,Array)
        from=reshape(from,1,:)
        l=length(from)
        for i=1:l
            try
                from[i]=convert(Int,from[i])
            catch e
                msg="Could not process FROM argument passed to constructor of LinkType: "
                msg*="the variable passed seems to be a collection, "
                msg*="but I could not convert element $i to Int."
                throw(ArgumentError(msg))
            end
        end
    else
        try
            from=convert(Int,from)
        catch
            msg="Could not process FROM argument passed to constructor of LinkType: "
            msg*="the argument must be either an Integer, \"*\", or a Tuple or Array Integers."
            throw(ArgumentError(msg))
        end
    end

    if isa(to,Array)
        to=reshape(to,1,:)
        l=length(to)
        for i=1:l
            try
                to[i]=convert(Int,to[i])
            catch e
                msg="Could not process TO argument passed to constructor of LinkType: "
                msg*="the variable passed seems to be a collection, "
                msg*="but I could not convert element $i to Int."
                throw(ArgumentError(msg))
            end
        end
    else
        try
            tp=convert(Int,to)
        catch
            msg="Could not process TO argument passed to constructor of LinkType: "
            msg*="the argument must be either an Integer, \"*\", or a Tuple or Array of Integers."
            throw(ArgumentError(msg))
        end
    end
    (from,to,dir)
end


text:
# ---
# title: 740. Delete and Earn
# id: problem740
# author: Indigo
# date: 2021-05-14
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/delete-and-earn/description/>
# hidden: true
# ---
#
# Given an array `nums` of integers, you can perform operations on the array.
#
# In each operation, you pick any `nums[i]` and delete it to earn `nums[i]`
# points. After, you must delete **every** element equal to `nums[i] - 1` or
# `nums[i] + 1`.
#
# You start with 0 points. Return the maximum number of points you can earn by
# applying such operations.
#
# **Example 1:**
#
#
#
#     Input: nums = [3, 4, 2]
#     Output: 6
#     Explanation:
#     Delete 4 to earn 4 points, consequently 3 is also deleted.
#     Then, delete 2 to earn 2 points. 6 total points are earned.
#
#
#
#
# **Example 2:**
#
#
#
#     Input: nums = [2, 2, 3, 3, 3, 4]
#     Output: 9
#     Explanation:
#     Delete 3 to earn 3 points, deleting both 2's and the 4.
#     Then, delete 3 again to earn 3 points, and 3 again to earn 3 points.
#     9 total points are earned.
#
#
#
#
# **Note:**
#
#   * The length of `nums` is at most `20000`.
#   * Each element `nums[i]` is an integer in the range `[1, 10000]`.
#
#
#
#
## @lc code=start
using LeetCode

function delete_and_earn(nums::Vector{Int})
    if length(nums) == 0
        return 0
    elseif length(nums) == 1
        return nums[1]
    end
    cntr = fill(0, maximum(nums))
    for num in nums
        cntr[num] += num
    end
    fst, scd = cntr[1], max(cntr[2], cntr[1])
    for i in 3:length(cntr)
        fst, scd = scd, max(fst + cntr[i], scd)
    end
    return scd
end
## @lc code=end


text:
@testset "Term coupling" begin
    @testset "LS Coupling" begin
        # Coupling two 2S terms should yields singlets and triplets of S,P,D
        @test couple_terms(Term(1,1//2,1), Term(1,1//2,1)) ==
            sort([Term(2,0,1), Term(2,1,1), Term(0,0,1), Term(1,1,1), Term(1,0,1), Term(0,1,1)])

        @test couple_terms(Term(0,1//2,1), Term(1,1//2,-1)) ==
            sort([Term(1,0,-1), Term(1,1,-1)])

        @test couple_terms(Term(0,1//2,-1), Term(1,1//2,-1)) ==
            sort([Term(1,0,1), Term(1,1,1)])

        function test_coupling(o1, o2)
            c1 = couple_terms(terms(o1), terms(o2))
            r = o1 + o2
            c2 = terms(r)
            @test c1 == c2
        end

        test_coupling(c"1s", c"2p")
        test_coupling(c"2s", c"2p")
        test_coupling(c"2p", c"3p")
        test_coupling(c"2p", c"3d")

        @test terms(c"1s") == [T"2S"]
        @test terms(c"1s2") == [T"1S"]
        @test terms(c"2p") == [T"2Po"]
        @test terms(c"[He]") == [T"1S"]
        @test terms(c"[Ne]") == [T"1S"]
        @test terms(c"[Ar]") == [T"1S"]
        @test terms(c"[Kr]") == [T"1S"]
        @test terms(c"[Xe]") == [T"1S"]
        @test terms(c"[Rn]") == [T"1S"]

        @test terms(c"1s ks") == [T"1S", T"3S"]
        @test terms(c"1s kp") == [T"1Po", T"3Po"]
    end

    @testset "Coupling of jj terms" begin
        @test couple_terms(1, 2) isa Vector{Int}
        @test couple_terms(1//2, 1//2) isa Vector{HalfInt}
        @test couple_terms(1.0, 1.5) isa Vector{HalfInt}

        @test couple_terms(HalfInteger(1//2), HalfInteger(0)) == [1//2]
        @test couple_terms(1//2, HalfInteger(0)) == [1//2]
        @test couple_terms(HalfInteger(1//2), 0) == [1//2]

        @test couple_terms(0, 1) == 1:1
        @test couple_terms(1//2,1//2) == 0:1
        @test couple_terms(1, 1) == 0:2
        @test couple_terms(1//1, 3//2) == 1//2:5//2

        @test couple_terms(0, 1.0) == 1:1
        @test couple_terms(1.0, 1.5) == 1//2:5//2

        @test intermediate_couplings([1, 2], 3) isa Vector{Vector{Int}}
        @test intermediate_couplings([1, 2], HalfInteger(3)) isa Vector{Vector{HalfInt}}
        @test intermediate_couplings([1, 2], HalfInteger(3//2)) isa Vector{Vector{HalfInt}}
        @test intermediate_couplings([1, 2.5], 3) isa Vector{Vector{HalfInt}}

        @test intermediate_couplings([0, 0], 0) == [[0,0,0]]
        @test intermediate_couplings([0, 0], 1) == [[1,1,1]]
        @test intermediate_couplings([1//2, 1//2], 0) == [[0,1//2,0],[0,1//2,1]]

        @test terms(rc"1s 2s") == [0,1]
    end
end


text:
using LinearAlgebra, StatsBase, ProgressMeter
using TensorOperations

struct Index
    dim::Array{Int64}
    name::Array{String}
end

mutable struct Tensor
    data::Array{Float64}
    index::Index
end

struct MPS
    site::Array{Tensor}
    siteIndex::Array{Integer}
    MPS(A,range) = new(A,collect(range))
end

function generate_MPS(func::Function; D=6, spin_deg=2, mps_size=10)::MPS

    makeTensor(theIndex) = Tensor( func(Float64,tuple(theIndex.dim[:]...)),theIndex )

    mps = MPS(
        vcat(
            vcat(
                   [makeTensor(Index([1,D,spin_deg],["ileft","i1","s1"]))],
                   [makeTensor(Index([D,D,spin_deg],["i$(i-1)","i$(i)","s$i"])) for i in 2:mps_size-1]
            ),
            [makeTensor(Index([D,1,spin_deg],["i$(mps_size-1)","iright","s$(mps_size)"]))]
        ),
        1:mps_size
    )
#     mps = MPS(
#         vcat(
#             vcat(
#                    [makeTensor(Index([D,D,spin_deg],["ileft","i1","s1"]))],
#                    [makeTensor(Index([D,D,spin_deg],["i$(i-1)","i$(i)","s$i"])) for i in 2:mps_size-1]
#             ),
#             [makeTensor(Index([D,D,spin_deg],["i$(mps_size-1)","iright","s$(mps_size)"]))]
#         ),
#         1:mps_size
#     )
    return mps
end

function generate_Heisenberg_Hamiltonian()
    σx = [0  1;1  0]
    iσy = [0 1;-1  0]
    σz = [1  0;0  -1]
    self_kron(A::Matrix) = tensorproduct(A,(11,19),A,(21,29),(11,21,19,29))
    H = 0.25*(self_kron(σx)-self_kron(iσy)+self_kron(σz))
    return H # J=1
end

function generate_Heisenberg_Hamiltonian(Jx,Jy,Jz)
    σx = [0  1;1  0]
    iσy = [0 1;-1  0]
    σz = [1  0;0  -1]
    self_kron(A::Matrix) = tensorproduct(A,(11,19),A,(21,29),(11,21,19,29))
    # bug??
    H = 0.25*(Jx*self_kron(σx)-Jy*self_kron(iσy)+Jz*self_kron(σz))
    return H
end

function generate_eye()
    eye = [1 0;0 1]
    self_kron(A::Matrix) = tensorproduct(A,(11,19),A,(21,29),(11,21,19,29))
    eye2 = 1.0*self_kron(eye)
    return eye2
end

function sparse_indices(H)::Dict
    C1 = CartesianIndices(H[:,:,1,1])
    C2 = CartesianIndices(H[1,1,:,:])
    D = Dict{Tuple{Int8,Int8},Array{Tuple{Int8,Int8,Float64},1}}()
    for i in C1
        tmp = Tuple{Int8,Int8,Float64}[]
        # adding diagonal index(j=i) first
        push!(tmp,(i.I...,H[i,i]))
        for j in C2
            if H[i,j] != 0.0 && i != j
                push!(tmp,(j.I...,H[i,j]))
            end
        end
        push!(D, i.I => tmp)
    end
    return D
end

# right regularization:
# BBB  B
# BBB RA
# BBB' A
# BBRA A
# BB' AA
function right_regularize!(mps::MPS)
    arrS = mps.site
    len = length(arrS)
    for i in len:-1:2
        A = arrS[i].data
        B = arrS[i-1].data
        l, r, s = arrS[i].index.dim
        F = lq( reshape(A,l,r*s) )
        if l <= r*s
            A[:] = reshape(Matrix(F.Q),l,r,s)
            # the "l r s" below are just simbols in the macro, they are not what we defined above.
            @tensor C[l,r',s] := B[l,r,s]*F.L[r,r']
        else
            A[:] = reshape(Matrix(I,l,r*s)*F.Q,l,r,s)
            L = 0.0*Matrix(I,l,l)
            L[:,1:r*s] = F.L
            @tensor C[l,r',s] := B[l,r,s]*L[r,r']
        end
        B[:] = C
    end
    arrS[1].data[:] = arrS[1].data / norm(arrS[1].data,2)
#     normalize!(arrS[1].data)
end

function exact_energy!(mps::MPS)::Real
    right_regularize!(mps)
    S = mps.site
    Et = 0.0
    A = S[1].data
    B = S[2].data
    # all real numbers in A, B, H.
    @tensoropt Ei[:] := A[l,Lu,s11]*B[Lu,r,s21]*H[s11,s21,s12,s22]*A[l,Ld,s12]*B[Ld,r,s22]
    Et += scalar(Ei)
    for i in 2:length(S)-1
        P = S[i-1].data
        A = S[i].data
        B = S[i+1].data
        l, r, s = S[i-1].index.dim
        M = tensorcopy(P,(1,2,3),(1,3,2))
        F = qr( reshape(M,l*s,r) )
        if l*s < r
            P2 = reshape( F.Q*Matrix(I,l*s,r), l,s,r )
            R = 0.0*Matrix(I,r,r)
            R[1:l*s,:] = F.R
            @tensor C[l',r,s] := R[l',l]*A[l,r,s]
        else
            P2 = reshape( Matrix(F.Q),l,s,r )
            R = F.R
            @tensor C[l',r,s] := R[l',l]*A[l,r,s]
        end
        tensorcopy!(P2,(1,3,2),P,(1,2,3))
        A[:] = C
        @tensoropt Ei[:] := A[l,Lu,s11]*B[Lu,r,s21]*H[s11,s21,s12,s22]*A[l,Ld,s12]*B[Ld,r,s22]
        Et += scalar(Ei)
    end
    return Et
end

function exact_energy(mps1::MPS)::Real
    mps = deepcopy(mps1)
    right_regularize!(mps)
    S = mps.site
    Et = 0.0
    A = S[1].data
    B = S[2].data
    # all real numbers in A, B, H.
    @tensoropt Ei[:] := A[l,Lu,s11]*B[Lu,r,s21]*H[s11,s21,s12,s22]*A[l,Ld,s12]*B[Ld,r,s22]
    Et += scalar(Ei)
    for i in 2:length(S)-1
        P = S[i-1].data
        A = S[i].data
        B = S[i+1].data
        l, r, s = S[i-1].index.dim
        M = tensorcopy(P,(1,2,3),(1,3,2))
        F = qr( reshape(M,l*s,r) )
        if l*s < r
            P2 = reshape( F.Q*Matrix(I,l*s,r), l,s,r )
            R = 0.0*Matrix(I,r,r)
            R[1:l*s,:] = F.R
            @tensor C[l',r,s] := R[l',l]*A[l,r,s]
        else
            P2 = reshape( Matrix(F.Q),l,s,r )
            R = F.R
            @tensor C[l',r,s] := R[l',l]*A[l,r,s]
        end
        tensorcopy!(P2,(1,3,2),P,(1,2,3))
        A[:] = C
        @tensoropt Ei[:] := A[l,Lu,s11]*B[Lu,r,s21]*H[s11,s21,s12,s22]*A[l,Ld,s12]*B[Ld,r,s22]
        Et += scalar(Ei)
    end
    return Et
end

function exact_Sz(mps1::MPS)::Real
    mps = deepcopy(mps1)
    right_regularize!(mps)
    S = mps.site
    len = length(S)
    Szt = 0.0
    A = S[1].data
#     B = S[2].data
    Ŝz = Matrix(0.5*Diagonal( (spin_deg-1):-2:-(spin_deg-1) ) )
    # all real numbers in A, B, H.
    @tensoropt Sz[:] := A[l,r,s1]*Ŝz[s1,s2]*A[l,r,s2]
    Szt += scalar(Sz)
    for i in 2:len
        P = S[i-1].data
        A = S[i].data
#         B = S[i+1].data
        l, r, s = S[i-1].index.dim
        M = tensorcopy(P,(1,2,3),(1,3,2))
        F = qr( reshape(M,l*s,r) )
        if l*s < r
            P2 = reshape( F.Q*Matrix(I,l*s,r), l,s,r )
            R = 0.0*Matrix(I,r,r)
            R[1:l*s,:] = F.R
            @tensor C[l',r,s] := R[l',l]*A[l,r,s]
        else
            P2 = reshape( Matrix(F.Q),l,s,r )
            R = F.R
            @tensor C[l',r,s] := R[l',l]*A[l,r,s]
        end
        tensorcopy!(P2,(1,3,2),P,(1,2,3))
        A[:] = C
        @tensoropt Sz[:] := A[l,r,s1]*Ŝz[s1,s2]*A[l,r,s2]
        Szt += scalar(Sz)
    end
    return Szt/len
end

function print_mat(A::AbstractArray)
    show(stdout,"text/plain",A)
    println(stdout,"\n")
end

# H = generate_Heisenberg_Hamiltonian()
# spH = sparse_indices(H)

mutable struct Sample
    state::Vector{Int8}
    inner::Float64   #W(S)
    energy::Float64  #E(S)
end

mutable struct Node
    spin::Int8
    value::Float64
    parent::Int32
end

function generate_BTree(s::Int)::Vector{Node}
    l = 2^(s+1)-2
    W = Vector{Node}(undef,l)
    for i in 1:2:l-1
        W[i] = Node(1,0.0,(i-1)/2)
    end
    for i in 2:2:l
        W[i] = Node(2,0.0,i/2-1)
    end
    return W
end

# configure next sample state in Markov chain
function next(tmp::Vector{Int8})
    l = length(tmp)
    new = deepcopy(tmp)
    for i in 1:2:l-1
        if tmp[i] != tmp[i+1]
            new[i] = tmp[i+1]
            new[i+1] = tmp[i]
        end
    end
    selected = StatsBase.sample(1:l,Int(ceil(l/4)),replace=false,ordered=true)
    new[selected] = 3 .- new[selected]
    return new
end

function state_inner(tmp::Sample,mps::MPS)
    l = length(tmp.state)
    if length(mps.site) != l
        throw(DomainError( (l,length(mps.site)), "length dismatch" ))
    end
    matchain = Array{Matrix{Float64},1}(undef,l)
    for i in 1:l
        spin = tmp.state[i]
        matchain[i] = mps.site[i].data[:,:,spin]
    end
    inner = matchain[1]
    for i in 2:l
        inner = inner*matchain[i]
    end
    return inner[1,1]
end

function state_inner(state::Vector{Int8},mps::MPS)
    l = length(state)
    if length(mps.site) != l
        throw(DomainError( (l,length(mps.site)), "length dismatch" ))
    end
    matchain = Array{Matrix{Float64},1}(undef,l)
    for i in 1:l
        spin = state[i]
        matchain[i] = mps.site[i].data[:,:,spin]
    end
    inner = matchain[1]
    for i in 2:l
        inner = inner*matchain[i]
    end
    return inner[1,1]
end

function calc_Es(sample::Sample, mps::MPS, spH::Dict)::Float64
    state = sample.state
    S = mps.site
    l = length(S)
    Es_diag = 0.0
    Es_offdiag = 0.0
    W = sample.inner
    for i in 1:l-1
        spin = (state[i], state[i+1])
        for (s1,s2,e) in spH[spin]
            if (s1, s2) == spin
                Es_diag += e
            else
                dual_state = deepcopy(state)
                dual_state[i:i+1] = [s1,s2]
                W′ = state_inner(dual_state,mps)
                Es_offdiag += W′*e
            end
        end
    end
    # add the following code to change OBC to PBC
#     if true
#         spin = (state[l], state[1])
#         for (s1,s2,e) in spH[spin]
#             if (s1, s2) == spin
#                 Es_diag += e
#             else
#                 dual_state = deepcopy(state)
#                 dual_state[l] = s1
#                 dual_state[1] = s2
#                 W′ = state_inner(dual_state,mps)
#                 Es_offdiag += W′*e
#             end
#         end
#     end
    # end code

    return Es_diag + Es_offdiag/W
end

# not a gradient actually, it just contracts MPS's matrix except the given one(with index n).
function single_site_grad(state, mps::MPS, H, n)
    S = mps.site
    l = length(state)
    if length(S) != l
        throw(DomainError( (l,length(mps.site)), "length dismatch" ) )
    end
    matchain = Array{Matrix{Float64},1}(undef,l)
    for i in 1:l
        spin = state[i]
        matchain[i] = S[i].data[:,:,spin]
    end
    lmat = 1.0
    rmat = 1.0
    for i in 1:n-1
        lmat = lmat*matchain[i]
    end
    for i in n+1:l
        rmat = rmat*matchain[i]
    end
    G = zeros(Float64,size(S[n].data))
    G[:,:,state[n]] = lmat' * rmat'
    return G
end

function tr2vec(Tree::Vector{Node},node_number)
    i = node_number
    j = mps_size
    chain = Vector{Int8}(undef, mps_size)
    while j != 0
        chain[end+1-j] = Tree[i].spin
        i = Tree[i].parent
        j -= 1
    end
    return chain
end

# from 2^k-1 to 2(2^k-1) are samples' <S'|H|S >. S: temporary(tmp) state.
function generate_samples!(Tree,tmp,H,mps)
    K = mps_size
    for k in 2:K
        for i in 2^k-1:2*(2^k-1)
            j = Tree[i].parent
            s11 = Tree[j].spin
            s12 = Tree[i].spin
            s21 = tmp.state[k-1]
            s22 = tmp.state[k]
            Tree[i].value = Tree[j].value +   H[s11,s12,s21,s22]
        end
    end
    samples = Vector{Sample}(undef,2^K)
    D = view(Tree,2^K-1:2*(2^K-1))
    for k in 1:2^K
        state = tr2vec(Tree,k+2^K-2)
        samples[k] = Sample(state,state_inner(state,mps),D[k].value)

    end
    return samples
end

function gradSampleGenerate(mps::MPS, H, chainlen)
    S = mps.site
    Ns = length(S)
    chain = rand(Int8.(1:spin_deg),Ns)
    tmp = Sample(chain,1e-5,0.0)
    new = tmp
    samples = Vector{Sample}(undef,chainlen)

    for i in 1:Int(ceil(chainlen/2.5))
        new = Sample(next(tmp.state), 0.0, 0.0)
        new.inner = state_inner(new, mps)
        if rand() < (new.inner/tmp.inner)^2
            tmp = new
        end
        # heating, do nothing here
    end

    samples = Vector{Sample}(undef,chainlen)
    Es_sum = 0.0   # ∑ E(S)
    grad = Vector{Array{Float64,3}}(undef,mps_size)
    for i in 1:mps_size
        # will be used to store ∑2/W(S)|S>
        grad[i] = zero(S[i].data)
    end
    delta = zero.(grad)

    tmp.energy = calc_Es(tmp, mps, spH)
    for i in 1:chainlen
        new = Sample(next(tmp.state), 0.0, 0.0)
        new.inner = state_inner(new, mps)
        if rand() < (new.inner/tmp.inner)^2
            tmp = new
            tmp.energy = calc_Es(tmp, mps, spH)
        end

        # tmp.energy: E(S)
        Es_sum += tmp.energy
        samples[i] = tmp
        for i in 1:mps_size
            B = 2/tmp.inner*single_site_grad(tmp.state, mps, H, i)
            delta[i] += B
            grad[i] += B*tmp.energy
        end
    end
    grad ./= chainlen
    delta ./= chainlen
    Es_sum /= chainlen
    grad = grad .- Es_sum.*delta
    return grad, Es_sum, samples
end

function exact_grad(mps::MPS, H)
    S = mps.site
    K = length(S)
    D = view(Tree,2^K-1:2*(2^K-1))
    # list all possible states in samples[], we only use samples[i].state
    samples = Vector{Sample}(undef,2^K)
    for k in 1:2^K
        state = tr2vec(Tree,k+2^K-2)
        samples[k] = Sample(state,0.0,0.0)
    end
    for tmp in samples
        # W(S)
        tmp.inner = state_inner(tmp, mps)
        # E(S)
        tmp.energy = calc_Es(tmp, mps, spH)
    end
    grad = gradFromSet(mps, H, samples)
    return grad
end

function gradFromSet(mps::MPS, mps_old::MPS, H, set)
    S_new = mps.site
    S_old = mps_old.site
    Es_sum = 0.0   # ∑ E(S)
    psum = 0.0 # Σ ρ'/ρ
    grad = Vector{Array{Float64,3}}(undef,mps_size)
    for i in 1:mps_size
        # ∑ 2/W(S)|S>
        grad[i] = zero(S_new[i].data)
    end
    delta = zero.(grad)
    for sp in set
        rate = (state_inner(sp,mps)/sp.inner)^2
        psum += rate

        energy1 = calc_Es(sp, mps, spH)
        Es_sum += energy1*rate

        for i in 1:mps_size
            B = 2/sp.inner*single_site_grad(sp.state, mps, H, i)
            delta[i] += B*rate
            grad[i] += B*sp.energy*rate
        end
    end
    grad ./= psum
    delta ./= psum
    Es_sum /= psum
    grad = grad .- Es_sum.*delta
    return grad
end

function gradFromSet(mps::MPS, H, set)
    S = mps.site
    Es_sum = 0.0   # ∑ E(S)
    grad = Vector{Array{Float64,3}}(undef,mps_size)
    for i in 1:mps_size
        # ∑ 2/W(S)|S>
        grad[i] = zero(S[i].data)
    end
    delta = zero.(grad)
    for sp in set
        Es_sum += sp.energy

        for i in 1:mps_size
            B = 2/sp.inner*single_site_grad(sp.state, mps, H, i)
            delta[i] += B
            grad[i] += B*sp.energy
        end
    end
    chainlen = length(set)
    grad ./= chainlen
    delta ./= chainlen
    Es_sum /= chainlen
    grad = grad .- Es_sum.*delta
    return grad
end

# OBC
function mkSGD(mps1::MPS, H; step=20,chainlen=200,η=0.005, progress=false)
    mps = deepcopy(mps1)
    val = Vector{Float64}(undef,step+1)
    if progress
        progress1 = Progress(step)
    end
    energy1 = 0.0
    val[1] = exact_energy!(mps)
    for j in 1:step
        grad, energy1, _ = gradSampleGenerate(mps, H, chainlen)
        g = grad
        for i in 1:mps_size
#             mps.site[i].data[:] = mps.site[i].data - η*g[i]
            mps.site[i].data[:] = mps.site[i].data - η*(1/(i/500+1))*g[i]
        end
        val[j+1] = exact_energy(mps)
#         val[j] = energy1
        progress && next!(progress1)
    end
    return val, mps
end

# OBC
function mkAdam(mps1::MPS, H; step=20,chainlen=200,η=0.005, β1=0.9, β2=0.999, ϵ=1e-8, progress=false)
    mps = deepcopy(mps1)
    val = Vector{Float64}(undef,step+1)
#     val[1] = exact_energy!(mps)
    if progress
        progress1 = Progress(step)
    end
    energy1 = 0.0
    grad, _, _ = gradSampleGenerate(mps, H, 2)
    g = zero.(grad)
    v = zero.(g)
    for j in 1:step
        grad, energy1, _ = gradSampleGenerate(mps, H, chainlen)
        for i in 1:mps_size
            g[i] = β1*g[i]+(1-β1)*grad[i]
            v[i] = β2*v[i]+(1-β2)*grad[i].^2
            ηt = η*sqrt(1-β2^j)/(1-β1^j)
#             mps.site[i].data[:] = mps.site[i].data - η*g[i]
            mps.site[i].data[:] = mps.site[i].data - ηt*(g[i]./(sqrt.(v[i]).+ϵ*sqrt(1-β2^j)))
        end
        val[j+1] = exact_energy(mps)
#         val[j] = energy1
        progress && next!(progress1)
    end
    return val, mps
end

function mkReweightSVRG(mps1::MPS, H; outer=20,inner=10,long=200,short=40,η=0.005, progress=false)
    mps = deepcopy(mps1)
    val = Vector{Float64}(undef,inner*outer+1)
    val[1] = exact_energy!(mps)
    if progress
        progress1 = Progress(step)
    end
    energy1 = 0.0
    for i in 1:outer
        mps_old = deepcopy(mps)
        grad, energy1, set = gradSampleGenerate(mps, H, long)
        for j in 1:inner

#             _, subset = gradSampleGenerate(mps,H,short)
            partlink = StatsBase.sample(1:long,short,replace=false,ordered=true)
            subset = set[partlink]

            g = grad-gradFromSet(mps_old,H,subset)+gradFromSet(mps_old,mps,H,subset)
#             g = g - ψ'*g*ψ
#             mps.site = mps.site .- η.*g
            for i in 1:mps_size
                mps.site[i].data[:] = mps.site[i].data - η*g[i]
            end
            val[(i-1)*inner+j+1] = exact_energy(mps)
#             val[(i-1)*inner+j+1] = energy1
        end
        progress && next!(progress1)
    end
    return val, mps
end

function tgReweightSVRG(mps1::MPS, H; outer=2,inner=10,long=200,short=40,η=0.005)
    mps = deepcopy(mps1)
    val = Vector{Float64}(undef,inner*outer+1)
    val[1] = exact_energy!(mps)
    progress1 = Progress(outer)
    energy1 = 0.0
    for i in 1:outer
        mps_old = deepcopy(mps)
        _, energy1, set = gradSampleGenerate(mps, H, long)
        grad = exact_grad(mps, H)
        for j in 1:inner
#             _, subset = gradSampleGenerate(mps,H,short)
            partlink = StatsBase.sample(1:long,short,replace=false,ordered=true)
            subset = set[partlink]
            g = grad-gradFromSet(mps_old,H,subset)+gradFromSet(mps_old,mps,H,subset)
#             g = g - ψ'*g*ψ
#             mps.site = mps.site .- η.*g
            for i in 1:mps_size
                mps.site[i].data[:] = mps.site[i].data - η*g[i]
            end
            val[(i-1)*inner+j+1] = exact_energy(mps)
#             val[(i-1)*inner+j+1] = energy1
        end
        next!(progress1)
    end
    return val, mps
end


text:
########

function get_rmse_in_time(operator_data_vec1::OperatorVectorArray, operator_data_vec2::OperatorVectorArray)
    t_count, M, K = size(operator_data_vec1)
    c = 0
    for t_i in 1:t_count
        c += norm(operator_data_vec1[t_i, :, :] - operator_data_vec2[t_i, :, :])^2
    end
    return sqrt(c / t_count)
end

function get_rmse_in_time(operator_data_vec1::OperatorVectorArray, operator_vec2::OperatorVector)
    operator_data_vec2 = operator_recast(operator_vec2)
    return get_rmse_in_time(operator_data_vec1, operator_data_vec2)
end

function get_rmse_in_time(operator_vec1::OperatorVector, operator_data_vec2::OperatorVectorArray)
    operator_data_vec1 = operator_recast(operator_vec1)
    return get_rmse_in_time(operator_data_vec1, operator_data_vec2)
end

function get_rmse_in_time(operator_vec1::OperatorVector, operator_vec2::OperatorVector)
    operator_data_vec1 = operator_recast(operator_vec1)
    operator_data_vec2 = operator_recast(operator_vec2)
    return get_rmse_in_time(operator_data_vec1, operator_data_vec2)
end

#######

function compare_rho(rho::Array, rho_ref::Array; relative=false, smooth_const=1e-9)
    N, M, K = size(rho)
    rho_sum = zeros(Float64, M, K)

    for t_i in 1:N
        if relative
            rho_abs = abs.(rho_ref[t_i, :, :])
            rho_d = rho_abs + smooth_const*ones(size(rho_abs))
            rho_sum[:, :] += abs.(rho[t_i, :, :] - rho_ref[t_i, :, :]) ./ rho_d
        else
            rho_sum[:, :] += abs.(rho[t_i, :, :] - rho_ref[t_i, :, :])
        end
    end
    rho_sum /= N
    return rho_sum
end

function compare_rho(rho::OperatorVectorArray, rho_ref::OperatorVector; relative=false, smooth_const=1e-9)
    rho_ref_array = operator_recast(rho_ref)
    return compare_rho(rho, rho_ref_array; relative=relative, smooth_const=smooth_const)
end

function compare_rho(rho::OperatorVector, rho_ref::OperatorVectorArray; relative=false, smooth_const=1e-9)
    rho_array = operator_recast(rho)
    return compare_rho(rho_array, rho_ref; relative=relative, smooth_const=smooth_const)
end

function compare_rho(rho::OperatorVector, rho_ref::OperatorVector; relative=false, smooth_const=1e-9)
    rho_array = operator_recast(rho)
    rho_ref_array = operator_recast(rho_ref)
    return compare_rho(rho_array, rho_ref_array; relative=relative, smooth_const=smooth_const)
end

######

function compare_rho_in_time(rho::OperatorVectorArray, rho_ref::OperatorVectorArray; smooth_const=1e-9)
    N, M, K = size(rho)
    rho_rel = zeros(Float64, N, M, K)

    for t_i in 1:N
        rho_abs = abs.(rho_ref[t_i, :, :])
        rho_d = rho_abs + smooth_const*ones(size(rho_abs))
        rho_rel[t_i, :, :] = abs.(rho[t_i, :, :] - rho_ref[t_i, :, :]) ./ rho_d
    end
    return rho_rel
end

function compare_rho_in_time(rho::OperatorVectorArray, rho_ref::OperatorVector; smooth_const=1e-9)
    rho_ref_array = operator_recast(rho_ref)
    return compare_rho_in_time(rho, rho_ref_array; smooth_const=smooth_const)
end

function compare_rho_in_time(rho::OperatorVector, rho_ref::OperatorVectorArray; smooth_const=1e-9)
    rho_array = operator_recast(rho)
    return compare_rho_in_time(rho_array, rho_ref; smooth_const=smooth_const)
end

function compare_rho_in_time(rho::OperatorVector, rho_ref::OperatorVector; smooth_const=1e-9)
    rho_ref_array = operator_recast(rho_ref)
    rho_array = operator_recast(rho)
    return compare_rho_in_time(rho_array, rho_ref_array; smooth_const=smooth_const)
end

max_stars_repo_name: detrin/OpenQuantumSystems

max_issues_repo_name: detrin/OpenQuantumSystems

max_forks_repo_name: detrin/OpenQuantumSystems


text:
using Test

using EllipticCurveCrypto

@testset "Field" begin
    f_prime = PrimeField{Int}(2 ^ 13 - 1)
    f_real = RealField{Float64}()
    f_rational = RationalField{Rational{BigInt}}()

    rand_field_elem(f::PrimeField{Int}) = f_from(f, rand(Int))
    rand_field_elem(_f::RealField{Float64}) = rand()
    rand_field_elem(_f::RationalField{Rational{BigInt}}) = rationalize(BigInt, rand())

    for f in [f_prime, f_real, f_rational]
        for _ in 1:10
            a = rand_field_elem(f)
            b = rand_field_elem(f)
            c = rand_field_elem(f)
            n = rand(0:10)

            @test begin
                f_add(f, a, b) == f_add(f, b, a) &&
                f_sub(f, a, b) == f_add(f, a, f_neg(f, b)) &&
                f_mul(f, a, b) == f_mul(f, b, a) &&
                f_mul(f, f_add(f, a, b), c) ≈ f_add(f, f_mul(f, a, c), f_mul(f, b, c)) &&
                f_add(f, a, f_neg(f, a)) == f_zero(f) &&
                f_add(f, f_add(f, a, b), c) == f_add(f, a, f_add(f, b, c))
            end

            if b != f_zero(f)
                @test f_mul(f, a, b) ≈ f_div(f, a, f_inv(f, b))
            end

            # Power
            acc = f_one(f)
            for i in 1:n
                acc = f_mul(f, acc, a)
            end
            @test acc ≈ f_pow(f, a, n)
        end
    end
end

@testset "Elliptic Curve - Smoke" begin
    # From <https://github.com/ashutosh1206/Crypton/blob/master/Elliptic-Curves>
    p = 89953523493328636138979614835438769105803101293517644103178299545319142490503
    a = 89953523493328636138979614835438769105803101293517644103178299545319142490500
    b = 28285296545714903834902884467158189217354728250629470479032309603102942404639
    ec = EllipticCurve(PrimeField(p), a, b)
    @test is_on_curve(ec, ec_zero(ec))

    function test_on_curve(x, y)
        @test is_on_curve(ec, Point(x, y))
        @test !is_on_curve(ec, Point(x + 1, y))
        @test !is_on_curve(ec, Point(x, y + 1))
    end

    x = 23292248698455723400586711958086172976663011248265958495236470661012001058264
    y = 11070572596091667736222067788474181728427994902322868933693472122188841639450
    test_on_curve(x, y)

    x = 19111966210181710119684383106309085309943311037354698052067443650291872435169
    y = 21947619442035360436009014982909893928179456172279297680575771006406194394019
    test_on_curve(x, y)

    x = 85505267126160074714029649657702688395481735744072668398817875879304329931353
    y = 15732931607183825022714544854432024236050897772129106298237207697160928704178
    test_on_curve(x, y)
end

@testset "Elliptic Curve - Randomized" begin
    f = PrimeField{Int}(37)
    a = -5
    b = 8
    ec = EllipticCurve(f, a, b)

    function rand_point(ec::EllipticCurve{T}, order::Int, g::Point{T})::Point{T} where T
        k = rand(0:(order - 1))
        ec_scalar_mul(ec, g, k)
    end

    z = ec_zero(ec)
    g = Point{Int}(1, 2)
    for _ in 1:20
        p1 = rand_point(ec, 37, g)
        p2 = rand_point(ec, 37, g)
        p3 = rand_point(ec, 37, g)

        @test begin
            ec_add(ec, p1, z) == p1 &&
            ec_add(ec, p1, ec_neg(ec, p1)) == z &&
            ec_add(ec, ec_add(ec, p1, p2), p3) == ec_add(ec, p1, ec_add(ec, p2, p3)) &&
            ec_add(ec, p1, p2) == ec_add(ec, p2, p1)
        end

        # Power
        n = rand(0:10)
        p = p1
        acc = z
        for i in 1:n
            acc = ec_add(ec, acc, p)
        end
        @test acc == ec_scalar_mul(ec, p, n)
        @test acc == ec_scalar_mul_ml(ec, p, n)
    end
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
#ДАНО: Робот находится в произвольной клетке ограниченного прямоугольного поля с внутренними перегородками.
#РЕЗУЛЬТАТ: Робот — в исходном положении в центре прямого креста из маркеров, расставленных вплоть до внешней рамки.
function mark_krest(r)
    for i = 0:3
        mark_line(r,i)
    end
    putmarker!(r)
end
function mark_line(r,HS)
    y = 0
    flag = true
    while !isborder(r,HorizonSide(HS % 4)) || flag
        if isborder(r,HorizonSide(HS % 4))==0
            move!(r,HorizonSide(HS % 4))
            putmarker!(r)
            y += 1
        end
        if isborder(r,HorizonSide(HS))
            flag, y = round(r,HS,y)
            putmarker!(r)
        end
        if flag == false
            break
        end
    end
    return(r,HS,y)
end
function return(r,HS,y)
    x = 0
    flag = true
    while isborder(r,HorizonSide(HS)) && !isborder(r,HorizonSide((HS + 1) % 4))
        move!(r,HorizonSide((HS + 1) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(HS)) && isborder(r,HorizonSide((HS + 1) % 4))==1
        flag = false
    end
    if  isborder(r,HorizonSide(HS))==0
        move!(r,HorizonSide(HS))
        y += 1
    end
    while !isborder(r,HorizonSide(HS)) && isborder(r,HorizonSide((HS + 3) % 4))
        move!(r,HorizonSide(HS))
        y += 1
    end
    for i = 1:x
        move!(r,HorizonSide((HS+3) % 4))
    end
    return flag,y
end
mark_krest(r)


text:
using GAIO
using StaticArrays
using Test

@testset "exported functionality" begin
    # a toy example with an easy attractor to test how well all the components
    # fit together
    f(x) = (x[1], x[2] * 0.5)
    test_points = [(-1.0, -1.0), (-1.0, 1.0), (1.0, -1.0), (1.0, 1.0)]
    center = (0.0, 0.0)
    radius = (1.0, 1.0)
    domain = Box(center, radius)
    g = PointDiscretizedMap(f, domain, test_points)
    partition = BoxPartition(domain)
    n = 10
    dims = (32, 32)
    partition_at_depth_n = BoxPartition(domain, dims)
    rga = relative_attractor(g, partition[:], steps = n)
    unstable = unstable_set!(g, partition_at_depth_n[:])
    # ground truths attractor and unstable set
    x_axis = [SVector(0, x) for x in range(-1, 1, length=100)]
    y_axis = [SVector(x, 0) for x in range(-1, 1, length=100)]
    gt_rga = partition_at_depth_n[y_axis]
    gt_unstable = partition_at_depth_n[x_axis]
    # make sure that the algorithms cover the ground truth (we won't have equality)
    @testset "relative attractor" begin
        @test length(gt_rga) > 0
        @test length(intersect(rga, gt_rga)) == length(gt_rga)
    end
    @testset "unstable set" begin
        @test length(gt_unstable) > 0
        @test length(intersect(unstable, gt_unstable)) == length(gt_unstable)
    end
end


text:
using SpikingNetworks

parseParameters"""
τ = 10
σ=80
Vₜ=0.8
Vᵣ=0
N=100
W=0.1*(ones(N,N)-diagm(ones(N)))"""

parseEquations"""
dv/dt = (Iₑ(t)-v+σ*ξ(t)+W⊙sigmoid(v-0.5))/τ
Iₑ(t)=0.5"""

net=init"""
N=N
v=rand(Vᵣ:0.01:Vₜ,N)"""

spike"v.>Vₜ"
reset"v[i]=Vᵣ"

net.dt=0.1
@time run!(net,100/net.dt)
time,rate=activity(net)
using MatlabPlot
mplot(time,rate)


text:
struct ModifiedEquinoctialElements <: AbstractElements
    p::Float64 # Parameter
    f::Float64 # These
    g::Float64 # Don't
    h::Float64 # Have
    k::Float64 # Names
    l::Float64 # True Longitude
end

# Faster than defaults
parameter(m::ModifiedEquinoctialElements) = m.p
true_longitude(m::ModifiedEquinoctialElements) = m.l

# Conversion from Keplerian Elements
function ModifiedEquinoctialElements(ke::KeplerianElements)
    ecc = eccentricity(ke)
    inc = inclination(ke)
    aop = argument_of_periapsis(ke)
    raan = right_ascension(ke)
    sr, cr = sincos(raan)
    sar, car = sincos(aop + raan)
    ti = tan(inc / 2)
    ModifiedEquinoctialElements(
        parameter(ke),
        ecc * car,
        ecc * sar,
        ti * cr,
        ti * sr,
        true_longitude(ke)
    )
end

# Interface
semi_major_axis(m::ModifiedEquinoctialElements) = parameter(m) / (1 - m.f^2 - m.g^2)
eccentricity(m::ModifiedEquinoctialElements) = sqrt(m.f^2  + m.g^2)
inclination(m::ModifiedEquinoctialElements) = atan(2*sqrt(m.h^2 + m.k^2), 1-m.h^2-m.k^2)
right_ascension(m::ModifiedEquinoctialElements) = atan(m.k, m.h)
true_anomaly(m::ModifiedEquinoctialElements) = m.l - atan(m.g / m.f)

function argument_of_periapsis(m::ModifiedEquinoctialElements)
    num = m.g * m.h - m.f * m.k
    den = m.f * m.h + m.g * m.k
    return atan(num, den)
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, MagnitudeDistributions, Distributions, Random

makedocs(sitename="MagnitudeDistributions.jl",
         pages=[
           "Home"=>"index.md",
           "Magnitude of completeness"=>"completeness.md",
           "Frequency-magnitude distributions"=>"distributions.md",
           "How-to"=>["Magnitude of completeness"=>"howto/mc.md",
                      "b-value"=>"howto/b.md"],
           "Utils"=>"utils.md",
              ])

deploydocs(
    repo = "github.com/riccmin/MagnitudeDistributions.jl.git",
    devbranch = "main"
)


text:
"""
Construct active power DeviceRangeConstraintInputs for specific types.
"""
function make_active_power_constraints_inputs(
    ::Type{T},
    ::Type{U},
    ::Type{V},
    ::Union{Nothing, AbstractAffectFeedForward},
    ::Bool,
    ::Bool,
) where {T <: PSY.Device, U <: AbstractDeviceFormulation, V <: PM.AbstractPowerModel}
    error("make_active_power_constraints_inputs is not implemented for types $T / $U / $V")
end

"""
Default implementation to add active power constraints.

Users of this function must implement a method for
[`make_active_power_constraints_inputs`](@ref) for their specific types.
Users may also implement custom active_power_constraints! methods.
"""
function active_power_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, U},
    ::Type{V},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.Device, U <: AbstractDeviceFormulation, V <: PM.AbstractPowerModel}
    use_parameters = model_has_parameters(psi_container)
    use_forecasts = model_uses_forecasts(psi_container)
    @assert !(use_parameters && !use_forecasts)
    inputs = make_active_power_constraints_inputs(
        T,
        U,
        V,
        feedforward,
        use_parameters,
        use_forecasts,
    )
    device_range_constraints!(psi_container, devices, model, feedforward, inputs)
end


text:
@polly function kernel_durbin(r, y)
    n = size(r,1)
    z = zeros(eltype(r),n)
    y[1] = -r[1]
    beta = one(eltype(y))
    alpha = -r[1]

    for k = 2:n
        beta = (1 - alpha * alpha) * beta
        sum = zero(eltype(y))
        for i = 1:(k-1)
            sum += r[k-i] * y[i]
        end
        alpha = - (r[k] + sum)/beta

        for i = 1:(k-1)
            z[i] = y[i] + alpha * y[k-i]
        end
        for i = 1:(k-1)
            y[i] = z[i]
        end
        y[k] = alpha
    end
end

let
    n = 2000

    r = zeros(Float32, n)
    y = zeros(Float32, n)

    for i = 1:n
        r[i] = (n+1-i)
    end

    SUITE["durbin"] = @benchmarkable kernel_durbin(r, y) setup = (r = copy($r); y = copy($y))
end


text:
using MeasureBase: logdensityof, Likelihood
export RadioLikelihood, logdensityof, MultiRadioLikelihood
using LinearAlgebra

struct RadioLikelihood{T,A} <: MB.AbstractMeasure
    lklhds::T
    ac::A
end

struct MultiRadioLikelihood{L} <: MB.AbstractMeasure
    lklhds::L
end

"""
    `MultiRadioLikelihood(lklhd1, lklhd2, ...)`
Combines multiple likelihoods into one object that is useful for fitting multiple days/frequencies.

```julia
lklhd1 = RadioLikelihood(dcphase1, dlcamp1)
lklhd2 = RadioLikelihood(dcphase2, dlcamp2)

lklhd = MultiRadioLikelihood(lklhd1, lklhd2)
```

"""
MultiRadioLikelihood(lklhds::RadioLikelihood...) = MultiRadioLikelihood(lklhds)

function MB.logdensityof(lklhds::MultiRadioLikelihood, m)
    sum(Base.Fix2(logdensityof, m), lklhds.lklhds)
end

"""
    `RadioLikelihood(data1, data2, ...)`
Forms a radio likelihood from a set of data products. These data products must share
the same array data/configuration. If you want to form a likelihood from multiple arrays
such as when fitting different wavelengths or days, you can combine them using
`MultiRadioLikelihood`

```julia
lklhd1 = RadioLikelihood(dcphase1, dlcamp1)
lklhd2 = RadioLikelihood(dcphase2, dlcamp2)

lklhd = MultiRadioLikelihood(lklhd1, lklhd2)
```
"""
function RadioLikelihood(data::EHTObservation...)
    ls = Tuple(map(makelikelihood, data))
    acs = arrayconfig.(data)
    #@argcheck acs[1] == acs[2]
    RadioLikelihood{typeof(ls), typeof(acs[1])}(ls, acs[1])
end

function RadioLikelihood(data::MT.Likelihood...)
    return RadioLikelihood{typeof(data), Nothing}(data, nothing)
end

function Base.show(io::IO, d::RadioLikelihood{T}) where {T}
    println(io, "RadioLikelihood{$T}")
    println(io, "\tNumber of data products: ", length(d.lklhds))
end

"""
    `logclosure_amplitudes(vis, ac::ArrayConfiguration)`

Compute the log-closure amplitudes for a set of visibilities and an array configuration

# Notes
This uses a closure design matrix for the computation.
"""
function logclosure_amplitudes(vis::AbstractArray{<:Complex}, ac::ArrayConfiguration)
    lva = log.(abs.(vis))
    return ac.designmat*lva
end

"""
    `closure_phases(vis, ac::ArrayConfiguration)`

Compute the closure phases for a set of visibilities and an array configuration

# Notes
This uses a closure design matrix for the computation.
"""
function closure_phases(vis::AbstractArray{<:Complex}, ac::ArrayConfiguration)
    ph = angle.(vis)
    return ac.designmat*ph
end

amplitudes(vis::AbstractArray{<:Complex}, ac::ArrayConfiguration) = abs.(vis)
phase(vis::AbstractArray{<:Complex}) = angle.(vis)

function MB.logdensityof(d::RadioLikelihood, m::ComradeBase.AbstractModel)
    ac = d.ac
    vis = visibilities(m, ac)
    return logdensityof(d, vis)
end

function MB.logdensityof(d::RadioLikelihood, vis::AbstractArray)
    # We use a for loop here since Zygote plays nice with this
    acc = logdensityof(first(d.lklhds), vis)
    @inbounds for l in d.lklhds[begin+1:end]
        acc += logdensityof(l, vis)
    end
    return acc
end

#function MB.logdensityof(d::RadioLikelihood, m::ComradeBase.AbstractModel)
#    acc = logdensityof(first(d.lklhds), m)
#    @inbounds for i in 2:length(d.lklhds)
#        acc += logdensityof(d.lklhds[i], m)
#    end
#    return acc
#end

function makelikelihood(data::Comrade.EHTObservation{<:Real, <:Comrade.EHTVisibilityDatum})
    errinv = inv.(data[:error])
    vis = StructArray{Complex{eltype(data[:visr])}}((data[:visr],data[:visi]))
    ℓ = Likelihood(vis) do (μ, )
        ComplexNormal{(:μ, :τ)}(μ, errinv)
    end
    return ℓ
end

function makelikelihood(data::Comrade.EHTObservation{<:Real, <:Comrade.EHTVisibilityAmplitudeDatum})
    τ = inv.(data[:error])
    amp = getdata(data, :amp)
    ℓ = Likelihood(amp) do (μ, )
        AmpNormal{(:μ, :τ)}(abs.(μ), τ)
    end
    return ℓ
end

function makelikelihood(data::Comrade.EHTObservation{<:Real, <:Comrade.EHTLogClosureAmplitudeDatum})
    dmat = data.config.designmat
    τ  = inv.(data[:error])
    #covm = σvis*dmat*transpose(σvis)

    #C = Cholesky(Symmetric(covm))
    #σ = C.L

    f = Base.Fix2(logclosure_amplitudes, data.config)
    amp = data[:amp]
    ℓ = Likelihood(amp) do (μ,)
        AmpNormal{(:μ, :τ)}(f(μ), τ)
    end
    #return AmpNormal{(:μ, :τ)}(amp, τ)
    return ℓ
end

function makelikelihood(data::Comrade.EHTObservation{<:Real, <:Comrade.EHTClosurePhaseDatum})
    τ = inv.(getdata(data, :error)).^2
    f = Base.Fix2(closure_phases, data.config)
    phase = data[:phase]
    ℓ = Likelihood(phase) do (μ,)
        CPVonMises{(:μ, :κ)}(f(μ), τ)
    end

    return ℓ
end


text:
push!(LOAD_PATH, "/Users/megorov/Desktop/projects/stanford/sniper/egorov/sniperRepo/new-version")

using Collada_
using UrbanMap_
using SniperPOMDP_
using SniperMDP_
using ValueIterations
using Helpers
using RayCaster_
using Simulator
using SimServers
using PolicyAssignment

map = UrbanMap(1, [0.5], [0.5], [0.42], [0.22], 10, 10)
mapPath = "/Users/megorov/Desktop/projects/stanford/sniper/egorov/sniperRepo/new-version/maps/demo_map_border2D.dae"
coll = ColladaObjects(mapPath)
#map = UrbanMap(coll, 15, 15) # 15x15 map

#t   = [(2.0,8.0), (2.0,1.0), (5.0,2.0), (7.0,3.0)]
t   = [(2.0,8.0), (2.0,1.0)]
#t   = [(7.0,7.0), (2.0,14.0), (8.0,4.0), (11.0,13.0)]
mdp = SniperMDP(map, nSnipers = 1, nMonitors = 2, targets=t)

p1 = load_policy("policies/mdp-1.policy")
p2 = load_policy("policies/mdp-2.policy")
p3 = load_policy("policies/mdp-3.policy")
p4 = load_policy("policies/mdp-4.policy")

#policies = merge(p1.policy,p2.policy,p3.policy,p4.policy)
policies = merge(p1.policy,p2.policy)
#values   = merge(p1.V,p2.V,p3.V,p4.V)
values   = merge(p1.V,p2.V)

policy = AssignmentPolicy(2, policies, values, mdp)

# create the simulation type
#init  = 8_888_888_803 # sniper at (3,1), all monitors at (9,9) for 10x10 4 agent world
#init = 97 # sniper at (7,7), all monitors at (1,1) for 15x15 4 agent world
#init = 20_583_164_347 # sniper at (7,7), all monitors at (1,i) for 15x15 4 agent world
init = 480_872
sim  = MDPSimulation(mdp, init, 100)

# set up the server
socket = 4444
server = SniperServer(socket)

# begin simulation and send results to client
#simAndSend(server, sim, mdp, policy)

# begin batch sim
n = 10
rewards = zeros(n)
shot    = zeros(Int64, n)
obs     = zeros(Int64, n)
for i = 1:n
    res = simulate(policy, mdp, 100, initState=init)
    rewards[i] = sum(res.rewards)
    shot[i]    = res.nKills
    obs[i]     = res.nObs
end


text:
NaiveNASlib.nin(t::FluxLayer, l) = throw(ArgumentError("Not implemented for $t"))
NaiveNASlib.nout(t::FluxLayer, l) = throw(ArgumentError("Not implemented for $t"))

NaiveNASlib.nin(::FluxParLayer, l) = [size(weights(l), indim(l))]
NaiveNASlib.nout(::FluxParLayer, l) = size(weights(l), outdim(l))
NaiveNASlib.nout(::FluxDepthwiseConv, l) = size(weights(l), outdim(l)) * nin(l)[]

NaiveNASlib.nin(::FluxParInvLayer, l) = [nout(l)]

NaiveNASlib.nout(::FluxDiagonal, l) = length(weights(l))
NaiveNASlib.nout(::FluxParInvLayer, l::LayerNorm) = nout(l.diag)
NaiveNASlib.nout(::FluxParNorm, l) = l.chs

NaiveNASlib.nout(::FluxRecurrent, l) = div(size(weights(l), outdim(l)), outscale(l))

outscale(l) = outscale(layertype(l))
outscale(::FluxRnn) = 1
outscale(::FluxLstm) = 4
outscale(::FluxGru) = 3

indim(l) = indim(layertype(l))
outdim(l) = outdim(layertype(l))
actdim(l) = actdim(layertype(l))
actrank(l) = actrank(layertype(l))

indim(t::FluxLayer) = throw(ArgumentError("Not implemented for $t"))
outdim(t::FluxLayer) = throw(ArgumentError("Not implemented for $t"))
actdim(t::FluxLayer) = throw(ArgumentError("Not implemented for $t"))
actrank(t::FluxLayer) = throw(ArgumentError("Not implemented for $t"))

indim(::Flux2D) = 2
outdim(::Flux2D) = 1
actdim(::Flux2D) = 1
actrank(::Flux2D) = 1

indim(::FluxDiagonal) = 1
outdim(::FluxDiagonal) = 1
actdim(::FluxDiagonal) = 1
actrank(::FluxDiagonal) = 1

indim(::FluxRecurrent) = 2
outdim(::FluxRecurrent) = 1
actdim(::FluxRecurrent) = 1
actrank(::FluxRecurrent) = 2

indim(::FluxConvolutional{N}) where N = 2+N
outdim(::FluxConvolutional{N}) where N = 1+N
actdim(::FluxConvolutional{N}) where N = 1+N
actrank(::FluxConvolutional{N}) where N = 1+N
indim(::Union{FluxConv{N}, FluxCrossCor{N}}) where N = 1+N
outdim(::Union{FluxConv{N}, FluxCrossCor{N}}) where N = 2+N
# Note: Absence of bias mean that bias is of type Flux.Zeros which mostly behaves like a normal array, mostly...
weights(l) = weights(layertype(l), l)
bias(l) = bias(layertype(l), l)

weights(::FluxDense, l) = l.weight
bias(::FluxDense, l) = l.bias

weights(::FluxConvolutional, l) = l.weight
bias(::FluxConvolutional, l) = l.bias

weights(::FluxDiagonal, l) = l.α
bias(::FluxDiagonal, l) = l.β

weights(lt::FluxRecurrent, l::Flux.Recur) = weights(lt, l.cell)
bias(lt::FluxRecurrent, l::Flux.Recur) = bias(lt, l.cell)
weights(::FluxRecurrent, cell) = cell.Wi
bias(::FluxRecurrent, cell) = cell.b

hiddenweights(l) = hiddenweights(layertype(l), l)
hiddenweights(lt::FluxRecurrent, l::Flux.Recur) = hiddenweights(lt, l.cell)
hiddenweights(::FluxRecurrent, cell) = cell.Wh

hiddenstate(l) = hiddenstate(layertype(l), l)
hiddenstate(lt::FluxRecurrent, l::Flux.Recur) = hiddenstate(lt, l.cell)
hiddenstate(::FluxRecurrent, cell) = cell.state0
hiddenstate(::FluxLstm, cell::Flux.LSTMCell) = [h for h in cell.state0]

state(l) = state(layertype(l), l)
state(::FluxRecurrent, l) = l.state
state(::FluxLstm, l) = [h for h in l.state]


text:
####################################################
function makecptfromgrd(G::GMT.GMTgrid; factor_lims=0.8, sigdigits_lims=2,
	                                    T=[], kwargs...)
	isempty(T) && ( T=round.(factor_lims.*extrema(G.z), sigdigits=sigdigits_lims) )
	return GMT.makecpt(; T=T, kwargs...)
end
####################################################

####################################################
function gmttopo(G::GMT.GMTgrid; factor_lims=0.8, sigdigits_lims=2,
	                             C=:geo, T=[], D::Bool=true, J="", R="", kwargs...)
	## cpt
	cptout = false
	if !isa(C, GMT.GMTcpt)
		cptout = true
		C = makecptfromgrd(G; factor_lims=factor_lims, sigdigits_lims=sigdigits_lims, C=C, T=T, D=D)
	end
    # options
	isempty(J) && ( J=getJ("X10", axesratio(G)) )
	isempty(R) && ( R=getR(G) )

    ## plot
    GMT.grdimage(G; C=C, J=J, R=R, Q=true, kwargs...)

    ## return
	if cptout; return C; else return nothing; end
end
####################################################
gmttopo(topo::VisClaw.Topo; kwargs...) = gmttopo(geogrd(topo); kwargs...)
####################################################

####################################################
"""
    gmttoporange!(geo::VisClaw.AbstractTopo; kwargs...)

plot a range of topo/bath using GMT
"""
function gmttoporange!(geo::VisClaw.AbstractTopo; kwargs...)
	# set square
	xp = [geo.x[1],  geo.x[1]  , geo.x[end], geo.x[end], geo.x[1]]
	yp = [geo.y[1],  geo.y[end], geo.y[end], geo.y[1]  , geo.y[1]]
	# plot
	GMT.plot!(xp, yp; marker=:none, kwargs...)
end
####################################################

####################################################
"""
    gmtcoastline(topo::VisClaw.Topo; kwargs...)
	gmtcoastline(G::GMT.GMTgrid; kwargs...)
    gmtcoastline!(topo::VisClaw.Topo; kwargs...)
	gmtcoastline!(G::GMT.GMTgrid; kwargs...)

plot coastlines from topography and bathymetry data using GMT
"""
gmtcoastline!(topo::VisClaw.Topo; kwargs...) = GMT.grdcontour!(geogrd(topo); C="-1e10,0,1e10", kwargs...)
####################################################
gmtcoastline!(G::GMT.GMTgrid; kwargs...) = GMT.grdcontour!(G; C="-1e10,0,1e10", kwargs...)
####################################################

####################################################
"""
$(@doc gmtcoastline!)
"""
gmtcoastline(topo::VisClaw.Topo; kwargs...) = GMT.grdcontour(geogrd(topo); C="-1e10,0,1e10", kwargs...)
####################################################
gmtcoastline(G::GMT.GMTgrid; kwargs...) = GMT.grdcontour(G; C="-1e10,0,1e10", kwargs...)
####################################################


text:
function VegaStream()
    VegaStream(nothing,nothing,nothing)
end

function VegaDataRef()
    VegaDataRef(nothing,nothing,nothing,nothing)
end

function VegaValueRef()
    VegaValueRef(nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing)
end

function VegaPadding()
    VegaPadding(80,80,80,80)
end

function VegaMarkPropertySet()
    VegaMarkPropertySet(nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing)
end

function VegaAxis()
    VegaAxis("x","x",nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,"front",false,nothing)
end

function VegaFormat()
    VegaFormat(nothing,nothing,nothing,nothing,nothing,nothing)
end

function VegaData()
    VegaData("table",nothing,nothing,nothing,nothing,nothing,nothing)
end

function VegaMarkProperties()
    VegaMarkProperties(nothing,nothing,nothing,nothing)
end

function VegaMarkFrom()
    VegaMarkFrom(nothing,nothing,nothing,nothing)
end

function VegaSignal()
    VegaSignal(nothing,nothing,false,nothing,nothing,nothing)
end

function VegaScale()
    VegaScale("x","linear",nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing)
end

function VegaMark()
    VegaMark(nothing,nothing,"symbol",nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing)
end

function VegaLegend()
    VegaLegend(nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing)
end

function VegaScene()
    VegaScene(nothing,nothing,nothing,nothing,nothing,nothing,nothing)
end

function VegaContainer()
    VegaContainer(nothing,nothing,nothing,nothing,nothing)
end

function VegaPredicate()
    VegaPredicate(nothing,nothing,nothing,nothing,nothing,nothing)
end

function VegaVisualization()
    VegaVisualization("Vega Visualization",450,450,"auto",nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing,nothing)
end

max_stars_repo_name: jpthor/vega


text:
# Printing -----------------------------

function printpath{T, S}(tn::TreeNode{T, S}, path::Vector{S} = Vector{S}(); sep = "/", io::IO = STDOUT)

    push!(path, tn.id)

    println(io, join(path, sep))

    if !isleaf(tn)
        for (key, chnode) in tn.children
            printpath(chnode, path; sep = sep, io = io)
        end
    end

    pop!(path)

    nothing
end

function printpath{T, S}(tn::TreeNode{T, S}; sep = "/", io::IO = STDOUT)

    path = Vector{S}()

    printpath(tn, path; sep = sep, io = io)
end

"""
# Usage
Prints an UTF-8 formatted representation of the `tree`.
# Examples

```julia
julia> print(tree)
```
"""
function print{T, S}(tn::TreeNode{T, S}; leading::AbstractString = "", io::IO = STDOUT)

    println(io, tn.id)

    if !isleaf(tn)
        last = length(tn.children)
        i::Int = 0
        for (key, chnode) in tn.children
            i += 1
            if (i != last)
                to_print    = leading * "\u251C\u2500 " # ├─
                new_leading = leading * "\u2502  " # │
            else
                to_print    = leading * "\u2514\u2500 " # └─
                new_leading = leading * "   "
            end

            print(io, to_print)
            print(chnode; leading = new_leading, io = io)
        end
    end
end

function print{T, S}(tn::Nullable{TreeNode{T, S}}; args...)
    isnull(tn) ? nothing : print(get(tn); args...)
end


text:
#=
a ball-packing type problem: choose centers and radii of m 2-norm balls into a box [0, 1]ⁿ in
n-dimensional space, maximizing the sum of the radii

constraints are nonconvex, so use SOS2 or PWL techniques with extended formulations
=#

struct BallPacking <: ExampleInstance
    m::Int # number of balls
    n::Int # dimension of space
    pwl::PWLSOS2 # type of piecewise linear SOS2 formulation to use
    num_pts::Int # number of points per piecewise linearization
end

function build(inst::BallPacking)
    (m, n, num_pts) = (inst.m, inst.n, inst.num_pts)
    @assert m >= 2
    @assert n >= 2
    @assert isodd(num_pts) # ensures 0 is included
    pts = collect(range(-1.0, 1.0, length = num_pts))
    @assert 0.0 in pts
    f_pts = abs2.(pts)

    # build model
    model = JuMP.Model()
    JuMP.@variable(model, C[1:m, 1:n]) # centers
    JuMP.@variable(model, 0 <= R[1:m] <= 0.5) # radii
    JuMP.@objective(model, Max, sum(R))

    JuMP.@constraints(model, begin
        # balls inside box
        [i in 1:m, k in 1:n], C[i, k] >= R[i]
        [i in 1:m, k in 1:n], C[i, k] <= 1 - R[i]
        # break some symmetry TODO more
        [i in 1:(m - 1)], C[i, 1] <= C[i + 1, 1]
    end)

    # for each pair of balls i < j, need Rᵢ + Rⱼ ≤ ‖Cᵢ - Cⱼ‖
    for i in 1:m, j in (i + 1):m
        # EF variables and linear constraint
        λ = JuMP.@variable(model, [1:n])
        Rij = JuMP.@variable(model)
        JuMP.@constraint(model, Rij == R[i] + R[j])
        JuMP.@constraint(model, Rij <= sum(λ))

        # 3-dim cone constraints; note Cᵢ - Cⱼ ∈ [-1, 1]ⁿ, ‖Cᵢ - Cⱼ‖ ≤ √n, Rij ∈ [0, 1]
        for k in 1:n
            # auxiliary variables and SOS2 formulation
            σ = JuMP.@variable(model, [1:length(pts)], lower_bound = 0)
            PWL_SOS2(inst.pwl, model, σ, 1.0)

            # data constraints
            JuMP.@constraints(model, begin
                dot(σ, f_pts) == λ[k]
                sum(σ) == Rij
                dot(σ, pts) == C[i, k] - C[j, k]
            end)
        end
    end

    # save for use in tests
    model.ext[:R] = R
    model.ext[:C] = C

    return model
end

function test_extra(inst::BallPacking, model::JuMP.Model)
    stat = JuMP.termination_status(model)
    @test stat in OPT_OR_LIMIT
    (stat == MOI.OPTIMAL) || return

    R = JuMP.value.(model.ext[:R])
    C = JuMP.value.(model.ext[:C])
    # check feasibility for linear constraints
    tol = eps()^0.4
    @test all(-tol .<= R .<= 0.5 + tol)
    @test all(-tol .<= C .<= 1 + tol)
    # check near-feasibility for nonconvex constraints
    tol_loose = eps()^0.1
    approxij(i, j) = (R[i] + R[j] <= norm(C[i, :] - C[j, :]) + tol_loose)
    @test all(approxij(i, j) for i in 1:(inst.m) for j in 1:(i - 1))
    return
end


text:
const DUMMY_PRECISION = 1e-12

Base.isfinite(a::AbstractVec) = all(isfinite.(a))
Base.isinf(a::AbstractVec) = any(isinf.(a))

function StaticArrays.similar_type(::Type{V}, ::Type{F}, size::Size{N}) where {V<:AbstractVec, F, N}
    if size == Size(V) && eltype(V) == F
        return V
    else # convert back to SArray
        return SArray{Tuple{N...},F,length(size),prod(size)}
    end
end

function StaticArrays.similar_type(::Type{V}, ::Type{T}, size::Size{N}) where {V<:VecE2, T, N}
    if size == Size(V)
        return VecE2{T}
    else # convert back to SArray
        return SArray{Tuple{N...},T,length(size),prod(size)}
    end
end

function StaticArrays.similar_type(::Type{V}, ::Type{T}, size::Size{N}) where {V<:VecE3, T, N}
    if size == Size(V)
        return VecE3{T}
    else # convert back to SArray
        return SArray{Tuple{N...},T,length(size),prod(size)}
    end
end

function StaticArrays.similar_type(::Type{V}, ::Type{T}, size::Size{N}) where {V<:VecSE2, T, N}
    if size == Size(V)
        return VecSE2{T}
    else # convert back to SArray
        return SArray{Tuple{N...},T,length(size),prod(size)}
    end
end

Base.abs(a::AbstractVec) = error("abs(v::AbstractVec) has been removed. Use norm(v) to get the norm; abs.(v) to get the element-wise absolute value.")


text:
@testset "at the beggining of equiSpacedPosition" begin
    @test equiSpacedPosition(1, [1,2,3,4,5,6,7,8]) == 1
end

@testset "inside of equiSpacedPosition 1" begin
    @test equiSpacedPosition(3.5, [1,2,3,4,5,6,7,8]) == 3
end

@testset "inside of equiSpacedPosition 2" begin
    @test equiSpacedPosition(2, [1,2,3,4,5,6,7,8]) == 2
end

@testset "at the end of equiSpacedPosition" begin
    @test equiSpacedPosition(8, [1,2,3,4,5,6,7,8]) == 7
end

@testset "the beginning of evalBSplineBasis" begin
    @test all(evalBSplineBasis(1.,buildKnots(1.,10., 12)) .≈ vcat([1.0/6.0, 4.0/6.0, 1.0/6.0, 0.], zeros(Float64,8)))
end

@testset "the end of evalBSplineBasis" begin
    @test all((evalBSplineBasis(10.,buildKnots(1.,10., 12))-vcat(zeros(Float64, 8), [0., 1.0/6.0, 4.0/6.0, 1.0/6.0])).<1e-16)
end

@testset "test evalBSpline" begin
    @test abs(evalBSpline(9.,buildKnots(1.,10., 12), ones(12))-1.0)<=1e-12
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule HelFEM_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("HelFEM")
JLLWrappers.@generate_main_file("HelFEM", UUID("584ea77b-75d0-51ad-9ad6-7fb559b95557"))
end  # module HelFEM_jll


text:
#=
Benchmark pure Julia hash function implementations.
=#
using LSMH
using BenchmarkTools, Random, Statistics

function MurmurHash3_x86_32(s::String, seed::UInt32)
    ccall(:memhash32_seed, UInt32, (Ptr{UInt8}, Csize_t, UInt32), pointer(s),
          sizeof(s), seed % UInt32)
end
MurmurHash3_x86_32(s::String) = MurmurHash3_x86_32(s, UInt32(0))

string_lengths = [5, 10, 50, 100, 500, 1000]

A = zeros(length(string_lengths), 4)

for i = 1:length(string_lengths)
    l = string_lengths[i]
    murmur32(randstring(l))
    xxhash32(randstring(l))
    a = @benchmark murmur32($(randstring(l)))
    b = @benchmark MurmurHash3_x86_32($(randstring(l)))
    c = @benchmark xxhash32($(randstring(l)))
    A[i,1] = l
    A[i,2] = median(a.times)
    A[i,3] = median(b.times)
    A[i,4] = median(c.times)
end

A


text:
using SafeTestsets

if Threads.nthreads() > 1
    @show Threads.nthreads()
else
    @warn("Tests using only 1 thread, multithreaded code will not be tested.")
end

@safetestset "test_CleanTable.jl" begin include("test_CleanTable.jl") end
@safetestset "test_polish_names.jl" begin include("test_polish_names.jl") end
@safetestset "test_compact.jl" begin include("test_compact.jl") end
@safetestset "test_delete_const.jl" begin include("test_delete_const.jl") end
@safetestset "test_row_as_names.jl" begin include("test_row_as_names.jl") end
@safetestset "test_reinfer_schema.jl" begin include("test_reinfer_schema.jl") end
@safetestset "test_drop_missing.jl" begin include("test_drop_missing.jl") end
@safetestset "test_add_index.jl" begin include("test_add_index.jl") end
@safetestset "test_returning_original.jl" begin include("test_returning_original.jl") end
@safetestset "test_utils.jl" begin include("test_utils.jl") end

max_stars_repo_name: TheRoniOne/Cleaner

max_issues_repo_name: TheRoniOne/Cleaner

max_forks_repo_name: TheRoniOne/Cleaner


text:
#include("ConeMantle.jl")
#include("EllipsoidMantle.jl")
#include("EllipticalSurface.jl")
#include("Plane.jl")
#include("Polygon.jl")
#include("TorusMantle.jl")

@recipe function f(s::AbstractSurfacePrimitive; n_arc = 40, n_vert_lines = 2)
    seriestype --> :csg
    if haskey(plotattributes, :seriestype) && plotattributes[:seriestype] == :csg
        @series begin
            label := ""
            linewidth := 0
            linecolor --> :white
            mesh(s, n_arc)
        end
        @series begin
            label --> get_label_name(s)
            linecolor --> :black
            fillalpha := 1
            linewidth --> 1.5
            mesh(s, n_arc, n_vert_lines)
        end
    elseif haskey(plotattributes, :seriestype) && plotattributes[:seriestype] == :mesh3d
        label --> get_label_name(s)
        linecolor --> :white
        mesh(s, n_arc)
    elseif haskey(plotattributes, :seriestype) && plotattributes[:seriestype] == :wireframe
        label --> get_label_name(s)
        seriescolor --> :black
        fillalpha := 1
        linewidth --> 2
        mesh(s, n_arc, n_vert_lines)
    end
end

@recipe function f(vp::AbstractVector{<:AbstractSurfacePrimitive})
    linecolor --> :black
    @series begin
        label --> "Faces"
        show_normal --> false
        vp[1]
    end
    if length(vp) > 1
        for p in vp[2:end]
            @series begin
                show_normal --> false
                label := nothing
                p
            end
        end
    end
end


text:
@testset "Lattices" begin
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x^2-2
  K, a = number_field(f,"a")
  D = matrix(K, 3, 3, [1//64, 0, 0, 0, 1//64, 0, 0, 0, 1//64])
  gensL = [[32, 0, 0], [944*a+704, 0, 0], [16, 16, 0], [72*a+96, 72*a+96, 0], [4*a, 4*a+8, 8], [20*a+32, 52*a+72, 32*a+40]]
  L = @inferred quadratic_lattice(K, gensL, gram = D)
  D = matrix(K, 3, 3, [1//64, 0, 0, 0, 1//64, 0, 0, 0, 1//64])
  gensM = [[32, 0, 0], [720*a+448, 0, 0], [16, 16, 0], [152*a+208, 152*a+208, 0], [4*a+24, 4*a, 8], [116*a+152, 20*a+32, 32*a+40]]
  M = @inferred quadratic_lattice(K, gensM, gram = D)
  @test norm(volume(M))*discriminant(K)^rank(L) == abs(det(restrict_scalars(M)))

  p = prime_decomposition(base_ring(L), 2)[1][1]
  @test @inferred is_locally_isometric(L, M, p)
  @test @inferred Hecke.is_locally_isometric_kirschmer(L, M, p)
  @test is_rationally_isometric(L, M)
  @test is_rationally_isometric(L, M, infinite_place(K, 1))

  q = quadratic_space(K, diagonal_matrix(Hecke.diagonal_of_rational_span(L)))
  @test is_isometric(ambient_space(L),q)
  gensL = generators(L, minimal=true)
  L1 = lattice(ambient_space(L), matrix(gensL))
  @test L1 == L

  L = quadratic_lattice(K, gensM, gram = 9*8*identity_matrix(K,3))
  p = prime_decomposition(base_ring(L), 2)[1][1]
  fl = false
  while !fl && is_integral(norm(L)) #for safety
    fl, Lover = Hecke.is_maximal_integral(L)
    @test ambient_space(Lover) === ambient_space(L)
    L = Lover
  end

  @test Hecke.is_maximal_integral(L, p)[1]
  @test !is_modular(L)[1]
  OK = maximal_order(K)
  p3 = prime_ideals_over(OK, 3)[1]
  @test is_modular(L, p3)[1]
  @test norm(volume(L))*discriminant(OK)^rank(L) == abs(det(restrict_scalars(L)))

  @test ambient_space(dual(L)) === ambient_space(L)
  @test ambient_space(Hecke.lattice_in_same_ambient_space(L,pseudo_matrix(L))) === ambient_space(L)

  B = identity_matrix(K, 3)
  Bp = pseudo_matrix(B)
  # test lazy creation of the ambient space.
  L = quadratic_lattice(K,Bp, gram=D)
  @test dim(ambient_space(L)) == 3
  @test degree(L) == ncols(L.pmat.matrix)
  @test sprint(show, L) isa String
  @test gram_matrix(ambient_space(L)) == D

  L = quadratic_lattice(K,B, gram=D)
  @test dim(ambient_space(L)) == 3
  @test sprint(show, L) isa String
  @test gram_matrix(ambient_space(L)) == D

  # printing
  @test sprint(show, L) isa String

  # Smoke test for genus symbol
  Qx, x = PolynomialRing(FlintQQ, "x")
  K, a = NumberField(x^2 - 2, "a")
  L = @inferred quadratic_lattice(K, identity_matrix(K, 10), gram = 35*identity_matrix(K, 10))
  P = prime_decomposition(maximal_order(K), 5)[1][1]
  #GM = @inferred Hecke._genus_symbol_kirschmer(L, P)
  #@test GM._genus_symbol_kirschmer(L, P).data == [(10, 1, 1)]

  #P = prime_decomposition(maximal_order(K), 2)[1][1]
  #GS = @inferred Hecke._genus_symbol_kirschmer(L, P)
  #@test GS.data[1] == [10]
  #@test GS.data[2] == [0]
  #@test GS.data[3] == [2]
  #@test GS.data[4] == [35]
  #@test GS.data[5] == []

  Qx, x = FlintQQ["x"];
  K, a = number_field(x^2 - 2, cached = true);
  Kt, t = K["t"];
  L, b = number_field(t^2 + 11, "b", check = true)
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  P = prime_decomposition(maximal_order(L), p)[1][1]
  H = @inferred Hecke.lattice(quadratic_space(K, 2 * identity_matrix(K, 3)), pseudo_matrix(identity_matrix(K, 3), [p, p, p]))
  #@test Hecke._genus_symbol_kirschmer(H, fmpz(2)) == Any[(3, 4, true, 4, -64)]
  @test @inferred is_locally_isometric(H, H, p)

  H = Hecke.lattice(hermitian_space(L, L(elem_in_nf(uniformizer(p))) * identity_matrix(L, 3)), pseudo_matrix(identity_matrix(L, 3), [P, P, P]))
  #@test Hecke._genus_symbol_kirschmer(H, p) == Any[(3, 3, false)]
  @test @inferred is_locally_isometric(H, H, p)
  GH = @inferred genus(H, p)
  GH2 = @inferred genus(HermLat, L, p, [(3, 3, -1)])
  @test GH == GH2

  V = @inferred hermitian_space(L, L(a) * identity_matrix(L, 3))
  M = @inferred Hecke.maximal_integral_lattice(V)
  GM = @inferred genus(M, p)
  GM2 = @inferred genus(HermLat, L, p, [(0, 2, 1), (1, 1, -1)])
  @test GM == GM2
  q = prime_decomposition(maximal_order(K), 11)[1][1]
  @test Hecke.genus(M, q) == genus(HermLat, L, q, [(-1, 2, 1), (0, 1, 1)])

  L, _ = number_field(t^2 - gen(K) * t + 1)
  V = hermitian_space(L, L(a) * identity_matrix(L, 3))
  M = @inferred Hecke.maximal_integral_lattice(V)
  @test Hecke.genus(M, p) == genus(HermLat, L, p, [(-2, 2, 1, 0), (0, 1, -1, 0)])

  V = hermitian_space(L, L(10) * identity_matrix(L, 3))
  M = @inferred Hecke.maximal_integral_lattice(V)
  @test Hecke.genus(M, p) == genus(HermLat, L, p, [(-2, 2, 1, 0), (0, 1, 1, 0)])

  K, a = CyclotomicRealSubfield(8, "a")
  Kt, t = K["t"]
  E, b = number_field(t^2 - a * t + 1, "b")
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  P = prime_decomposition(maximal_order(E), p)[1][1]
  ME = maximal_order(E)
  pm = pseudo_matrix(matrix(E, 3, 3, [1, 0, 0, b, 1, 0, 0, 0, 1]), [one(E)*ME, inv(P)^2, one(E)*ME])
  V = @inferred hermitian_space(E, identity_matrix(E, 3))
  L = @inferred lattice(V, pm)
  o = @inferred Hecke.automorphism_group_order(L)
  @test o == 1536

  @test ambient_space(P*L) === ambient_space(L)

  @test is_sublattice(L, P*L)
  @test !is_sublattice(P * L, L)
  @test issubset(P*L, L)
  @test !issubset(L, P*L)
  @test is_integral(L) == is_sublattice(dual(L), L)
  VV = hermitian_space(E, identity_matrix(E, 3), cached = false)
  LL = lattice(VV, pm)
  @test L != LL
  @test !issubset(L, LL)

  # intersections for modules of non-full rank not yet implemented
  K, a = rationals_as_number_field()
  OK = maximal_order(K)
  q = quadratic_space(K, K[1 0; 0 1])
  L = fractional_ideal(OK, K(1//2))*lattice(q)
  S = lattice(q, matrix(generators(L)[1:1]))
  @test_broken @inferred intersect(L, S)
  @test_broken issubset(orthogonal_complement(L,S), L)

  E8 = root_lattice(:E,8)
  L = Hecke._to_number_field_lattice(E8)
  @test L == dual(L)

  K, a = CyclotomicRealSubfield(8, "a")
  Kt, t = K["t"]
  E, b = number_field(t^2 - a * t + 1, "b")
  V = hermitian_space(E, gram_matrix(root_lattice(:E, 8)))
  L = lattice(V)
  @test L == dual(L)
  R = @inferred fixed_ring(L)
  @test R === base_ring(base_ring(L))
  @test is_maximal(R)

  L = root_lattice(:E, 8)
  R = @inferred fixed_ring(L)
  @test R == ZZ
  @test R != base_ring(base_ring(L))
end

@testset "Misc" begin
  Qx, x = FlintQQ["x"]
  K, a = NumberField(x - 1, "a")
  Kt, t = K["t"]
  E, b = NumberField(t^2 - 2, "b")
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  u = @inferred Hecke._non_norm_rep(E, K, p)
  @test parent(u) === K
  @test @inferred !is_local_norm(E, u, p)
  @test valuation(u - 1, p) ==  normic_defect(E, u, p)
end

@testset "Jordan decomposition" begin
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [3, 2, 1, 2, 3, 1, 1, 1, 1]);
  gens = [[1, -1, 0], [1, -1, 0], [0, 1, -1], [0, 1, -1]]
  L = quadratic_lattice(K, gens, gram = D)
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  B, B, S = jordan_decomposition(L, p)
  @test length(S) == 1
end


text:
# implements Tries

# `shiftby` is equal to the number of bits required to represent index information
# for one level of the BitmappedTrie.
#
# Here, `shiftby` is 5, which means that the BitmappedTrie Arrays will be length 32.
const shiftby = 5
const trielen = 2^shiftby

abstract type BitmappedTrie{T} end

# Copy elements from one Array to another, up to `n` elements.
#
function copy_to(from::Array{T}, to::Array{T}, n::Int) where {T}
    for i=1:n
        to[i] = from[i]
    end
    to
end

# Copies elements from one Array to another of size `len`.
#
copy_to_len(from::Array{T}, len::Int) where {T} =
    copy_to(from, Array{T,1}(undef, len), min(len, length(from)))

mask(t::BitmappedTrie, i::Int) = (((i - 1) >>> shift(t)) & (trielen - 1)) + 1

Base.lastindex(t::BitmappedTrie) = length(t)

Base.length(t::BitmappedTrie) =
    error("$(typeof(t)) does not implement Base.length")
shift(t::BitmappedTrie) =
    error("$(typeof(t)) does not implement FunctionalCollections.shift")
maxlength(t::BitmappedTrie) =
    error("$(typeof(t)) does not implement FunctionalCollections.maxlength")
arrayof(t::BitmappedTrie) =
    error("$(typeof(t)) does not implement FunctionalCollections.arrayof")

function Base.isequal(t1::BitmappedTrie, t2::BitmappedTrie)
    length(t1)    == length(t2)    &&
    shift(t1)     == shift(t2)     &&
    maxlength(t1) == maxlength(t2) &&
    arrayof(t1)   == arrayof(t2)
end

==(t1::BitmappedTrie, t2::BitmappedTrie) = isequal(t1, t2)

# Dense Bitmapped Tries
# =====================

abstract type DenseBitmappedTrie{T} <: BitmappedTrie{T} end

# Why is the shift value of a DenseLeaf 5 instead of 0, and why does
# the shift value of a DenseNode start at 10?
#
# The PersistentVector implements a "tail" optimization, where it
# inserts appended elements into a tail array until that array is 32
# elements long, only then inserting it into the actual bitmapped
# vector trie. This significantly increases the performance of
# operations that touch the very end of the vector (last, append, pop,
# etc.) because you don't have to traverse the trie.
#
# However, it adds a small amount of complexity to the implementation;
# when you query the trie, you now get back a length-32 array instead
# of the actual element. This is why the DenseLeaf has a shift value
# of 5: it leaves an "extra" 5 bits for the PersistentVector to use to
# index into the array returned from the trie. (This also means that a
# DenseNode has to start at shiftby*2.)

struct DenseNode{T} <: DenseBitmappedTrie{T}
    arr::Vector{DenseBitmappedTrie{T}}
    shift::Int
    length::Int
    maxlength::Int
end

struct DenseLeaf{T} <: DenseBitmappedTrie{T}
    arr::Vector{T}
end
DenseLeaf{T}() where {T} = DenseLeaf{T}(T[])

arrayof(    node::DenseNode) = node.arr
shift(      node::DenseNode) = node.shift
maxlength(  node::DenseNode) = node.maxlength
Base.length(node::DenseNode) = node.length

arrayof(    leaf::DenseLeaf) = leaf.arr
shift(          ::DenseLeaf) = shiftby
maxlength(  leaf::DenseLeaf) = trielen
Base.length(leaf::DenseLeaf) = length(arrayof(leaf))

function promoted(node::DenseBitmappedTrie{T}) where T
    DenseNode{T}(DenseBitmappedTrie{T}[node],
                 shift(node) + shiftby,
                 length(node),
                 maxlength(node) * trielen)
end

function demoted(node::DenseNode{T}) where T
    if shift(node) == shiftby * 2
        DenseLeaf{T}(T[])
    else
        DenseNode{T}(DenseBitmappedTrie{T}[],
                     shift(node) - shiftby,
                     0,
                     round(Int, maxlength(node) / trielen))
    end
end

function witharr(node::DenseNode{T}, arr::Array, lenshift::Int=0) where T
    DenseNode{T}(arr, shift(node), length(node) + lenshift, maxlength(node))
end
witharr(leaf::DenseLeaf{T}, arr::Array) where {T} = DenseLeaf{T}(arr)

function append(leaf::DenseLeaf, el)
    if length(leaf) < maxlength(leaf)
        newarr = copy_to_len(arrayof(leaf), 1 + length(leaf))
        newarr[end] = el
        witharr(leaf, newarr)
    else
        append(promoted(leaf), el)
    end
end
function append(node::DenseNode{T}, el) where T
    if length(node) == 0
        child = append(demoted(node), el)
        witharr(node, DenseBitmappedTrie{T}[child], 1)
    elseif length(node) < maxlength(node)
        if length(arrayof(node)[end]) == maxlength(arrayof(node)[end])
            newarr = copy_to_len(arrayof(node), 1 + length(arrayof(node)))
            newarr[end] = append(demoted(node), el)
            witharr(node, newarr, 1)
        else
            newarr = arrayof(node)[:]
            newarr[end] = append(newarr[end], el)
            witharr(node, newarr, 1)
        end
    else
        append(promoted(node), el)
    end
end
push(leaf::DenseLeaf, el) = append(leaf, el)
push(node::DenseNode, el) = append(node, el)

Base.getindex(leaf::DenseLeaf, i::Int) = arrayof(leaf)[mask(leaf, i)]
Base.getindex(node::DenseNode, i::Int) = arrayof(node)[mask(node, i)][i]

function assoc(leaf::DenseLeaf{T}, i::Int, el) where T
    newarr = arrayof(leaf)[:]
    newarr[mask(leaf, i)] = el
    DenseLeaf{T}(newarr)
end
function assoc(node::DenseNode, i::Int, el)
    newarr = arrayof(node)[:]
    idx = mask(node, i)
    newarr[idx] = assoc(newarr[idx], i, el)
    witharr(node, newarr)
end

peek(bt::DenseBitmappedTrie) = bt[end]

# Pop is usually destructive, but that doesn't make sense for an immutable
# structure, so `pop` is defined to return a Trie without its last
# element. Use `peek` to access the last element.
#
pop(leaf::DenseLeaf) = witharr(leaf, arrayof(leaf)[1:end-1])
function pop(node::DenseNode)
    newarr = arrayof(node)[:]
    newarr[end] = pop(newarr[end])
    witharr(node, newarr, -1)
end

# Sparse Bitmapped Tries
# ======================

abstract type SparseBitmappedTrie{T} <: BitmappedTrie{T} end

struct SparseNode{T} <: SparseBitmappedTrie{T}
    arr::Vector{SparseBitmappedTrie{T}}
    shift::Int
    length::Int
    maxlength::Int
    bitmap::Int
end
SparseNode(T::Type) = SparseNode{T}(SparseBitmappedTrie{T}[], shiftby*7, 0, trielen^7, 0)

struct SparseLeaf{T} <: SparseBitmappedTrie{T}
    arr::Vector{T}
    bitmap::Int
end
SparseLeaf{T}() where {T} = SparseLeaf{T}(T[], 0)

arrayof(    n::SparseNode) = n.arr
shift(      n::SparseNode) = n.shift
maxlength(  n::SparseNode) = n.maxlength
Base.length(n::SparseNode) = n.length

arrayof(    l::SparseLeaf) = l.arr
shift(       ::SparseLeaf) = 0
maxlength(  l::SparseLeaf) = trielen
Base.length(l::SparseLeaf) = length(arrayof(l))

function demoted(n::SparseNode{T}) where T
    shift(n) == shiftby ?
    SparseLeaf{T}(T[], 0) :
    SparseNode{T}(SparseBitmappedTrie{T}[],
                  shift(n) - shiftby,
                  0,
                  round(Int, maxlength(n) / trielen), 0)
end

bitpos(  t::SparseBitmappedTrie, i::Int) = 1 << (mask(t, i) - 1)
hasindex(t::SparseBitmappedTrie, i::Int) = t.bitmap & bitpos(t, i) != 0
index(   t::SparseBitmappedTrie, i::Int) =
    1 + count_ones(t.bitmap & (bitpos(t, i) - 1))

function update(l::SparseLeaf{T}, i::Int, el::T) where T
    hasi = hasindex(l, i)
    bitmap = bitpos(l, i) | l.bitmap
    idx = index(l, i)
    if hasi
        newarr = arrayof(l)[:]
        newarr[idx] = el
    else
        newarr = vcat(arrayof(l)[1:idx-1], [el], arrayof(l)[idx:end])
    end
    (SparseLeaf{T}(newarr, bitmap), !hasi)
end
function update(n::SparseNode{T}, i::Int, el::T) where T
    bitmap = bitpos(n, i) | n.bitmap
    idx = index(n, i)
    if hasindex(n, i)
        newarr = arrayof(n)[:]
        updated, inc = update(newarr[idx], i, el)
        newarr[idx] = updated
    else
        child, inc = update(demoted(n), i, el)
        newarr = vcat(arrayof(n)[1:idx-1], [child], arrayof(n)[idx:end])
    end
    (SparseNode{T}(newarr,
                   n.shift,
                   inc ? n.length + 1 : n.length,
                   n.maxlength, bitmap),
     inc)
end

Base.get(n::SparseLeaf, i::Int, default) =
    hasindex(n, i) ? arrayof(n)[index(n, i)] : default
Base.get(n::SparseNode, i::Int, default) =
    hasindex(n, i) ? get(arrayof(n)[index(n, i)], i, default) : default

function initial_state(t::SparseBitmappedTrie)
    t.length == 0 && return Int[]
    ones(Int, 1 + round(Int, t.shift / shiftby))
end

function Base.iterate(t::SparseBitmappedTrie, state = initial_state(t))
    if isempty(state)
        return nothing
    else
        item = directindex(t, state)
        while true
            index = pop!(state)
            node = directindex(t, state)
            if length(node) > index
                push!(state, index + 1)
                return item, vcat(state, ones(Int, 1 + round(Int, t.shift / shiftby) -
                                                   length(state)))
            elseif node === arrayof(t)
                return item, Int[]
            end
        end
    end
end

function directindex(t::SparseBitmappedTrie, v::Vector{Int})
    isempty(v) && return arrayof(t)
    local node = arrayof(t)
    for i=v
        node = node[i]
        node = isa(node, SparseBitmappedTrie) ? arrayof(node) : node
    end
    node
end


text:
function outer_indices!(u)
    uu = Array(u)
    idx = outer_indices!(uu)
    return cu(idx)
end

function outer_indices!(u::Vector)
    M_out  = length(u)
    usum   = zero(eltype(u))

    # compute cumulative sum, overwriting u
    @inbounds for i in 1:M_out
        usum += u[i]
        u[i] = usum
    end

    # shift cumulative sums to the right by one
    @inbounds for i in M_out:-1:2
        u[i] = u[i-1]
    end
    u[1] = 0

    # sample descending sequence of sorted random numbers
    # Algorithm by:
    # Bentley & Saxe, ACM Transactions on Mathematical Software, Vol 6, No 3
    # September 1980, Pages 359--364
    CurMax = one(eltype(u))
    idx = zeros(Int, M_out)
    bindex = M_out # bin index
    @inbounds for i in M_out:-1:1
        CurMax *= exp(log(rand(eltype(u))) / i)
        # scale random numbers (this is equivalent to normalizing u)
        rsample = CurMax * usum
        # checking bindex >= 1 is redundant since
        # ucum[1] = 0
        while rsample < u[bindex]
            bindex -= 1
        end
        idx[i] = bindex
    end
    return idx
end

function outer_resample!(state::BinomialState, model::BinomialGridModel, u)
    idx = outer_indices!(u)
    state.n .= state.n[idx,:]
    state.k .= state.k[idx,:]
    model.Nind .= model.Nind[idx]
    model.pind .= model.pind[idx]
    model.qind .= model.qind[idx]
    model.σind .= model.σind[idx]
    model.τind .= model.τind[idx]
    return state, model
end


text:
@doc raw"""
    block_matrix_logdet(A, B, C, D)

Log-determinant of a block matrix using the determinant lemma.

```math
\det\left(
    \begin{bmatrix}
        \mathbf{A} & \mathbf{B} \\
        \mathbf{C} & \mathbf{D}
    \end{bmatrix}
\right)
= \det(A) \det(D - CA^{-1}B)
= \det(D) \det(A - BD^{-1}C)
```

Here we assume that `A` and `D` are invertible, and moreover are easy to invert
(for example, if they are diagonal).
We use this to chose one or the other of the two formulas above.
"""
function block_matrix_logdet(
    A::AbstractMatrix, B::AbstractMatrix,
    C::AbstractMatrix, D::AbstractMatrix
)
    @assert size(A, 1) == size(B, 1)
    @assert size(C, 1) == size(D, 1)
    @assert size(A, 2) == size(C, 2)
    @assert size(B, 2) == size(D, 2)

    if length(A) ≥ length(D)
        return logdet(A) + logdet(D - C * inv(A) * B)
    else
        return logdet(D) + logdet(A - B * inv(D) * C)
    end
end

@doc raw"""
    block_matrix_invert(A, B, C, D)

Inversion of a block matrix, using the formula:

```math
\begin{bmatrix}
    \mathbf{A} & \mathbf{B} \\
    \mathbf{C} & \mathbf{D}
\end{bmatrix}^{-1}
=
\begin{bmatrix}
    \left(\mathbf{A} - \mathbf{B} \mathbf{D}^{-1} \mathbf{C}\right)^{-1} & \mathbf{0} \\
    \mathbf{0} & \left(\mathbf{D} - \mathbf{C} \mathbf{A}^{-1} \mathbf{B}\right)^{-1}
\end{bmatrix}
\begin{bmatrix}
    \mathbf{I} & -\mathbf{B} \mathbf{D}^{-1} \\
    -\mathbf{C} \mathbf{A}^{-1} & \mathbf{I}
\end{bmatrix}
```

Assumes that `A` and `D` are square and invertible.
"""
function block_matrix_invert(
    A::AbstractMatrix, B::AbstractMatrix,
    C::AbstractMatrix, D::AbstractMatrix
)
    @assert size(A, 1) == size(A, 2) == size(B, 1) == size(C, 2)
    @assert size(D, 1) == size(D, 2) == size(B, 2) == size(C, 1)

    CinvA = C * inv(A)
    BinvD = B * inv(D)

    M = [
        inv(A - BinvD * C)  zero(B)
        zero(C)  inv(D - CinvA * B)
    ]

    N = [
        I(size(B, 1))  -BinvD
        -CinvA  I(size(C, 1))
    ]

    return M * N
end


text:
"""
    sandwich!(m::WSVarLmmModel)

Calculate the sandwich estimator of the asymptotic covariance of the parameters,
based on values `m.Hββ`, `m.Hττ`, `m.HτLγ`, `m.HLγLγ`, `m.data[i].∇β`,
`m.data[i].∇τ`, and `m.data[i].∇Lγ`. `m.vcov` is updated by the sandwich
estimator and returned.
"""
function sandwich!(m::WSVarLmmModel{T}) where T <: BlasReal
    p, q, l = m.p, m.q, m.l
    minv    = inv(m.m)
    # form A matrix in the sandwich formula
    fill!(m.Ainv, 0)
    m.Ainv[          1:p,                 1:p      ] = m.Hββ
    m.Ainv[    (p + 1):(p + l),     (p + 1):(p + l)] = m.Hττ
    m.Ainv[    (p + 1):(p + l), (p + l + 1):end    ] = m.HτLγ
    m.Ainv[(p + l + 1):end,     (p + l + 1):end    ] = m.HLγLγ
    copytri!(m.Ainv, 'U')
    lmul!(minv, m.Ainv)
    # form B matrix in the sandwich formula
    fill!(m.B, 0)
    for obs in m.data
        copyto!(m.ψ, 1    , obs.∇β)
        copyto!(m.ψ, p + 1, obs.∇τ)
        offset = p + l + 1
        @inbounds for j in 1:q, i in j:q
            m.ψ[offset] = obs.∇Lγ[i, j]
            offset += 1
        end
        BLAS.syr!('U', T(1), m.ψ, m.B)
    end
    copytri!(m.B, 'U')
    lmul!(minv, m.B)
    # calculuate A (pseudo)-inverse
    Aeval, Aevec = eigen(Symmetric(m.Ainv))
    @inbounds for j in 1:size(m.Ainv, 2)
        ej = Aeval[j]
        invsqrtej = ej > sqrt(eps(T)) ? inv(sqrt(ej)) : T(0)
        for i in 1:size(m.Ainv, 1)
            Aevec[i, j] *= invsqrtej
        end
    end
    mul!(m.Ainv, Aevec, transpose(Aevec))
    # calculate vcov
    mul!(Aevec , m.Ainv, m.B   ) # use Avec as scratch space
    mul!(m.vcov,  Aevec, m.Ainv)
    m.vcov .*= minv
end


text:
"""
	$(SIGNATURES)

Output of a Cobb-Douglas. Returns a vector of length 1 for xM with 1 row.
"""
function output_cobb_douglas(xM :: AbstractMatrix{F1},
    expV :: AbstractVector{F1}) where F1 <: Real

    T, n = size(xM);
    outputV = ones(F1, T);
    for j = 1 : n
        @views outputV .*= (xM[:,j] .^ expV[j]);
    end
    return outputV
end

# Returns scalar output.
function output_cobb_douglas(xV :: AbstractVector{F1},
    expV :: AbstractVector{F1}) where F1 <: Real

    n = length(xV);
    output = one(F1);
    for j = 1 : n
        output *= xV[j] ^ expV[j];
    end
    return output
end

function mproducts_cobb_douglas(xM :: AbstractMatrix{F1},
    expV :: AbstractVector{F1}) where F1 <: Real

    yV = output_cobb_douglas(xM, expV);
    mpM = (yV * expV') ./ xM;
    return mpM
end

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


text:
# BBI B+-Tree
# ===========

const BTREE_MAGIC = 0x78CA8C91

# Supplemental Table 8.
struct BTreeHeader
    magic::UInt32
    block_size::UInt32
    key_size::UInt32
    val_size::UInt32
    item_count::UInt64
    reserved::UInt64
end

function Base.write(stream::IO, header::BTreeHeader)
    return write(
        stream,
        header.magic,    header.block_size, header.key_size,
        header.val_size, header.item_count, header.reserved)
end

# Supplemental Table 9.
struct BTreeNode
    isleaf::UInt8
    reserved::UInt8
    count::UInt16
end

function Base.write(stream::IO, node::BTreeNode)
    return write(stream, node.isleaf, node.reserved, node.count)
end

function isleaf(node::BTreeNode)
    return node.isleaf != 0
end

function Base.read(io::IO, ::Type{BTreeNode})
    return BTreeNode(read(io, UInt8), read(io, UInt8), read(io, UInt16))
end

# disk-serialized B+ tree
struct BTree{T<:IO}
    stream::T
    offset::UInt64
    header::BTreeHeader
end

function BTree(stream::IO, offset::Integer)
    read32() = read(stream, UInt32)
    read64() = read(stream, UInt64)
    seek(stream, offset)
    magic = read32()
    if magic != BTREE_MAGIC
        error("invalid B+-tree magic bytes")
    end
    return BTree(stream, UInt64(offset), BTreeHeader(magic, read32(), read32(), read32(), read64(), read64()))
end

# Load the chromosome list.
function chromlist(tree::BTree)
    list = Tuple{String,UInt32,UInt32}[]
    # traverse tree
    key = Vector{UInt8}(undef, tree.header.key_size)
    stack = UInt64[tree.offset + sizeof(BTreeHeader)]
    while !isempty(stack)
        offset = pop!(stack)
        seek(tree.stream, offset)
        node = read(tree.stream, BTreeNode)
        for i in 1:node.count
            if isleaf(node)
                read!(tree.stream, key)
                i = findfirst(isequal(0x00), key)
                chromname = i === nothing ? String(copy(key)) : String(key[1:i-1])
                chromid = read(tree.stream, UInt32)
                chromsize = read(tree.stream, UInt32)
                push!(list, (chromname, chromid, chromsize))
            else
                skip(tree.stream, tree.header.key_size)
                offset = read(tree.stream, UInt64)
                push!(stack, offset)
            end
        end
    end
    return list
end

function write_btree(stream::IO, chromlist::Vector{Tuple{String,UInt32,UInt32}})
    # This function stores all chromosomes in the root node as a leaf because it
    # is simple to implement.
    blksize = length(chromlist)
    keysize = Base.maximum(sizeof(name) for (name, _) in chromlist)
    valsize = 8
    n = 0

    # write header
    n += write(stream, BTreeHeader(BTREE_MAGIC, blksize, keysize, valsize, length(chromlist), 0))

    # write the root node format
    n += write(stream, BTreeNode(0x01, 0x00, length(chromlist)))

    # write the root node
    key = Vector{UInt8}(undef, keysize)
    for (name, id, len) in sort(chromlist, by=x->x[1])  # sort by name
        fill!(key, 0x00)
        @assert sizeof(name) ≤ keysize
        Mem.copy(key, name, sizeof(name))
        n += write(stream, key, id, len)
    end

    return n
end

# Add a unique ID for each chromosome; chromlist is a tuple of (chrom name, chrom length).
function add_chrom_ids(chromlist::Vector{Tuple{String,UInt32}})
    chromlist = sort(chromlist, by=first)  # sort by chromosome name
    return [(name, UInt32(id - 1), UInt32(len)) for (id, (name, len)) in enumerate(chromlist)]
end

function add_chrom_ids(chromlist)
    return add_chrom_ids([(String(name), UInt32(len)) for (name, len) in chromlist])
end


text:
module FESetShellT3Module

using FinEtools
using FinEtools.MeshQuadrilateralModule
import FinEtools.FESetModule: cat, subset, nodesperelem
using FinEtools.MatrixUtilityModule: complete_lt!
using LinearAlgebra: norm, Transpose, mul!

struct FESetShellT3 <: AbstractFESet2Manifold{3}
end

end # module


text:
# Based on code written by Ole Kröger (@Wikunia).
# https://github.com/Wikunia/ConstraintSolver.jl/blob/0f63f1601f007852d6cc5939ab3768dd4c4b237c/src/MOI_wrapper/Bridges/strictly_greater_than.jl

struct StrictlyGreaterThan2StrictlyLessThanBridge{T <: Real} <: MOIBC.AbstractBridge
    con::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, CP.Strictly{MOI.LessThan{T}, T}}
end

function MOIBC.bridge_constraint(
    ::Type{StrictlyGreaterThan2StrictlyLessThanBridge{T}},
    model,
    f::MOI.VariableIndex,
    s::CP.Strictly{MOI.GreaterThan{T}, T},
) where {T}
    return MOIBC.bridge_constraint(
        StrictlyGreaterThan2StrictlyLessThanBridge{T},
        model,
        MOI.ScalarAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{StrictlyGreaterThan2StrictlyLessThanBridge{T}},
    model,
    f::MOI.ScalarAffineFunction{T},
    s::CP.Strictly{MOI.GreaterThan{T}, T},
) where {T}
    con = MOI.add_constraint(
        model,
        -f,
        CP.Strictly(MOI.LessThan(- s.set.lower))
    )

    return StrictlyGreaterThan2StrictlyLessThanBridge(con)
end

function MOI.supports_constraint(
    ::Type{StrictlyGreaterThan2StrictlyLessThanBridge{T}},
    ::Union{Type{MOI.VariableIndex}, Type{MOI.ScalarAffineFunction{T}}},
    ::Type{CP.Strictly{MOI.GreaterThan{T}, T}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{StrictlyGreaterThan2StrictlyLessThanBridge{T}}) where {T}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{StrictlyGreaterThan2StrictlyLessThanBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, CP.Strictly{MOI.LessThan{T}, T}),
    ]
end

function MOI.get(
    ::StrictlyGreaterThan2StrictlyLessThanBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, CP.Strictly{MOI.LessThan{T}, T},
    },
) where {T}
    return 1
end

function MOI.get(
    b::StrictlyGreaterThan2StrictlyLessThanBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, CP.Strictly{MOI.LessThan{T}, T},
    },
) where {T}
    return [b.con]
end