text:
@polly function kernel_correlation(float_n, data, corr, mean, stddev)
    n,m = size(data)
    eps = 0.1;

    for j = 1:m
        mean[j] = zero(eltype(mean))
        for i = 1:n
            mean[j] += data[i,j]
        end
        mean[j] /= float_n
    end

    for j = 1:m
        stddev[j] = zero(eltype(stddev))
        for i = 1:n
            stddev[j] += (data[i,j] - mean[j]) * (data[i,j] - mean[j])
        end
        stddev[j] /= float_n
        @fastmath stddev[j] = sqrt(stddev[j])
        stddev[j] = stddev[j] <= eps ? one(eltype(stddev)) : stddev[j]
    end

    # Center and reduce the column vectors.
    for i = 1:n, j = 1:m
        data[i,j] -= mean[j]
        @fastmath data[i,j] /= sqrt(float_n) * stddev[j]
    end

    # Calculate the m * m correlation matrix.
    for i = 1:(m-1)
        corr[i,i] = one(eltype(corr))
        for j = (i+1):m
            corr[i,j] = one(eltype(corr))
            for k = 1:n
                corr[i,j] += (data[k,i] * data[k,j])
            end
            corr[j,i] = corr[i,j]
        end
    end

    corr[m-1,m-1] = one(eltype(corr))
end

let
    m = 1200
    n = 1400

    data = zeros(Float32,n,m)
    corr = zeros(Float32,m,m)
    mean = zeros(Float32,m)
    stddev = zeros(Float32,m)

    float_n = Float32(n);

    for i = 1:n, j = 1:m
        data[i,j] = (i*j)/m + i
    end

    SUITE["correlation"] = @benchmarkable kernel_correlation($float_n, data, corr, mean, stddev) setup=(data=copy($data); corr=copy($corr); mean=copy($mean); stddev=copy($stddev))
end


text:
using TensorToolbox
using Random
using LinearAlgebra
Random.seed!(3)
include("LD.jl")

function main()
    T = rand( 100,100, 100)
    R = TLD(T)

    T = rand( 100,100, 100)
    R_newton = TLD_newton(T)

    println(" reconst_error by grad decent ",  norm(T - R) )
    println(" reconst_error by newton method ", norm(T - R_newton) )
end

main()

ext: jl

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
using Printf

"""
    LATTICE

    Defines a lattice type. Supported lattice types include 'primitive',
    'body_centered' and 'face_centered' lattices.
"""
module LATTICE
    @enum TYPE begin
        primitive       = 1
        body_centered   = 2
        face_centered   = 3
    end
end

"""
    generate_template(cell_sizes::Vector{Float64}, type::LATTICE.TYPE)

        Generate a template for the given lattice 'type', with the cell size
        defined in 'cell_sizes' (should be a 3x1 Vector{Int64, 1}, with cell
        sizes in each of the 3 dimensions, in nm). Supported lattice types
        include 'primitive', 'body_centered' and 'face_centered' lattices.
        Returns an Array{Float64, 2}.

    Example:

        generate_template([1.0, 1.0, 1.0], LATTICE.primitive)
"""
function generate_template(cell_sizes::Vector{Float64}, type::LATTICE.TYPE)::Array{Float64, 2}
    hf = cell_sizes ./ 2
    atoms = [0.0 0.0 0.0]
    if type == LATTICE.body_centered
        atoms = vcat(atoms, [hf[1] hf[2] hf[3]])
    elseif type == LATTICE.face_centered
        atoms = vcat(atoms, [hf[1] hf[2] 0.0; hf[1] 0.0 hf[3]; 0.0 hf[2] hf[3]])
    end
    return atoms
end

"""
    calculate_n_atoms(rep::Vector{Int64}, type::LATTICE.TYPE)

        Calculate the number of atoms in a lattice of a given 'type', for the
        given amount of repetitions 'rep' (should be a 3x1 Vector{Int64, 1},
        with repetitions in each of the 3 dimensions). Supported lattice types
        include 'primitive', 'body_centered' and 'face_centered' lattices.
        Returns an Int64.

    Example:

        calculate_n_atoms([10, 10, 10], LATTICE.primitive)
"""
function calculate_n_atoms(rep::Vector{Int64}, type::LATTICE.TYPE)::Int64

    @assert size(rep) == (3,) "rep argument must be a 3x1 Array{Int64, 1}"

    if type == LATTICE.primitive
        n_atoms = (rep[1] + 1) * (rep[2] + 1) * (rep[3] + 1)
    elseif type == LATTICE.body_centered
        n_atoms = (rep[1] + 1) * (rep[2] + 1) * (rep[3] + 1)
        n_atoms += rep[1]*rep[2]*rep[3]
    elseif type == LATTICE.face_centered
        n_atoms = (rep[1] + 1) * (rep[2] + 1) * (rep[3] + 1)
        n_atoms += rep[1]*rep[2]*rep[3] * 3
        n_atoms += rep[1]*rep[2] + rep[2] * rep[3] + rep[3] * rep[1]
    end
    return n_atoms
end

"""
    generate_cubic_lattice(sizes::Float64, rep::Int64, type::LATTICE.TYPE)

        Generates a cubic lattice with all dimensions lenght set to 'sizes' (in
        nm). Each dimension contains 'rep' repeats of the particle template, of
        a specific 'type'. Supported lattice types include 'primitive',
        'body_centered' and 'face_centered' lattices.

        Returns a Matrix{Float64} in Array of Structures (AoS) format.
        Use 'convert' function to change the content type.
        Ex: convert(Matrix{Float32}, array)
        Use 'collect' function with the transpose to change the format to SoA.
        Ex: collect(array')

    Example:

        generate_cubic_lattice(10.0, 10, LATTICE.primitive)
"""
function generate_cubic_lattice(sizes::Float64, rep::Int64, type::LATTICE.TYPE)::Matrix{Float64}

    xyz = generate_lattice([sizes, sizes, sizes], [rep, rep, rep], type)
    return xyz
end

"""
    generate_lattice(sizes::Vector{Float64, 1}, rep::Vector{Int64, 1}, type::LATTICE.TYPE)

        Generates a lattice with dimensions lenght set to 'sizes' (should be a
        3x1 array with each dimanesion size, in nm). Each dimension contains
        'rep' repeats of the particle template of a specific 'type' (should be a
        3x1 array with each dimension repetitions). Supported lattice types
        include 'primitive', 'body_centered' and 'face_centered' lattices.

        Returns a Matrix{Float64} in Array of Structures (AoS) format. Use
        'convert function to change the content type.
        Ex: convert(Matrix{Float32}, array)

    Example:

        generate_lattice([10.0, 10.0, 5.0], [10, 10, 5], LATTICE.primitive)
"""
function generate_lattice(sizes::Vector{Float64}, rep::Vector{Int64},
    type::LATTICE.TYPE)::Matrix{Float64}

    @assert size(sizes) == (3,) "sizes argument must be a 3x1 Array{Float64, 1}"
    @assert size(rep) == (3,) "rep argument must be a 3x1 Array{Int64, 1}"

    # Generate the template
    cs = (sizes./rep)
    template = generate_template(cs, type)

    # Create empty list of atoms
    n_atoms  = calculate_n_atoms(rep, type)
    xyz = zeros(3, n_atoms)

    # Fill xyz with repetitions of template
    c = 1
    for i in 0:rep[1]
        for j in 0:rep[2]
            for k in 0:rep[3]
                for index in 1:size(template, 1)
                    new_pos = template[index, :] .+ [i*cs[1], j*cs[2], k*cs[3]]
                    if all(new_pos .<= sizes)
                        xyz[1:3, c] .= new_pos
                        c +=1
                    end
                end
            end
        end
    end

    return xyz
end

"""
    print_xyz(coords::Matrix{T}, [io::Union{String, Base.TTY} = Base.stdout]) where {T <: AbstractFloat}

        Print the contents of a Matrix{T}, where T is a Float and the format can
        be either AoS or SoA, in a XYZ readable format. By default, frints to
        stdout, but a file name can be given in argumento 'io' to print to file
        instead.

    Example:

        print_xyz(coords)
        print_xyz(coords, "filename.xyz")
"""
function print_xyz(coords::Matrix{T}, io::Union{String, Base.TTY} = Base.stdout) where {T <: AbstractFloat}

    @assert min(size(coords)...) == 3 "Minimum size of matrix must be 3."
    AoS = false
    if argmin(collect(size(coords))) == 1
        AoS = true
    end

    if !isa(io, Base.TTY)
        io = open(io, "w")
    end

    for i in 1:max(size(coords)...)
        if AoS
            xyz = coords[:, i]
        else
            xyz = coords[i, :]
        end
        println(io, @sprintf "%-5s %8.2f %8.2f %8.2f" "H" xyz[1] xyz[2] xyz[3])
    end

    if !isa(io, Base.TTY)
        println(" Printed coordinates to $(io.name)")
        close(io)
    end
end


text:
using Clang.cindex
using Base.Test

top = cindex.parse("cxx/cbasic.h")

###############################################################################
# Test tokenizer interface

CONST = tokenize(cindex.search(
                    top,
                    x-> (isa(x, cindex.MacroDefinition) && name(x) == "CONST"))[1])
CONSTADD = tokenize(cindex.search(
                    top,
                    x-> (isa(x, cindex.MacroDefinition) && name(x) == "CONSTADD"))[1])

@test (typeof(CONST[1]) == cindex.Identifier && typeof(CONST[2]) == cindex.Literal)
@test (CONSTADD[1].text == "CONSTADD" && CONSTADD[2].text == "CONST" &&
       CONSTADD[3].text == "+" && CONSTADD[4].text == "2")


text:
import Base: *, eltype, size
import LinearAlgebra: adjoint, mul!, rmul!

struct Householder{T,S<:StridedVector}
    v::S
    τ::T
end
struct HouseholderBlock{T,S<:StridedMatrix,U<:StridedMatrix}
    V::S
    T::UpperTriangular{T,U}
end

size(H::Householder) = (length(H.v), length(H.v))
size(H::Householder, i::Integer) = i <= 2 ? length(H.v) : 1

eltype(H::Householder{T})      where T = T
eltype(H::HouseholderBlock{T}) where T = T

adjoint(H::Householder{T})      where {T} = Adjoint{T,typeof(H)}(H)
adjoint(H::HouseholderBlock{T}) where {T} = Adjoint{T,typeof(H)}(H)

function lmul!(H::Householder, A::StridedMatrix)
    m, n = size(A)
    length(H.v) == m - 1 || throw(DimensionMismatch(""))
    v = view(H.v, 1:m - 1)
    τ = H.τ
    for j = 1:n
        va = A[1,j]
        Aj = view(A, 2:m, j)
        va += dot(v, Aj)
        va = τ*va
        A[1,j] -= va
        axpy!(-va, v, Aj)
    end
    A
end

function rmul!(A::StridedMatrix, H::Householder)
    m, n = size(A)
    length(H.v) == n - 1 || throw(DimensionMismatch(""))
    v = view(H.v, :)
    τ = H.τ
    a1 = view(A, :, 1)
    A1 = view(A, :, 2:n)
    x = A1*v
    axpy!(one(τ), a1, x)
    axpy!(-τ, x, a1)
    rankUpdate!(A1, x, v, -τ)
    A
end

function lmul!(adjH::Adjoint{<:Any,<:Householder}, A::StridedMatrix)
    H = parent(adjH)
    m, n = size(A)
    length(H.v) == m - 1 || throw(DimensionMismatch(""))
    v = view(H.v, 1:m - 1)
    τ = H.τ
    for j = 1:n
        va = A[1,j]
        Aj = view(A, 2:m, j)
        va += dot(v, Aj)
        va = τ'va
        A[1,j] -= va
        axpy!(-va, v, Aj)
    end
    A
end

# FixMe! This is a weird multiplication method and at least its behavior needs to be explained
function lmul!(H::HouseholderBlock{T}, A::StridedMatrix{T}, M::StridedMatrix{T}) where T
    V = H.V
    mA, nA = size(A)
    nH = size(V, 1)
    blocksize = min(nH, size(V, 2))
    nH == mA || throw(DimensionMismatch(""))

    # Reflector block is split into a UnitLowerTriangular top part and rectangular lower part
    V1 = LinearAlgebra.UnitLowerTriangular(view(V, 1:blocksize, 1:blocksize))
    V2 = view(V, blocksize+1:mA, 1:blocksize)

    # Split A to match the split in the reflector block
    A1 = view(A, 1:blocksize, 1:nA)
    A2 = view(A, blocksize+1:mA, 1:nA)

    # Copy top part of A (A1) to temporary work space
    copyto!(M, A1)
    # Multiply UnitLowerTriangular V1 and A1 in-place in M (M = V1'A1)
    lmul!(V1', M)
    # Add V2'A2 to V1'A1 (M := M + V2'A2)
    mul!(M, V2', A2, one(T), one(T))
    # Multiply the elementary block loading T (UpperTriangular) and M in-place
    lmul!(H.T, M)

    # A2 := A2 - V2*M
    mul!(A2, V2, M, -one(T), one(T))
    # A1 := A1 - V1*M but since V1 is UnitLowerTriangular we do it in two steps:
    ## 1. M  := -V1*M
    ## 2. A1 := A1 + M
    lmul!(V1, M, -one(T))
    axpy!(one(T), M, A1)

    return A
end
(*)(H::HouseholderBlock{T}, A::StridedMatrix{T}) where {T} =
        lmul!(H, copy(A), similar(A, (min(size(H.V)...), size(A, 2))))

function lmul!(adjH::Adjoint{T,<:HouseholderBlock{T}}, A::StridedMatrix{T}, M::StridedMatrix) where T
    H = parent(adjH)
    V = H.V
    mA, nA = size(A)
    nH = size(V, 1)
    blocksize = min(nH, size(V, 2))
    nH == mA || throw(DimensionMismatch(""))

    # Reflector block is split into a UnitLowerTriangular top part and rectangular lower part
    V1 = LinearAlgebra.UnitLowerTriangular(view(V, 1:blocksize, 1:blocksize))
    V2 = view(V, blocksize+1:mA, 1:blocksize)

    # Split A to match the split in the reflector block
    A1 = view(A, 1:blocksize, 1:nA)
    A2 = view(A, blocksize+1:mA, 1:nA)

    # Copy top part of A (A1) to temporary work space
    copyto!(M, A1)
    # Multiply UnitLowerTriangular V1 and A1 in-place in M (M = V1'A1)
    lmul!(V1', M)
    # Add V2'A2 to V1'A1 (M := M + V2'A2)
    mul!(M, V2', A2, one(T), one(T))
    # Multiply the elementary block loading T (UpperTriangular) and M in-place
    lmul!(H.T', M)

    # A2 := A2 - V2*M
    mul!(A2, V2, M, -one(T), one(T))
    # A1 := A1 - V1*M but since V1 is UnitLowerTriangular we do it in two steps:
    ## 1. M  := -V1*M
    ## 2. A1 := A1 + M
    lmul!(V1, M, -one(T))
    axpy!(one(T), M, A1)

    return A
end
(*)(adjH::Adjoint{T,<:HouseholderBlock{T}}, A::StridedMatrix{T}) where {T} =
        lmul!(adjH, copy(A), similar(A, (min(size(parent(adjH).V)...), size(A, 2))))

convert(::Type{Matrix}, H::Householder{T}) where {T} = lmul!(H, Matrix{T}(I, size(H, 1), size(H, 1)))
convert(::Type{Matrix{T}}, H::Householder{T}) where {T} = lmul!(H, Matrix{T}(I, size(H, 1), size(H, 1)))


text:
using ExploreFunction
using SafeTestsets

@safetestset "MCMCMapping.jl" begin
    include("MCMCMapping_test.jl")
end

@safetestset "VietorisRips.jl" begin
    include("VietorisRips_test.jl")
end


text:
using GLMakie,  FileIO
using Downloads: download
nasa = "https://eoimages.gsfc.nasa.gov/images/imagerecords/"
earth_link = nasa * "73000/73963/gebco_08_rev_bath_3600x1800_color.jpg"
earth = load(download(earth_link))
sun_link = "https://www.solarsystemscope.com/textures/download/2k_sun.jpg"
sun = load(download(sun_link))
n = 1024 ÷ 4 # 2048
U = LinRange(-pi, pi, 2n)
V = LinRange(-pi, pi, n)
x1 = [cos(u) + .5 * cos(u) * cos(v)      for u in U, v in V]
y1 = [sin(u) + .5 * sin(u) * cos(v)      for u in U, v in V]
z1 = [.5 * sin(v)                        for u in U, v in V]
x2 = [1 + cos(u) + .5 * cos(u) * cos(v)  for u in U, v in V]
y2 = [.5 * sin(v)                        for u in U, v in V]
z2 = [sin(u) + .5 * sin(u) * cos(v)      for u in U, v in V]

with_theme(theme_dark()) do
    fig, ax, = surface(x1, y1, z1; color = earth'[:,end:-1:1],
        shading=true, ambient = Vec3f0(0.85, 0.85, 0.85), backlight = 1.5f0,
        axis=(type=Axis3, aspect = :data, perspectiveness = 0.5,
        azimuth = -1.56, elevation = 0.58),
        figure = (;resolution = (1200,800)))
    surface!(ax, x2, y2, z2,  color = sun'[:,end:-1:1], shading= true,
        ambient = Vec3f0(0.85, 0.85, 0.85), backlight = 1.5f0)
    hidedecorations!(ax; grid=false)
    save("surface_tori.png", fig)
    display(fig)
end


text:
module MPSIMD

using SIMD

struct Internal{F,I,N}
end

export IntN
struct IntN{F,I,N} #<: Signed
    # TODO: add a flag memoizing whether the representation is normalized

    digits::Vec{N,F}

    function IntN(::Internal{F1, I1, N1},
                  digits::Vec{N1, F1}) where {F1, I1, N1}
        @assert F1 <: AbstractFloat
        @assert I1 <: Signed
        @assert N1 isa Integer
        @assert N1 > 0
        new{F1, I1, N1}(digits)
    end
    function IntN(::Internal{F1, I1, N1}) where {F1, I1, N1}
        @assert F1 <: AbstractFloat
        @assert I1 <: Signed
        @assert N1 isa Integer
        @assert N1 > 0
        new{F1, I1, N1}(Vec{N1, F1}(0))
    end
end

function IntN{F,I,N}() where {F,I,N}
    IntN(Internal{F,I,N}())
end

function Base.show(io::IO, x::IntN{F,I,N}) where {F,I,N}
    print(io, "($F,$I,$N)[")
    join(io, [string(Integer(x.digits[n]), base=16) for n in N:-1:1], ", ")
    print(io, "]")
end

export normalize
function normalize(x::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    nbits = 8 * sizeof(I)
    carry = F(0)
    xs = x.digits
    rs = zero(typeof(xs))
    for n in 1:N
        rn = xs[n] + carry
        carry = round(rn * ldexp(F(1), -nbits))
        rn = rn - carry * ldexp(F(1), nbits)
        rs = setindex(rs, rn, n)
    end
    @assert carry == 0
    # @assert carry == 0 || carry == -1
    # rs = setindex(rs, rs[N] + carry * ldexp(F(1), nbits), N)
    IntN(Internal{F,I,N}(), rs)
end

function IntN{F,I,N}(x::IntN{F,I,N}) where {F,I,N}
    x
end

function IntN{F,I,N}(i::J) where {J<:Integer, F,I,N}
    nbits = 8 * sizeof(I)
    if Base.hastypemax(J) &&
        typemin(I) <= typemin(J) && typemax(J) <= typemax(I)
        ds = setindex(IntN{F,I,N}().digits, F(i), 1)
        return IntN(Internal{F,I,N}(), ds)
    end
    ds = IntN{F,I,N}().digits
    for n in 1:N
        i0 = i
        i = (i >> nbits) + (i < 0)
        d = F(i0 - (i << nbits))
        ds = setindex(ds, d, n)
    end
    @assert i == 0 || i == -1
    ds = setindex(ds, ds[N] + (i << nbits), N)
    IntN(Internal{F,I,N}(), ds)
end

function (::Type{J})(x::IntN{F,I,N}) where {J<:Integer, F,I,N}
    nbits = 8 * sizeof(I)
    r = J(0)
    for n in N:-1:1
        r <<= nbits
        r += J(x.digits[n])
    end
    r
end

function Base.cmp(x::IntN{F,I,N}, y::IntN{F,I,N})::Int where {F,I,N}
    J = SIMD.int_type(F)
    x1 = normalize(x)
    y1 = normalize(y)
    # r = cmp(x1.digits, y1.digits)
    r = vifelse(x1.digits == y1.digits, Vec{N, J}(0),
                vifelse(x1.digits < y1.digits, Vec{N, J}(-1), Vec{N, J}(1)))
    for n in N:-1:1
        r[n] != 0 && return r[n]
    end
    0
end

function Base. ==(x::IntN{F,I,N}, y::IntN{F,I,N})::Bool where {F,I,N}
    cmp(x, y) == 0
end

function Base.isless(x::IntN{F,I,N}, y::IntN{F,I,N})::Bool where {F,I,N}
    cmp(x, y) < 0
end

function Base.signbit(x::IntN{F,I,N})::Bool where {F,I,N}
    cmp(x, zero(x)) < 0
end

function Base.flipsign(x::IntN{F,I,N},
                       y::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    ifelse(signbit(y), -x, x)
end

function Base.abs(x::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    flipsign(x, x)
end

function Base.copysign(x::IntN{F,I,N},
                       y::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    flipsign(abs(x), y)
end

function Base.zero(::Type{IntN{F,I,N}})::IntN{F,I,N} where {F,I,N}
    IntN{F,I,N}()
end

function Base.one(::Type{IntN{F,I,N}})::IntN{F,I,N} where {F,I,N}
    IntN{F,I,N}(1)
end

function Base.zero(::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    zero(IntN{F,I,N})
end

function Base.one(::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    one(IntN{F,I,N})
end

function shift_up(x::Vec{N,F})::Vec{N,F} where {N,F}
    # @assert x[N] == 0
    shufflevector(x, zero(typeof(x)), Val{ntuple(n -> n==1 ? 2*N-1 : n-2, N)})
end

function shift_up(x::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    IntN(Internal{F,I,N}(), shift_up(x.digits))
end

export normalize1
function normalize1(::Type{I}, x::Vec{N,F}) where {I,N,F}
    nbits = 8 * sizeof(I)
    rhi = round(x * ldexp(F(1), -nbits))
    rlo = x - rhi * ldexp(F(1), nbits)
    rlo + shift_up(rhi)
end

function normalize1(x::IntN{F,I,N}) where {F,I,N}
    IntN(Internal{F,I,N}(), normalize1(I, x.digits))
end

function Base. +(x::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    x
end

function Base. -(x::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    IntN(Internal{F,I,N}(), -x.digits)
end

export fastadd
@inline function fastadd(x::IntN{F,I,N},
                  y::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    IntN(Internal{F,I,N}(), x.digits + y.digits)
end

export fastsub
@inline function fastsub(x::IntN{F,I,N},
                  y::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    IntN(Internal{F,I,N}(), x.digits - y.digits)
end

function Base. +(x::IntN{F,I,N},
                 y::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    normalize1(fastadd(x, y))
end

function Base. -(x::IntN{F,I,N},
                 y::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    normalize1(fastsub(x, y))
end

export fastmul
@inline function fastmul(x::IntN{F,I,N}, y::I)::IntN{F,I,N} where {F,I,N}
    fastmul(x, F(y))
end

@inline function fastmul(x::IntN{F,I,N}, y::F)::IntN{F,I,N} where {F,I,N}
    nbits = 8 * sizeof(I)
    nbits2 = nbits ÷ 2
    yhi = round(y * ldexp(F(1), -nbits2)) * ldexp(F(1), nbits2)
    ylo = y - yhi
    rlo = x.digits * ylo
    rhi = normalize1(I, x.digits * yhi)
    r = rlo + rhi
    IntN(Internal{F,I,N}(), r)
end

@inline function fastmul(x::IntN{F,I,N}, y::IntN{F,I,N}) where {F,I,N}
    r = zero(x)
    for n in N:-1:1
        r = fastadd(shift_up(r), fastmul(x, y.digits[n]))
    end
    r
end

function Base. *(x::IntN{F,I,N}, y::I)::IntN{F,I,N} where {F,I,N}
    normalize1(fastmul(x, y))
end

function Base. *(x::IntN{F,I,N}, y::IntN{F,I,N})::IntN{F,I,N} where {F,I,N}
    normalize1(fastmul(x, y))
end

end


text:
"""
    Tetrahedron3

An object `t` of type [`Tetrahedron3`](@ref) is an oriented tetrahedron in the
three-dimensional Euclidean space ``Ε^3``.

It is defined by four vertices ``p_0``, ``p_1``, ``p_2`` and ``p_3``. The
orientation of a tetrahedron is the orientation of its four vertices. That means
it is positive when ``p_3`` is on the positive side of the plane defined by
``p_0``, ``p_1`` and ``p_2``.

The tetrahedron itself splits the space ``Ε^3`` in a *positive* and a *negative*
side.

The boundary of a tetrahedron splits the space in two open regions, a bounded
one and an unbounded one.
"""
Tetrahedron3

"""
    Tetrahedron3(p₀::Point3, p₁::Point3, p₂::Point3, p₃::Point3)

Introduces a tetrahedron `t` with vertices `p₀`, `p₁`, `p₂` and `p₃`.
"""
Tetrahedron3(p₀::Point3, p₁::Point3, p₂::Point3, p₃::Point3)

"""
    ==(t₁::Tetrahedron3, t₂::Tetrahedron3)

Test for equality: two tetrahedra `t₁` and `t₂` are equal, iff `t₁` and `t₂`
have the same orientation and their sets (not sequences) of vertices are equal.
"""
==(t₁::Tetrahedron3, t₂::Tetrahedron3)

"""
    vertex(t::Tetrahedron3, i::Integer)

Returns the `i`ᵗʰ vertex modulo 4 of `t`.
"""
vertex(t::Tetrahedron3, i::Integer)

"""
    is_degenerate(t::Tetrahedron3)

Tetrahedron `t` is degenerate, if the vertices are coplanar.
"""
is_degenerate(t::Tetrahedron3)

"""
    orientation(t::Tetrahedron3)

Returns the orientation of `t`.
"""
orientation(t::Tetrahedron3)

"""
    oriented_side(t::Tetrahedron3, p::Point2)

Returns [`ON_ORIENTED_BOUNDARY`](@ref), or [`POSITIVE_SIDE`](@ref), or the
constant [`ON_NEGATIVE_SIDE`](@ref), determined by the position of point `p`.

!!! info "Precondition"

    `t` is not degenerate.
"""
oriented_side(t::Tetrahedron3, p::Point2)

"""
    bounded_side(t::Tetrahedron3, p::Point2)

Returns [`ON_BOUNDARY`](@ref), [`ON_BOUNDED_SIDE`](@ref), or else
[`ON_UNBOUNDED_SIDE`](@ref), depending on where point `p` is.

!!! info "Precondition"

    `t` is not degenerate.
"""
bounded_side(t::Tetrahedron3, p::Point2)

"""
    has_on_positive_side(t::Tetrahedron3, p::Point2)

Returns `true`, iff `p` lies on the positive side of `t`.
"""
has_on_positive_side(t::Tetrahedron3, p::Point2)

"""
    has_on_negative_side(t::Tetrahedron3, p::Point2)

Returns `true`, iff `p` lies on the negative side of `t`.
"""
has_on_negative_side(t::Tetrahedron3, p::Point2)

"""
    has_on_boundary(t::Tetrahedron3, p::Point2)

Returns `true`, iff `p` lies on the boundary of `t`.
"""
has_on_boundary(t::Tetrahedron3, p::Point2)

"""
    has_on_bounded_side(t::Tetrahedron3, p::Point2)

Returns `true`, iff `p` lies properly inside of `t`.
"""
has_on_bounded_side(t::Tetrahedron3, p::Point2)

"""
    has_on_unbounded_side(t::Tetrahedron3, p::Point2)

Returns `true`, iff `p` lies properly outside of `t`.
"""
has_on_unbounded_side(t::Tetrahedron3, p::Point2)

"""
    volume(t::Tetrahedron3)

Returns the signed volume of `t`.
"""
volume(t::Tetrahedron3)

"""
    bbox(t::Tetrahedron3)

Returns a bounding box containing `t`.
"""
bbox(t::Tetrahedron3)

"""
    transform(t::Tetrahedron3, at::AffTransformation3)

Returns the tetrahedron obtained by applying `at` on the four vertices of `t`.
"""
transform(t::Tetrahedron3, at::AffTransformation3)


text:
using DDR2import

#TODO Network
#TODO Flight Flight Points


text:
using Test, JuLIP, StaticArrays, Printf, LinearAlgebra
using JuLIP.Potentials
using JuLIP.Testing
using JuLIP.Potentials: evaluate_d, evaluate_dd

##

h2("Testing pair potential hessian")

at = bulk(:Cu, cubic=true)
pp = lennardjones(r0=rnn(:Cu))

h3("test the potential (scalar) itself")
rr = 0.9*rnn(:Cu) .+ 3.0*rand(100)
ddJ = [@DD pp(r) for r in rr]
ddJh = [((@D pp(r + 1e-4)) - (@D pp(r - 1e-4))) / 2e-4  for r in rr]
println( @test (@show maximum(abs.(ddJ .- ddJh))) < 1e-4 )

##

h3("test `hess`; with two test vectors")
for R in ( [0.0,-3.61,-3.61], [-1.805,-1.805,-3.61] )
   r = norm(R)
   @show r, R
   E0 = pp(r)
   # f0 = grad(pp, r, JVecF(R))
   f0 = evaluate_d(pp, r) * R/r
   size(f0)
   # df0 = hess(pp, r, JVecF(R))
   df0 = evaluate_dd(pp, r) * (R/r) * (R/r)' +
         (evaluate_d(pp, r)/r) * (I - (R/r) * (R/r)')
   dfh = zeros(3,3)
   fh = zeros(3)
   for p = 2:11
     h = 0.1^p
     for n = 1:3
       R[n] += h
      #  fh[n] = (evaluate(pp, norm(R)) - E0)
       dfh[:, n] = (evaluate_d(pp, norm(R)) * R/norm(R) - f0) / h
         # (grad(pp, r, JVecF(R)) - f0) / h
       R[n] -= h
     end
     @printf("%1.1e | %4.2e \n", h, norm(dfh - df0, Inf))
   end
   @printf("-------|--------- \n")
end

##

h3("full finite-difference test for pairpot `hessian`")
at = at * 2
set_pbc!(at, false)
set_calculator!(at, pp)
println(@test fdtest_hessian( x->JuLIP.gradient(at, x), x->hessian(at, x), dofs(at) ))

##

h2("Testing EAM hessian")
# setup a geometry
at = bulk(:Fe, cubic=true) * 2
set_pbc!(at, false)
eam = eam_Fe
set_calculator!(at, eam)
rattle!(at, 0.1)

h3("test a single stencil")
r = []
R = []
for (idx, _j, R1) in sites(at, cutoff(eam))
   if idx == 3
      global r = norm(R1)
      global R = R1
      break
   end
end

# evaluate site gradient and hessian
dVs = evaluate_d(eam, R)
hVs = evaluate_dd(eam, R)
# and convert them to vector form
dV = mat(dVs)[:]
hV = zeros(3*size(hVs,1), 3*size(hVs,2))
for i = 1:size(hVs,1), j = 1:size(hVs,2)
   hV[3*(i-1).+(1:3), 3*(j-1).+(1:3)] = hVs[i,j]
end
matR = mat(R)

errs = []
for p = 2:9
   h = 0.1^p
   hVh = fill(0.0, size(hV))
   for n = 1:length(matR)
      matR[n] += h
      r = norm.(R)
      dVh = mat(evaluate_d(eam, R))[:]
      hVh[:, n] = (dVh - dV) / h
      matR[n] -= h
   end
   push!(errs, norm(hVh - hV, Inf))
   @printf("%1.1e | %4.2e \n", h, errs[end])
end
println(@test /(extrema(errs)...) < 1e-3)

##

h3("full finite-difference test ...")
h3(" ... EAM forces")
println(@test fdtest( x -> energy(at, x), x -> JuLIP.gradient(at, x), dofs(at) ))
h3(" ... EAM hessian")
println(@test fdtest_hessian( x->gradient(at, x), x->hessian(at, x), dofs(at) ))

##

h2("Testing Stillinger-Weber hessian")

# setup a geometry
at = bulk(:Si, cubic=true) * 2
rattle!(at, 0.02)
set_pbc!(at, false)
sw = StillingerWeber()
set_calculator!(at, sw)

h3("full finite-difference test ...")
h3(" ... SW forces")
println(@test fdtest( x -> energy(at, x), x -> JuLIP.gradient(at, x), dofs(at) ))
h3(" ... SW hessian")
println(@test fdtest_hessian( x->JuLIP.gradient(at, x), x->hessian(at, x), dofs(at) ))


text:
#=
Running order:
1. refresh_norsk_dt_210423()
2. cv_210510()
=#
using DataFrames, JLD2, CSV, GZip, LinearAlgebra, SparseArrays, Statistics, CodecZlib
import RD2:dat_dir
"""
    function read_n_drp(drp)
---
This is to read data sent by Mario on 2021-04-23.
More traits were added.
"""
function read_n_drp(drp)
    @info "  - Norwegian data training set"
    np = DataFrame(ID          = String[],
                   Milk        = Float64[],
                   Milk_ERC    = Float64[],
                   Fat         = Float64[],
                   Fat_ERC     = Float64[],
                   Protein     = Float64[],
                   Protein_ERC = Float64[],
                   SCS         = Float64[],
                   SCS_ERC     = Float64[],
                   x10         = Float64[], # NoInsHeifer DRP
                   x11         = Float64[], # NoInsHeifer ERC
                   x12         = Float64[], # NoInsCow DRP
                   x13         = Float64[], # NoInsCow ERC
                   x14         = Float64[], # calving interval DRP
                   x15         = Float64[], # calving interval ERC
                   N           = Int[])
    for line in eachline(drp)
        v = split(line)
        length(v) != 16 && continue # because two lines have only 2 columns in trnn set
        id = v[1]
        y = parse.(Float64, v[2:end-1])
        n = parse(Int, v[end])
        push!(np, [id; y; n])
    end
    np
end

function refresh_norsk_dt_210423()
    @info "Read phnotypes, EBV and CV sets"
    @load "$dat_dir/jld/drp-training-2021-03-11.jld" dts gts nts
    nts = read_n_drp(joinpath(dat_dir, "mario/20210423/norway_DRPs.txt"))
    @save joinpath(dat_dir, "jld/drp-training-2021-04-23.jld") {compress=true} dts gts nts

    @info "Insert new breeding values of Norwegian data"
    df = CSV.File(joinpath(dat_dir, "data-2020-07-02/norsk_RDC_REL_RBV.csv")) |> DataFrame
    @load "$dat_dir/jld/ebv-all.jld" dbv gbv nbv
    df.ID = string.(df.GenoId)
    df1 = select(df, :ID,
                 :r_Idx_KgProtein305d1_3=>:Protein,
                 :r_Idx_KgFat305d1_3=>:Fat)
    df1[df1.Fat .== ".", :Fat] .= "999"
    df1[df1.Protein .== ".", :Protein] .= "999"
    df1.Fat = parse.(Float64, df1.Fat)
    df1.Protein = parse.(Float64, df1.Protein)
    nbv = innerjoin(nbv, df1, on=:ID)
    @save joinpath(dat_dir, "jld/ebv_202105.jld") {compress=true} dbv gbv nbv
end

function ccwcv(yt, wt, yv, Gi, X)
    nt = length(yt)
    nv = length(yv)
    nx = size(X)[2] + 1

    #Ri = diagm(wt)
    XR = X .* wt
    Z  = sparse([I zeros(nt, nv)])
    ZR = Z .* wt
    lhs = [ XR'X XR'Z
            ZR'X ZR'Z + Gi ] # as va = ve = 1, h2=0.5.
    rhs = [ XR'yt
            ZR'yt ]
    bv = lhs \ rhs
    cc = cor(yv, bv[nt+nx:end])
    return cc, bv[nt+nx:end]
end

function cccv(yt, yv, Gi, X)
    nt = length(yt)
    nv = length(yv)
    nx = size(X)[2] + 1

    Z = sparse([I zeros(nt, nv)])
    lhs = [ X'X X'Z
            Z'X Z'Z + Gi ]      # as va = ve = 1, h2=0.5.
    rhs = [ X'yt
            Z'yt ]
    bv = lhs \ rhs
    cc = cor(yv, bv[nt+nx:end])
    return cc, bv[nt+nx:end]
end

function cv_G(ts, vs, G, trait, erc, X)
    id = ts[:, :ix]
    yt = ts[:, trait]
    wt = ts[:, erc]
    nt = length(id)

    yv = vs[:, trait]
    append!(id, vs.ix)
    nv = length(yv)

    gi = inv(G[id, id])
    c1, _ = ccwcv(yt, wt, yv, gi, X)
    c2, _ = cccv(yt, yv, gi, X)
    return c1, c2
end

function make_X(len...)
    nc = length(len)
    nr = sum(len)
    X = zeros(nr, nc)
    ac = 0
    for i in 1:nc
        X[ac+1:ac+len[i], i] .=1
        ac += len[i]
    end
    X
end

function mk_cv_tbl(dts, dvs, gts, gvs, nts, nvs, G, trait, erc)
    # target is each country in turn
    rst = zeros(6, 7)           # 6 = w/nw by (DGN); 7 = D G N DG DN GN DGN

    @info "Prediction of one country with one of 3 coutries"
    # D → D
    X = make_X(nrow(dts))
    rst[1, 1], rst[4, 1] = cv_G(dts, dvs, G, trait, erc, X)
    # G → D
    X = make_X(nrow(gts))
    rst[1, 2], rst[4, 2] = cv_G(gts, dvs, G, trait, erc, X)
    # N → D
    X = make_X(nrow(nts))
    rst[1, 3], rst[4, 3] = cv_G(nts, dvs, G, trait, erc, X)
    # D → G
    X = make_X(nrow(dts))
    rst[2, 1], rst[5, 1] = cv_G(dts, gvs, G, trait, erc, X)
    # G → G
    X = make_X(nrow(gts))
    rst[2, 2], rst[5, 2] = cv_G(gts, gvs, G, trait, erc, X)
    # N → G
    X = make_X(nrow(nts))
    rst[2, 3], rst[5, 3] = cv_G(nts, gvs, G, trait, erc, X)
    # D → G
    X = make_X(nrow(dts))
    rst[3, 1], rst[6, 1] = cv_G(dts, nvs, G, trait, erc, X)
    # G → G
    X = make_X(nrow(gts))
    rst[3, 2], rst[6, 2] = cv_G(gts, nvs, G, trait, erc, X)
    # N → G
    X = make_X(nrow(nts))
    rst[3, 3], rst[6, 3] = cv_G(nts, nvs, G, trait, erc, X)

    @info "Prediction of one country with two of 3 coutries"
    dgts = vcat(dts, gts)
    dnts = vcat(dts, nts)
    gnts = vcat(gts, nts)
    # D + G → D
    X = make_X(nrow(dts), nrow(gts))
    rst[1, 4], rst[4, 4] = cv_G(dgts, dvs, G, trait, erc, X)
    # D + N → D
    X = make_X(nrow(dts), nrow(nts))
    rst[1, 5], rst[4, 5] = cv_G(dnts, dvs, G, trait, erc, X)
    # G + N → D
    X = make_X(nrow(gts), nrow(nts))
    rst[1, 6], rst[4, 6] = cv_G(gnts, dvs, G, trait, erc, X)
    # D + G → G
    X = make_X(nrow(dts), nrow(gts))
    rst[2, 4], rst[5, 4] = cv_G(dgts, gvs, G, trait, erc, X)
    # D + N → G
    X = make_X(nrow(dts), nrow(nts))
    rst[2, 5], rst[5, 5] = cv_G(dnts, gvs, G, trait, erc, X)
    # G + N → G
    X = make_X(nrow(gts), nrow(nts))
    rst[2, 6], rst[5, 6] = cv_G(gnts, gvs, G, trait, erc, X)
    # D + G → N
    X = make_X(nrow(dts), nrow(gts))
    rst[3, 4], rst[6, 4] = cv_G(dgts, nvs, G, trait, erc, X)
    # D + N → N
    X = make_X(nrow(dts), nrow(nts))
    rst[3, 5], rst[6, 5] = cv_G(dnts, nvs, G, trait, erc, X)
    # G + N → N
    X = make_X(nrow(gts), nrow(nts))
    rst[3, 6], rst[6, 6] = cv_G(gnts, nvs, G, trait, erc, X)

    @info "Prediction of one country with all 3 coutries"
    dgnts = vcat(dts, gts, nts)
    # D + G + N → D
    X = make_X(nrow(dts), nrow(gts), nrow(nts))
    rst[1, 7], rst[4, 7] = cv_G(dgnts, dvs, G, trait, erc, X)
    # D + G + N → G
    X = make_X(nrow(dts), nrow(gts), nrow(nts))
    rst[2, 7], rst[5, 7] = cv_G(dgnts, gvs, G, trait, erc, X)
    # D + G + N → N
    X = make_X(nrow(dts), nrow(gts), nrow(nts))
    rst[3, 7], rst[6, 7] = cv_G(dgnts, nvs, G, trait, erc, X)

    round.(rst, sigdigits = 2)
end

function make_tv(ts, cs, bv, drp, erc)
    t1 = filter(row ->row.st == 't', cs)
    df1 = select(innerjoin(t1, ts, on = :ID), :ix, drp, erc)
    dft = filter(row -> row[drp] < 999, df1)
    t2  = filter(row -> row.st == 'v', cs)
    dfv = select(innerjoin(t2, bv, on = :ID), :ix, drp)
    return dft, dfv
end

function cv_210510()
    @load "$dat_dir/jld/cv-setup.jld" dcs gcs ncs
    @load "$dat_dir/jld/drp-training-2021-04-23.jld" dts gts nts
    @load "$dat_dir/jld/ebv_202105.jld"  dbv gbv nbv
    @load "$dat_dir/jld/GAID.jld"     G A ID

    begin                       # The milk trait
        dt, dv = make_tv(dts, dcs, dbv, :Milk, :Milk_ERC)
        gt, gv = make_tv(gts, gcs, gbv, :Milk, :Milk_ERC)
        nt, nv = make_tv(nts, ncs, nbv, :Milk, :Milk_ERC)
        rst = mk_cv_tbl(dt, dv, gt, gv, nt, nv, G, :Milk, :Milk_ERC)
        for row in eachrow(rst)
            println("| ", join(row, " | "), " |")
        end
    end
    begin                       # The protein trait
        dt, dv = make_tv(dts, dcs, dbv, :Protein, :Protein_ERC)
        gt, gv = make_tv(gts, gcs, gbv, :Protein, :Protein_ERC)
        nt, nv = make_tv(nts, ncs, nbv, :Protein, :Protein_ERC)
        rst = mk_cv_tbl(dt, dv, gt, gv, nt, nv, G, :Protein, :Protein_ERC)
        for row in eachrow(rst)
            println("| ", join(row, " | "), " |")
        end
    end
    begin                       # The fat trait
        dt, dv = make_tv(dts, dcs, dbv, :Fat, :Fat_ERC)
        gt, gv = make_tv(gts, gcs, gbv, :Fat, :Fat_ERC)
        nt, nv = make_tv(nts, ncs, nbv, :Fat, :Fat_ERC)
        rst = mk_cv_tbl(dt, dv, gt, gv, nt, nv, G, :Fat, :Fat_ERC)
        for row in eachrow(rst)
            println("| ", join(row, " | "), " |")
        end
    end
    begin                       # The SCS trait
        dt, dv = make_tv(dts, dcs, dbv, :SCS, :SCS_ERC)
        gt, gv = make_tv(gts, gcs, gbv, :SCS, :SCS_ERC)
        nt, nv = make_tv(nts, ncs, nbv, :SCS, :SCS_ERC)
        rst = mk_cv_tbl(dt, dv, gt, gv, nt, nv, G, :SCS, :SCS_ERC)
        for row in eachrow(rst)
            println("| ", join(row, " | "), " |")
        end
    end
end


text:
# functions related to student's T distribution

import .RFunctions:
    tdistpdf,
    tdistlogpdf,
    tdistcdf,
    tdistccdf,
    tdistlogcdf,
    tdistlogccdf,
    tdistinvcdf,
    tdistinvccdf,
    tdistinvlogcdf,
    tdistinvlogccdf


text:
module SimpleRegressionGP

    using Optim, Random, Distributions, Printf, LinearAlgebra

    using PyPlot, MiscUtil

    include("kernels.jl")
    include("simpleregressiongp.jl")

    export simpleregressiongp
end


text:
using DifferentialEquations
using StochasticDiffEq

include("../src/utils.jl")

function generateNoisyTetradData(u0, tsteps)
    numTrajs = size(u0)[2];
    tspan = (tsteps[1], tsteps[end]);
    c = reshape(createNoiseRateForcingHIT(15,15), (9,9));

    function tetradDrift(u, p, t)
        numDims, numTrajs = size(u);
        rhoPiece = zeros(9,numTrajs);
        vgtPiece = zeros(9,numTrajs);
        for i in 1:numTrajs
            ρ = reshape(u[1:9],(3,3));
            M = reshape(u[10:18],(3,3));
            k = inv(ρ);
            Π = (k*k')/tr(k*k');
            rhoPiece[:,i] = reshape(M'*ρ, 9);
            vgtPiece[:,i] = reshape(-M^2 + tr(M^2)*Π, 9);
        end
        return cat(rhoPiece, vgtPiece, dims=1);
    end

    function hitDiff(u,p,t)
        return [Matrix(I,9,9) zeros(9,9)
                zeros(9,9)    c];
    end
    W = WienerProcess(0.0, u0);

    sdeProb = SDEProblem(tetradDrift,
                         hitDiff,
                         u0[:,1],
                         tspan,
                         noise=W,
                         noise_rate_prototype=zeros(18,18));

    sol = solve(sdeProb,
                u0 = u0,
                EM(),
                tstops = tsteps,
                saveat = tsteps);

    return permutedims(Array(sol), [1,3,2]);
end

max_stars_repo_name: cmhyett/MLCoarseGrainedVGT

max_issues_repo_name: cmhyett/MLCoarseGrainedVGT

max_forks_repo_name: cmhyett/MLCoarseGrainedVGT


text:
# using Pkg; for p in ("Knet","ArgParse"); haskey(Pkg.installed(),p) || Pkg.add(p); end

"""

julia treenn.jl # to use with default options
julia treenn.jl -h # to see all options with default values

This example implements a binary tree-structured LSTM networks proposed
in 'Improved Semantic Representations From Tree-Structured Long Short-Term
Memory Networks', Kai Sheng Tai, Richard Socher, Christopher D. Manning,
arXiv technical report 1503.00075, 2015.

* Paper url: https://arxiv.org/pdf/1503.00075.pdf
* Project page: https://github.com/stanfordnlp/treelstm

"""
module TreeLSTM
using Knet, CUDA, ArgParse, Dates, Printf, Random

const UNK = "_UNK_"
t00 = now()
include(Knet.dir("data","treebank.jl"))

function main(args)
    s = ArgParseSettings()
    s.description = "Tree-structured LSTM network in Knet."

    @add_arg_table s begin
        # ("--usegpu"; action=:store_true; help="use GPU or not")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"); help="array type: Array for cpu, KnetArray for gpu")
        ("--embed"; arg_type=Int; default=128; help="word embedding size")
        ("--hidden"; arg_type=Int; default=128; help="LSTM hidden size")
        ("--timeout"; arg_type=Int; default=600; help="max timeout (in seconds)")
        ("--epochs"; arg_type=Int; default=3; help="number of training epochs")
        ("--minoccur"; arg_type=Int; default=0; help="word min occurence limit")
        ("--report"; arg_type=Int; default=1000; help="report period (in iters)")
        ("--seed"; arg_type=Int; default=-1; help="random seed")
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    println(o)
    o[:seed] > 0 && Knet.seed!(o[:seed])
    # atype = o[:atype] = !o[:usegpu] ? Array{Float32} : KnetArray{Float32}
    atype = eval(Meta.parse(o[:atype]))
    datadir = abspath(joinpath(@__DIR__, "../data/trees"))
    datadir = isdir(datadir) ? datadir : TREEBANK_DIR

    # read data
    trn, dev = load_treebank_data(datadir)

    # build vocabs
    l2i, w2i, i2l, i2w = build_treebank_vocabs(trn)
    nwords = length(w2i); nlabels = length(l2i)

    make_data!(trn, w2i, l2i)
    make_data!(dev, w2i, l2i)

    # build model
    w = initweights(atype, o[:hidden], nwords, nlabels, o[:embed])
    opt = optimizers(w, Adam)

    # main loop
    println("startup time: ", Int((now()-t00).value)*0.001); flush(stdout)
    all_time = sents = 0
    o[:timeout] = o[:timeout] <= 0 ? Inf : o[:timeout]
    for epoch = 1:o[:epochs]
        closs = cwords = 0
        shuffle!(trn)
        t0 = now()
        for k = progress(1:length(trn))
            sents += 1
            iter = (epoch-1)*length(trn) + k
            tree = trn[k]
            this_loss, this_words = train!(w,tree,opt)
            closs += this_loss
            cwords += this_words

            if o[:report] > 0 && iter % o[:report] == 0
                @printf("%f\n", closs/cwords); flush(stdout)
                closs = cwords = 0
            end
        end
        all_time += Int((now()-t0).value)*0.001

        good = bad = 0
        for tree in dev
            ind, nwords = predict(w, tree)
            ypred = i2l[ind]
            ygold = tree.label
            if ypred == ygold
                good += 1
            else
                bad += 1
            end
        end
        @printf(
            "acc=%.4f, time=%.4f, sent_per_sec=%.4f\n",
            good/(good+bad), all_time, sents/all_time); flush(stdout)

        all_time > o[:timeout] && return
    end
end

function initweights(atype, hidden, words, labels, embed, winit=0.01)
    w = Array{Any}(undef,9)
    w[1] = winit*randn(3*hidden, embed)
    w[2] = zeros(3*hidden, 1)
    w[3] = winit*randn(3*hidden, 2*hidden)
    w[4] = zeros(3*hidden, 1)
    w[5] = winit*randn(hidden, hidden)
    w[6] = winit*randn(hidden, hidden)
    w[7] = ones(hidden,1)
    w[8] = winit*randn(labels, hidden)
    w[9] = winit*randn(embed, words)
    return map(i->convert(atype, i), w)
end

function lstm(w,ind)
    x = w[end][:,ind]
    x = reshape(x, length(x), 1)
    hsize = size(x,1)
    gates = w[1] * x .+ w[2]
    i = sigm.(gates[1:hsize,:])
    o = sigm.(gates[1+hsize:2hsize,:])
    u = sigm.(gates[1+2hsize:3hsize,:])
    c = i .* u
    h = o .* tanh.(c)
    return (h,c)
end

function slstm(w,h1,h2,c1,c2)
    hsize = size(h1,1)
    h = vcat(h1,h2)
    gates = w[3] * h .+ w[4]
    i  = sigm.(gates[1:hsize,:])
    o  = sigm.(gates[1+hsize:2hsize,:])
    u  = sigm.(gates[1+2hsize:3hsize,:])
    f1 = sigm.(w[5] * h1 .+ w[7])
    f2 = sigm.(w[6] * h2 .+ w[7])
    c  = i .* u .+ f1 .* c1 .+ f2 .* c2
    h  = o .* tanh.(c)
    return (h,c)
end

let
    global traverse
    function traverse(w, tree)
        h,c,hs,ys = helper(w,tree, Any[], Any[])
        return hs,ys
    end

    function helper(w,t,hs,ys)
        h = c = nothing
        if length(t.children) == 1 && isleaf(t.children[1])
            l = t.children[1]
            h,c = lstm(w,l.data)
        elseif length(t.children) == 2
            t1,t2 = t.children[1], t.children[2]
            h1,c1,hs,ys = helper(w,t1,hs,ys)
            h2,c2,hs,ys = helper(w,t2,hs,ys)
            h,c = slstm(w,h1,h2,c1,c2)
        else
            error("invalid tree")
        end
        return (h,c,[hs...,h],[ys...,t.data])
    end
end

# treenn loss function
function loss(w, tree, values=[])
    hs, ygold = traverse(w, tree)
    ypred = w[end-1] * hcat(hs...)
    len = length(ygold)
    (lossval, cnt) = nll(ypred,ygold; average=false)
    @assert len == cnt
    push!(values, lossval); push!(values, len)
    return lossval/len
end

# tag given input sentence
function predict(w,tree)
    total = 0
    hs, ys = traverse(w, tree)
    ypred = w[end-1] * hs[end]
    ypred = convert(Array{Float32}, ypred)[:]
    return (argmax(ypred),length(ys))
end

lossgradient = grad(loss)

function train!(w,tree,opt)
    values = []
    gloss = lossgradient(w, tree, values)
    update!(w,gloss,opt)
    return (values...,)
end

splitdir(PROGRAM_FILE)[end] == "treenn.jl" && main(ARGS)

end # module


text:
#=
P31 (**) Determine whether a given integer number is prime.
Example:
* (is-prime 7)
T
=#

function is_prime(n)
    # Using built-in `isprime` function :)
    isprime(n)
end

@assert is_prime(2)
@assert is_prime(3)
@assert is_prime(5)
@assert is_prime(7)
@assert is_prime(11)
@assert is_prime(999983)  # Biggest prime number under 1,000,000

println("Tests passed: JL-31.jl")


text:
@testset "send_transmit.jl" begin
    @testset "Three short writes" begin
        isn = rand(UInt32(0):typemax(UInt32))

        sender = TCPSender(; fixed_isn=WrappingInt32(isn))
        fill_window!(sender)

        expect_seg(sender; no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_state(sender, JLSponge.SYN_ACKED)
        write_bytes!(sender, "ab")
        expect_seg(sender; data="ab", seqno=WrappingInt32(isn + 1))
        write_bytes!(sender, "cd")
        expect_seg(sender; data="cd", seqno=WrappingInt32(isn + 3))
        write_bytes!(sender, "abcd")
        expect_seg(sender; data="abcd", seqno=WrappingInt32(isn + 5))
        @test next_seqno(sender) == WrappingInt32(isn + 9)
        @test sender.bytes_in_flight == 8
        expect_state(sender, JLSponge.SYN_ACKED)
    end

    @testset "Many short writes, continuous acks" begin
        isn = rand(UInt32(0):typemax(UInt32))

        sender = TCPSender(; fixed_isn=WrappingInt32(isn))
        fill_window!(sender)

        expect_seg(sender; no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_state(sender, JLSponge.SYN_ACKED)
        max_block_size = 10
        n_rounds = 10000
        bytes_sent = 0

        for i in 1:n_rounds
            block_size = rand(1:max_block_size)
            data = join(rand('a':'z', block_size))
            @test next_seqno(sender) == WrappingInt32(isn + bytes_sent + 1)
            write_bytes!(sender, data)
            bytes_sent += block_size
            @test sender.bytes_in_flight == block_size
            expect_seg(sender; seqno=WrappingInt32(isn + 1 + bytes_sent - block_size),
                       data=data)
            expect_no_seg(sender)
            ack_received_test(sender, WrappingInt32(isn + 1 + bytes_sent))
        end
    end
    @testset "Many short writes, ack at end" begin
        isn = rand(UInt32(0):typemax(UInt32))

        sender = TCPSender(; fixed_isn=WrappingInt32(isn))
        fill_window!(sender)

        expect_seg(sender; no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 65000)
        expect_state(sender, JLSponge.SYN_ACKED)
        max_block_size = 10
        n_rounds = 10000
        bytes_sent = 0

        for i in 1:n_rounds
            block_size = rand(1:max_block_size)
            data = join(rand('a':'z', block_size))
            @test next_seqno(sender) == WrappingInt32(isn + bytes_sent + 1)
            write_bytes!(sender, data)
            bytes_sent += block_size
            @test sender.bytes_in_flight == bytes_sent
            expect_seg(sender; seqno=WrappingInt32(isn + 1 + bytes_sent - block_size),
                       data=data)
            expect_no_seg(sender)
        end
        @test sender.bytes_in_flight == bytes_sent
        ack_received_test(sender, WrappingInt32(isn + 1 + bytes_sent))
        @test sender.bytes_in_flight == 0
    end

    @testset "Window filling" begin
        isn = rand(UInt32(0):typemax(UInt32))

        sender = TCPSender(; fixed_isn=WrappingInt32(isn))
        fill_window!(sender)

        expect_seg(sender; no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 3)
        expect_state(sender, JLSponge.SYN_ACKED)
        write_bytes!(sender, "01234567")
        @test sender.bytes_in_flight == 3
        expect_seg(sender; data="012")
        expect_no_seg(sender)
        @test next_seqno(sender) == WrappingInt32(isn + 1 + 3)
        ack_received_test(sender, WrappingInt32(isn + 1 + 3), 3)
        @test sender.bytes_in_flight == 3
        expect_seg(sender; data="345")
        expect_no_seg(sender)
        @test next_seqno(sender) == WrappingInt32(isn + 1 + 6)
        ack_received_test(sender, WrappingInt32(isn + 1 + 6), 3)
        @test sender.bytes_in_flight == 2
        expect_seg(sender; data="67")
        expect_no_seg(sender)
        @test next_seqno(sender) == WrappingInt32(isn + 1 + 8)
        ack_received_test(sender, WrappingInt32(isn + 1 + 8), 3)
        @test sender.bytes_in_flight == 0
        expect_no_seg(sender)
    end

    @testset "Window filling" begin
        isn = rand(UInt32(0):typemax(UInt32))

        sender = TCPSender(; fixed_isn=WrappingInt32(isn))
        fill_window!(sender)

        expect_seg(sender; no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        ack_received_test(sender, WrappingInt32(isn + 1), 3)
        expect_state(sender, JLSponge.SYN_ACKED)
        write_bytes!(sender, "01")
        @test sender.bytes_in_flight == 2
        expect_seg(sender; data="01")
        expect_no_seg(sender)
        @test next_seqno(sender) == WrappingInt32(isn + 1 + 2)
        write_bytes!(sender, "23")
        @test sender.bytes_in_flight == 3
        expect_seg(sender; data="2")
        expect_no_seg(sender)
        @test next_seqno(sender) == WrappingInt32(isn + 1 + 3)
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
"""
	options = NewtonPar(tol = 1e-4,...)

Returns a variable containing parameters to affect the `newton` algorithm when solving `F(x) = 0`.

# Arguments (with default values):
- `tol = 1e-10`: absolute tolerance for `F(x)`
- `maxIter = 50`: number of Newton iterations
- `verbose = false`: display Newton iterations?
- `linsolver = DefaultLS()`: linear solver, must be `<: AbstractLinearSolver`
- `eigsolver = DefaultEig()`: eigen solver, must be `<: AbstractEigenSolver`

# Arguments only used in `newtonPALC`
- `linesearch = false`: use line search algorithm
- `alpha = 1.0`: alpha (damping) parameter for line search algorithm
- `almin  = 0.001 `: minimal vslue of the damping `alpha`

!!! tip "Mutating"
    For performance reasons, we decided to use an immutable structure to hold the parameters. One can use the package `Setfield.jl` to drastically simplify the mutation of different fields. See tutorials for more examples.
"""
@with_kw struct NewtonPar{T, L <: AbstractLinearSolver, E <: AbstractEigenSolver}
	tol::T			 = 1e-10
	maxIter::Int64 	 = 50
	alpha::T         = 1.0        # damping
	almin::T         = 0.001      # minimal damping
	verbose::Bool    = false
	linesearch::Bool = false
	linsolver::L 	 = DefaultLS()
	eigsolver::E 	 = DefaultEig()
end

####################################################################################################
"""
		newton(F, J, x0, options::NewtonPar; normN = norm, callback = (x, f, J, res, iteration, optionsN; kwargs...) -> true, kwargs...)

This is the Newton Solver for `F(x) = 0` with Jacobian `J` and initial guess `x0`. The function `normN` allows to specify a norm for the convergence criteria. It is important to set the linear solver `options.linsolver` properly depending on your problem. This linear solver is used to solve ``J(x)u = -F(x)`` in the Newton step. You can for example use `linsolver = Default()` which is the operator backslash: it works well for Sparse / Dense matrices. See [Linear solvers](@ref) for more informations.

# Arguments:
- `x -> F(x)` functional whose zeros are looked for. In particular, it is not **inplace**,
- `dF(x) = x -> J(x)` compute the jacobian of `F` at `x`. It is then passed to `options.linsolver`. The Jacobian `J(x)` can be a matrix or an out-of-place function.
- `x0` initial guess
- `options` variable holding the internal parameters used by the `newton` method
- `callback` function passed by the user which is called at the end of each iteration. Can be used to update a preconditionner for example. The `optionsN` will be `options` passed in order to change the linear / eigen solvers
- `kwargs` arguments passed to the callback. Useful when `newton` is called from `continuation`

Simplified calls are provided, for example when `J` is not passed. It is then computed with finite differences.

# Output:
- solution:
- history of residuals
- flag of convergence
- number of iterations
"""
function newton(Fhandle, Jhandle, x0, options::NewtonPar{T}; normN = norm, callback = (x, f, J, res, iteration, optionsN; kwargs...) -> true, kwargs...) where T
	# Extract parameters
	@unpack tol, maxIter, verbose, linesearch = options

	# Initialise iterations
	x = similar(x0); copyto!(x, x0) # x = copy(x0)
	f = Fhandle(x)
	d = similar(f); copyto!(d, f)	# d = copy(f)

	neval = 1
	res = normN(f)
	resHist = [res]
	it = 0

	# Displaying results
	verbose && displayIteration(it, neval, res)

	# Main loop
	while (res > tol) & (it < maxIter)
		J = Jhandle(x)
		d, _, itlinear = options.linsolver(J, f)

		# Update solution: x .= x .- d
		minus!(x, d)

		copyto!(f, Fhandle(x))
		res = normN(f)

		neval += 1
		push!(resHist, res)
		it += 1

		callback(x, f, J, res, it, options; kwargs...) == false && (it = maxIter)
		verbose && displayIteration(it, neval, res, itlinear)
	end
	(resHist[end] > tol) && @error("\n--> Newton algorithm failed to converge, residual = $(res[end])")
	return x, resHist, resHist[end] < tol, it
end

# simplified call to newton when no Jacobian is passed in which case we estimate it using finiteDifferences
function newton(Fhandle, x0, options:: NewtonPar{T};kwargs...) where T
	Jhandle = u -> finiteDifferences(Fhandle, u)
	return newton(Fhandle, Jhandle, x0, options; kwargs...)
end


text:
for N in [Float64, Rational{Int}, Float32]
    # random singleton
    rand(Singleton)

    # center/radius/high/low/generators
    c = N[0, 0]
    r = N[0, 0]
    s = Singleton(c)
    @test center(s) == c
    @test radius_hyperrectangle(s) == r
    @test high(s) == r
    @test low(s) == -r
    for i in 1:2
        @test radius_hyperrectangle(s, i) == r[i]
        @test high(s, i) == r[i]
        @test low(s, i) == -r[i]
    end
    @test collect(generators(s)) == Vector{Vector{N}}()
    @test genmat(s) == Matrix{N}(undef, 2, 0)
    @test ngens(s) == 0

    # 1D singleton
    s = Singleton(N[1])
    d = N[1]
    @test σ(d, s) == N[1]
    d = N[-1]
    @test σ(d, s) == N[1]
    d = N[0]
    @test σ(d, s) == N[1]

    # 2D singleton
    s = Singleton(N[1, 2])
    d = N[1, 0]
    @test σ(d, s) == N[1, 2]
    d = N[-1, 0.5]
    @test σ(d, s) == N[1, 2]
    d = N[0, 0]
    @test σ(d, s) == N[1, 2]

    # test that abstract vectors can be used as well
    S = Singleton(SingleEntryVector(1, 1000, 0.5))
    @test element(S) == SingleEntryVector(1, 1000, 0.5)

    # boundedness
    @test isbounded(s)

    # element function
    @test element(s) == s.element
    for i in 1:2
        @test element(s, i) == s.element[i]
    end

    # isempty
    @test !isempty(s)

    # isuniversal
    answer, w = isuniversal(s, true)
    @test !isuniversal(s) && !answer && w ∉ s

    # membership
    S = Singleton(N[1, 1])
    @test N[0.9, 1.1] ∉ S
    @test N[1, 1] ∈ S

    # check membership error message
    @test_throws ErrorException S.element ⊆ S

    # an_element function
    @test an_element(S) ∈ S

    # vertices_list
    @test vertices_list(S)[1] ∈ S

    # radius_hyperrectangle
    @test iszero(radius_hyperrectangle(S))

    # high and low
    @test high(S) == N[1, 1]
    @test low(S) == N[1, 1]

    # concrete linear map
    M = N[0 1; -1 0]
    @test element(linear_map(M, S)) == an_element(M * S)

    # translation
    @test translate(s, N[1, 2]) == Singleton(N[2, 4])

    # subset
    s1 = Singleton(N[0, 1])
    s2 = Singleton(N[0, 3])
    p1 = VPolygon([N[0,0], N[0, 2]])
    p2 = VPolygon([N[0,0], N[0, 2], N[2, 0]])
    b = BallInf(N[0, 1], N(1))
    @test s1 ⊆ p1 && ⊆(s1, p1, true)[1]
    subset, point = ⊆(s2, p2, true)
    @test s2 ⊈ p2 && !subset && point ∈ s2 && point ∉ p2
    @test s1 ⊆ s1 && ⊆(s1, s1, true)[1]
    subset, point = ⊆(s1, s2, true)
    @test s1 ⊈ s2 && !subset && point ∈ s1 && point ∉ s2
    subset, point = ⊆(s1, b, true)
    @test subset && s1 ⊆ b
    subset, point = ⊆(s2, b, true)
    @test s2 ⊈ b && !subset && point ∈ s2 && point ∉ b

    # intersection
    S1 = Singleton(N[1, 1])
    S2 = Singleton(N[0, 0])
    S3 = ZeroSet{N}(2)
    H = BallInf(N[1, 1], N(0.5))
    M = LinearMap(N[1 0; 0 1], H)
    @test intersection(S1, H) == intersection(H, S1) == S1
    @test intersection(S2, H) == intersection(H, S2) == EmptySet{N}(2)
    @test intersection(S3, H) == intersection(H, S3) == EmptySet{N}(2)
    @test is_intersection_empty(S1, S2) && is_intersection_empty(S1, S2, true)[1]
    intersection_empty, point = is_intersection_empty(S2, S3, true)
    @test !is_intersection_empty(S2, S3) && !intersection_empty &&
        point ∈ S2 && point ∈ S3
    for X in [H, M]
        intersection_empty, point = is_intersection_empty(S1, X, true)
        @test !is_intersection_empty(S1, X) && !intersection_empty &&
            point ∈ S1 && point ∈ X
        @test is_intersection_empty(S2, X) &&
            is_intersection_empty(S2, X, true)[1]
    end
end


text:
@testset "Sentiment" begin
    m = SentimentAnalyzer()

    d = StringDocument("a very nice thing that everyone likes")

    @test m(d) > 0.5

    d = StringDocument("a horrible thing that everyone hates")

    @test m(d) < 0.5

    # testing default behaviour of handle_unknown
    d = StringDocument("some sense and some nonSense")

    @test m(d) < 0.5

    # testing behaviour of words which are present in dictionary but do not have embedding assigned
    d = StringDocument("some sense and some duh")

    @test m(d) < 0.5

    # testing user given handle_unknown function
    d = Document("a Horrible thing that Everyone Hates")

    @test m(d, (x) -> [lowercase(x)]) < 0.5

    # Make it throw an error when unknown word encountered
    d = Document("some sense and some Hectic")

    @test_throws ErrorException m(d, (x) -> error("OOV word $x encountered"))

end


text:
include("leapfrog.jl")

include("simple_gauss_hmc.jl")
include("gdemo_hmc.jl")
include("gdemo_nuts.jl")


text:
function parse_map(s)
    lines = readlines(IOBuffer(s))
    is_seat = falses((length(lines) + 2, length(lines[1]) + 2))
    for i in 1:length(lines)
        for j in 1:length(lines[1])
            is_seat[i + 1, j + 1] = lines[i][j] == 'L'
        end
    end
    return is_seat
end

function count_neighbors(is_occupied::BitArray{2}, i::Int, j::Int)::Int
    count::Int = 0
    for i_offset in (-1, 0, 1)
        for j_offset in (-1, 0, 1)
            if i_offset == 0 && j_offset == 0
            elseif is_occupied[i + i_offset, j + j_offset]
                count += 1
            end
        end
    end
    return count
end

function run(s)
    # Your code here
    is_seat = parse_map(s)
    is_occupied = copy(is_seat)
    old_is_occupied = copy(is_seat)
    n_rows, n_cols = size(is_seat)

    has_changed::Bool = true
    new_free::Bool = false
    new_occupied::Bool = false
    n_neighbors::Int = 0
    while has_changed
        # is_seat .& .!is_occupied = not_occupied
        has_changed = false
        old_is_occupied .= is_occupied
        for i in 2:n_rows - 1
            for j in 2:n_cols - 1
                if is_seat[i, j]
                    n_neighbors = count_neighbors(old_is_occupied, i, j)
                    new_occupied = is_seat[i, j] & !old_is_occupied[i, j] & (n_neighbors == 0)
                    new_free = old_is_occupied[i, j] & (n_neighbors >= 4)
                    if new_occupied
                        has_changed = true
                        is_occupied[i, j] = true
                    elseif new_free
                        has_changed = true
                        is_occupied[i, j] = false
                    end
                end
            end
        end
    end
    return sum(is_occupied)
end

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

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

main()


text:
@testset "Modulo Test" begin
    @test ALifeBenchmark.mod_1(1, 10) == 1
    @test ALifeBenchmark.mod_1(0, 10) == 10
    @test ALifeBenchmark.mod_1(10, 10) == 10
    @test ALifeBenchmark.mod_1(15, 10) == 5
    @test ALifeBenchmark.mod_1(-2, 10) == 8
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using NbodyGradient, DelimitedFiles

#="""Setup NBG and output initial conditions for rebound and TTVFast"""
function setup_sim(ntransits, denom)
    elements = readdlm("elements_noprior_students.txt",',')
    elements[:,end] .= 0; # Set Ω to 0
    ic = ElementsIC(0.0, 8, elements)

    h = elements[2,2]/denom
    tmax = ntransits*elements[end,2]
    intr = Integrator(h, 0.0, tmax)

    s = State(ic)
    tt = TransitTiming(intr.tmax, ic)
    output_rebound_ICs(s)
    output_TTVFast_ICs(s)
    return s, tt, intr
end
=#
function setup_sim(tmax, denom)
    elements = readdlm("elements_noprior_students.txt",',')
    elements[:,end] .= 0; # Set Ω to 0
    ic = ElementsIC(0.0, 8, elements)

    h = elements[2,2]/denom
    intr = Integrator(h, 0.0, tmax)

    s = State(ic)
    tt = TransitTiming(tmax, ic)
    output_rebound_ICs(s)
    output_TTVFast_ICs(s)
    return s, tt, intr
end

function output_NBG_times(s,tt,intr)
    intr(s,tt,grad=false)
    times = tt.tt[2:end,1:3000]
    open("nbg_tts.txt", "w") do io
        writedlm(io, eachrow(times),',')
    end
    return times
end

function output_rebound_ICs(s::State)
    # Linear transformation to rebound coordinate system.
    A = [0 0 -1; 1 0 0; 0 -1 0]
    pos = copy(s.x)
    vel = copy(s.v)
    reb_pos = zeros(size(pos))
    reb_vel = zeros(size(vel))
    for i in 1:8
        reb_pos[:,i] .= A*pos[:,i]
        reb_vel[:,i] .= A*vel[:,i]
    end
    open("rebound_pos.txt", "w") do io
        writedlm(io, eachrow(reb_pos),',')
    end
    open("rebound_vel.txt", "w") do io
        writedlm(io, eachrow(reb_vel),',')
    end
end

function output_TTVFast_ICs(s::State)
    # Transform to ttvfast astrocentric coordinates
    xstar = copy(s.x[:,1])
    vstar = copy(s.v[:,1])
    ttvf_pos = copy(s.x[:,2:end])
    ttvf_vel = copy(s.v[:,2:end])
    for i in 1:7
        ttvf_pos[:,i] .-= xstar
        ttvf_vel[:,i] .-= vstar
    end
    A = [-1 0 0; 0 1 0; 0 0 -1]
    ttvf_pos = (A * ttvf_pos)'
    ttvf_vel = (A * ttvf_vel)'

    open("ttvfast_ics.txt", "w") do io
        writedlm(io, NbodyGradient.GNEWT) # set units
        writedlm(io, s.m[1])
        writedlm(io, eachrow([s.m[2:end] ttvf_pos ttvf_vel]), ' ')
    end
    run(`mv ttvfast_ics.txt ttvfast/`)
end

max_forks_repo_name: langfzac/nbg-papers


text:
####################### core taylor expansion based algorithm ##################
# still based on BarnesHutFactorization, should probably introduce new TaylorFactorization type
# use_com = true Taylor expands around the center of mass,
# improves on accuracy of BarnesHut if weight vectors have negative elements
# Generally, taylor! is ~2x faster and slightly more accurate than splitting_barneshut!
# when a signficiant degree of compression is taking place (i.e larger θ)
function taylor!(b::AbstractVector, F::BarnesHutFactorization{<:Number}, w::AbstractVector,
                    α::Number = 1, β::Number = 0, θ::Real = F.θ; use_com::Bool = true)
    size(F, 2) == length(w) || throw(DimensionMismatch("length of w does not match second dimension of F: $(length(w)) ≠ $(size(F, 2))"))
    eltype(b) == promote_type(eltype(F), eltype(w)) ||
            throw(TypeError("eltype of target vector b not equal to eltype of matrix-vector product: $(eltype(b)) and $(promote_type(eltype(F), eltype(w)))"))

    @. b = iszero(β) ? 0 : β * b # this avoids trouble if b is initialized with NaN's, e.g. thorugh "similar"
    f_orders(r²) = value_derivative(F.k, r²)
    sums_w = node_sums(w, F.Tree) # IDEA: could pre-allocate
    sums_w_r = weighted_node_sums(w, F.y, F.Tree)
    centers = use_com ? compute_centers_of_mass(F, w) : get_hyper_centers(F)
    @. sums_w_r -= sums_w * centers # need to center the moments
    @threads for i in eachindex(F.x) # exactly 4 * length(y) allocations?
        Fw_i = taylor_recursion(1, F.k, f_orders, F.x[i], F.y, w, sums_w, sums_w_r, θ, F.Tree, centers)::eltype(b) # x[i] creates an allocation
        b[i] += α * Fw_i
    end
    if !isnothing(F.D) # if there is a diagonal correction, need to add it
        mul!(b, F.D, w, α, 1)
    end

    return b
end

# uses taylor expansion as approximation (identical with bh_recursion if center is center of mass)
function taylor_recursion(index, k, f_orders, xi, y::AbstractVector, w::AbstractVector,
                      sums_w::AbstractVector, sums_w_r::AbstractVector, θ::Real, T::BallTree, centers)
    h = T.hyper_spheres[index]
    if isleaf(T.tree_data.n_internal_nodes, index) # do direct computation
        elty = promote_type(eltype(xi), eltype(eltype(y)), eltype(w))
        val = zero(elty)
        @inbounds @simd for i in get_leaf_range(T.tree_data, index) # loop vectorization?
            j = T.indices[i]
            val += k(xi, y[j]) * w[j]
        end
        return val

    elseif h.r < θ * euclidean(xi, centers[index]) # compress
        ri = difference(xi, centers[index])
        sum_abs2_ri = sum(abs2, ri)
        f0, f1 = f_orders(sum_abs2_ri) # contains first and second order derivative evaluation
        value = f0 * sums_w[index]
        first_order = -2*f1 * dot(ri, sums_w_r[index])
        # second_order ...
        return value + first_order

    else # recurse NOTE: parallelizing here is not as efficient as parallelizing over target points
        l = taylor_recursion(getleft(index), k, f_orders, xi, y, w, sums_w, sums_w_r, θ, T, centers)
        r = taylor_recursion(getright(index), k, f_orders, xi, y, w, sums_w, sums_w_r, θ, T, centers)
        return l + r
    end
end

# this represents the Array with values r[i_1]*r[i_2]*...*r[i_p] lazily
# it will be useful for fast approximate matrix multiplication based on
# higher order taylor expansions
struct PowersArray{T, X <: AbstractVector{T}, P} <: AbstractArray{T, P}
    x::X
end

function PowersArray(x, p::Int)
    p ≥ 1 || throw(DomainError("Can't construct PowersArray with p < 1"))
    PowersArray{eltype(x), typeof(x), p}(x)
end
power(A::PowersArray{<:Any, <:Any, P}) where P = P

Base.length(A::PowersArray) = length(A.x)^power(A)
Base.size(A::PowersArray) = tuple((size(A, 1) for _ in 1:power(A))...)
function Base.size(A::PowersArray, i::Int)
    1 ≤ i ≤ power(A) ? length(A.x) : 1
end

function Base.getindex(A::PowersArray{<:Any, <:Any, P}, indices::Vararg{Int, P}) where P
    r = one(eltype(A.x))
    for i in indices
        r *= A.x[i]
    end
    return r
end


text:
using ModelingToolkit
@variables x y
@parameters a b
loss = (a - x)^2 + b * (y - x^2)^2
sys1 = OptimizationSystem(loss,[x,y],[a,b],name=:sys1)
sys2 = OptimizationSystem(loss,[x,y],[a,b],name=:sys2)

@variables z
@parameters β
loss2 = sys1.x - sys2.y + z*β
combinedsys = OptimizationSystem(loss2,[z],[β],systems=[sys1,sys2],name=:combinedsys)

equations(combinedsys)
states(combinedsys)
parameters(combinedsys)

calculate_gradient(combinedsys)
calculate_hessian(combinedsys)
generate_function(combinedsys)
generate_gradient(combinedsys)
generate_hessian(combinedsys)


text:
module PandorasBoxMarioClearPipe

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/clearPipe" ClearPipe(x::Integer, y::Integer, texture::String="glass", surfaceSound::Integer=-1, hasPipeSolids::Bool=true, transportSpeed::Number=175.0, pipeWidth::Int=32)

const placements = Ahorn.PlacementDict(
    "Clear Pipe (Pandora's Box)" => Ahorn.EntityPlacement(
        ClearPipe,
        "point",
        Dict{String, Any}(),
        function(entity)
            x, y = Ahorn.position(entity)

            entity.nodes = [(x + 32, y)]
        end
    )
)

const textures = Dict{String, String}(
    "Glass" => "glass",
    "Green" => "green"
)

Ahorn.editingOptions(entity::ClearPipe) = Dict{String, Any}(
    "texture" => textures
)

Ahorn.nodeLimits(entity::ClearPipe) = 1, -1

const exitTextureQuads = Dict{String, Dict{String, Tuple{Int, Int}}}(
    "up" => Dict{String, Tuple{Int, Int}}(
        "left" => (0, 3),
        "middle" => (1, 3),
        "right" => (2, 3)
    ),
    "down" => Dict{String, Tuple{Int, Int}}(
        "left" => (0, 5),
        "middle" => (1, 5),
        "right" => (2, 5)
    ),
    "left" => Dict{String, Tuple{Int, Int}}(
        "top" => (0, 0),
        "middle" => (0, 1),
        "bottom" => (0, 2)
    ),
    "right" => Dict{String, Tuple{Int, Int}}(
        "top" => (2, 0),
        "middle" => (2, 1),
        "bottom" => (2, 2)
    )
)

const straightTextureQuads = Dict{String, Dict{String, Tuple{Int, Int}}}(
    "vertical" => Dict{String, Tuple{Int, Int}}(
        "left" => (0, 4),
        "middle" => (1, 4),
        "right" => (2, 4)
    ),
    "horizontal" => Dict{String, Tuple{Int, Int}}(
        "top" => (1, 0),
        "middle" => (1, 1),
        "bottom" => (1, 2)
    )
)

const cornerTextureQuads = Dict{String, Dict{String, Tuple{Int, Int}}}(
    "upLeft" => Dict{String, Tuple{Int, Int}}(
        "inner" => (3, 0),
        "outer" => (5, 2),
        "horizontal_middle" => (3, 1),
        "horizontal_wall" => (3, 2),
        "vertical_middle" => (4, 0),
        "vertical_wall" => (5, 0)
    ),
    "upRight" => Dict{String, Tuple{Int, Int}}(
        "inner" => (8, 0),
        "outer" => (6, 2),
        "horizontal_middle" => (8, 1),
        "horizontal_wall" => (8, 2),
        "vertical_middle" => (7, 0),
        "vertical_wall" => (6, 0)
    ),
    "downRight" => Dict{String, Tuple{Int, Int}}(
        "inner" => (8, 5),
        "outer" => (6, 3),
        "horizontal_middle" => (8, 4),
        "horizontal_wall" => (8, 3),
        "vertical_middle" => (7, 5),
        "vertical_wall" => (6, 5)
    ),
    "downLeft" => Dict{String, Tuple{Int, Int}}(
        "inner" => (3, 5),
        "outer" => (5, 3),
        "horizontal_middle" => (3, 4),
        "horizontal_wall" => (3, 3),
        "vertical_middle" => (4, 5),
        "vertical_wall" => (5, 5)
    ),
)

function getDirection(x1, y1, x2, y2)
    if x2 < x1 && y1 == y2
        return "left"

    elseif x2 > x1 && y1 == y2
        return "right"

    elseif y2 < y1 && x1 == x2
        return "up"

    elseif y2 > y1 && x1 == x2
        return "down"
    end

    return "none"
end

function getCornerType(direction, nextDirection)
    if direction == "right" && nextDirection == "up" || direction == "down" && nextDirection == "left"
        return "upLeft"

    elseif direction == "down" && nextDirection == "right" || direction == "left" && nextDirection == "up"
        return "upRight"

    elseif direction == "left" && nextDirection == "down" || direction == "up" && nextDirection == "right"
        return "downRight"

    elseif direction == "up" && nextDirection == "left" || direction == "right" && nextDirection == "down"
        return "downLeft"
    end

    return "unknown"
end

function getLength(x1, y1, x2, y2, width, startNode, endNode)
    length = sqrt((x1 - x2)^2 + (y1 - y2)^2)

    if startNode
        length += width / 2
    end

    if endNode
        length -= width / 2
    end

    return floor(Int, length)
end

function getStraightDrawingInfo(x1, y1, x2, y2, width, length, direction, startNode, endNode)
    if direction == "up"
        return x2 - width / 2, y2 - (endNode ? 0 : width / 2), width, length

    elseif direction == "right"
        return x1 + (startNode ? 0 : width / 2), y1 - width / 2, width, length

    elseif direction == "down"
        return x1 - width / 2, y1 + (startNode ? 0 : width / 2), width, length

    elseif direction == "left"
        return x2 - (endNode ? 0 : width / 2), y2 - width / 2, width, length
    end

    return x1, y1, width, length
end

function getCornerDrawingInfo(cornerType, direction, width, length)
    verticalWallX, horizontalWallY = -1, -1
    offsetX, offsetY = 0, 0

    if cornerType == "upLeft"
        verticalWallX = width - 8
        horizontalWallY = width - 8

    elseif cornerType == "upRight"
        verticalWallX = 0
        horizontalWallY = width - 8

    elseif cornerType == "downRight"
        verticalWallX = 0
        horizontalWallY = 0

    elseif cornerType == "downLeft"
        verticalWallX = width - 8
        horizontalWallY = 0
    end

    if direction == "right"
        offsetX = length - width

    elseif direction == "down"
        offsetY = length - width
    end

    return verticalWallX, horizontalWallY, offsetX, offsetY
end

# -1 and -9 in loops used to get a inclusive-exclusive range
function renderStraightHorizontalSection(ctx::Ahorn.Cairo.CairoContext, texture, x, y, width, length, direction, startNode, endNode)
    for ox in 0:8:length - 1
        top = straightTextureQuads["horizontal"]["top"]
        middle = straightTextureQuads["horizontal"]["middle"]
        bottom = straightTextureQuads["horizontal"]["bottom"]

        # Right exit
        if direction == "left" && ox == length - 8 && startNode || direction == "right" && ox == length - 8 && endNode
            top = exitTextureQuads["right"]["top"]
            middle = exitTextureQuads["right"]["middle"]
            bottom = exitTextureQuads["right"]["bottom"]

        # Left exit
        elseif direction == "left" && ox == 0 && endNode || direction == "right" && ox == 0 && startNode
            top = exitTextureQuads["left"]["top"]
            middle = exitTextureQuads["left"]["middle"]
            bottom = exitTextureQuads["left"]["bottom"]
        end

        Ahorn.drawImage(ctx, texture, x + ox, y, top[1] * 8, top[2] * 8, 8, 8)
        Ahorn.drawImage(ctx, texture, x + ox, y + width - 8, bottom[1] * 8, bottom[2] * 8, 8, 8)

        for oy in 8:8:width - 9
            Ahorn.drawImage(ctx, texture, x + ox, y + oy, middle[1] * 8, middle[2] * 8, 8, 8)
        end
    end
end

# -1 and -9 in loops used to get a inclusive-exclusive range
function renderStraightVerticalSection(ctx::Ahorn.Cairo.CairoContext, texture, x, y, width, length, direction, startNode, endNode)
    for oy in 0:8:length - 1
        left = straightTextureQuads["vertical"]["left"]
        middle = straightTextureQuads["vertical"]["middle"]
        right = straightTextureQuads["vertical"]["right"]

        # Down exit
        if direction == "up" && oy == length - 8 && startNode || direction == "down" && oy == length - 8 && endNode
            left = exitTextureQuads["down"]["left"]
            middle = exitTextureQuads["down"]["middle"]
            right = exitTextureQuads["down"]["right"]

        # Up exit
        elseif direction == "up" && oy == 0 && endNode || direction == "down" && oy == 0 && startNode
            left = exitTextureQuads["up"]["left"]
            middle = exitTextureQuads["up"]["middle"]
            right = exitTextureQuads["up"]["right"]
        end

        Ahorn.drawImage(ctx, texture, x, y + oy, left[1] * 8, left[2] * 8, 8, 8)
        Ahorn.drawImage(ctx, texture, x + width - 8, y + oy, right[1] * 8, right[2] * 8, 8, 8)

        for ox in 8:8:width - 9
            Ahorn.drawImage(ctx, texture, x + ox, y + oy, middle[1] * 8, middle[2] * 8, 8, 8)
        end
    end
end

# -1 in loops used to get a inclusive-exclusive range
function renderCorner(ctx, texture, x, y, width, pipeWidth, length, direction, cornerType, startNode, endNode)
    if endNode || width < pipeWidth || cornerType == "unknown"
        return
    end

    nonWallQuad = straightTextureQuads["horizontal"]["middle"]
    innerCornerQuad = cornerTextureQuads[cornerType]["inner"]
    outerCornerQuad = cornerTextureQuads[cornerType]["outer"]
    verticalWallQuad = cornerTextureQuads[cornerType]["vertical_wall"]
    horizontalWallQuad = cornerTextureQuads[cornerType]["horizontal_wall"]

    verticalWallX, horizontalWallY, offsetX, offsetY = getCornerDrawingInfo(cornerType, direction, pipeWidth, length)

    for ox in 0:8:pipeWidth - 1
        for oy in 0:8:pipeWidth - 1
            quad = nonWallQuad

            onVertical = ox == verticalWallX
            onHorizontal = oy == horizontalWallY
            innerCorner = ox == abs(verticalWallX - pipeWidth + 8) && oy == abs(horizontalWallY - pipeWidth + 8)

            if innerCorner
                quad = innerCornerQuad

            elseif onVertical && onHorizontal
                quad = outerCornerQuad

            elseif onVertical
                quad = verticalWallQuad

            elseif onHorizontal
                quad = horizontalWallQuad
            end

            Ahorn.drawImage(ctx, texture, x + ox + offsetX, y + oy + offsetY, quad[1] * 8, quad[2] * 8, 8, 8)
        end
    end
end

function renderStraightSection(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Entity, texture::String, px::Int, py::Int, nx::Int, ny::Int, nnx::Int, nny::Int, width::Int, startNode::Bool, endNode::Bool)
    direction = getDirection(px, py, nx, ny)
    nextDirection = getDirection(nx, ny, nnx, nny)
    cornerType = getCornerType(direction, nextDirection)
    pipeLength = getLength(px, py, nx, ny, width, startNode, endNode)

    vertical = direction == "up" || direction == "down"
    horizontal = direction == "left" || direction == "right"

    drawX, drawY, drawWidth, drawLength = getStraightDrawingInfo(px, py, nx, ny, width, pipeLength, direction, startNode, endNode)
    straightLength = !endNode && cornerType != "unknown" ? drawLength - drawWidth : drawLength
    drawingOffset = !endNode && (direction == "up" || direction == "left") ? drawWidth : 0

    if horizontal
        renderStraightHorizontalSection(ctx, texture, drawX + drawingOffset, drawY, drawWidth, straightLength, direction, startNode, endNode)
        renderCorner(ctx, texture, drawX, drawY, drawWidth, width, pipeLength, direction, cornerType, startNode, endNode)

    elseif vertical
        renderStraightVerticalSection(ctx, texture, drawX, drawY + drawingOffset, drawWidth, straightLength, direction, startNode, endNode)
        renderCorner(ctx, texture, drawX, drawY, drawWidth, width, pipeLength, direction, cornerType, startNode, endNode)
    end
end

function renderPipe(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Entity)
    type = get(entity, "texture", "glass")
    texture = "objects/pandorasBox/clearPipe/$type/pipe"

    hasSolids =  get(entity, "hasPipeSolids", true)
    width = Int(get(entity, "pipeWidth", 32))

    px, py = Ahorn.position(entity)
    nodes = get(entity, "nodes", [])
    nodesWithPosition = vcat([(px, py)], nodes)

    if hasSolids
        for (i, node) in enumerate(nodes)
            startNode, endNode = i == 1, i == length(nodes)
            nx, ny = node
            nnx, nny = endNode ? (-1, -1) : Int.(nodes[i + 1])

            renderStraightSection(ctx, entity, texture, px, py, nx, ny, nnx, nny, width, startNode, endNode)

            px, py = nx, ny
        end
    end

    Ahorn.drawLines(ctx, nodesWithPosition, Ahorn.colors.selection_selected_fc, thickness=2)
end

function Ahorn.selection(entity::ClearPipe)
    nodes = get(entity.data, "nodes", ())
    x, y = Ahorn.position(entity)

    res = Ahorn.Rectangle[Ahorn.Rectangle(x - 6, y - 6, 12, 12)]

    for node in nodes
        nx, ny = Int.(node)

        push!(res, Ahorn.Rectangle(nx - 6, ny - 6, 12, 12))
    end

    return res
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::ClearPipe, room::Maple.Room)
    renderPipe(ctx, entity)
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
using PyPlot
using LinearAlgebra
using SparseArrays
using Arpack
using ForwardDiff
using BenchmarkTools
using Formatting
using FileIO

function Ematrix(e, v, cond)
    if cond == "stress"
        E = (e / (1 - v * v)) * [
            1 v 0
            v 1 0
            0 0 (1 - v) / 2]
    elseif cond == "strain"
        E =
            (e * (1 - v) / ((1 + v) * (1 - 2 * v))) * [
                1 v / (v - 1) 0
                v / (v - 1) 1 0
                0 0 (1 - 2 * v) / (2 - 2 * v)]
    else
        println("no planestress or planestrain condition")
        exit()
    end
    return E
end
function quadra(uv,xs,ys)
    (u,v) = uv
    N1 = u*v*(u-1)*(v-1)/4; N5 = v*(1-u*u)*(v-1)/2;
    N2 = u*v*(u+1)*(v-1)/4; N6 = u*(u+1)*(1-v*v)/2;
    N3 = u*v*(u+1)*(v+1)/4; N7 = v*(1-u*u)*(v+1)/2;
    N4 = u*v*(u-1)*(v+1)/4; N8 = u*(u-1)*(1-v*v)/2;
            N9 = (1-u*u)*(1-v*v);
    x = dot(xs,[N1 N2 N3 N4 N5 N6 N7 N8 N9])
    y = dot(ys,[N1 N2 N3 N4 N5 N6 N7 N8 N9])
    return x,y
end
function Nuv(uv)
    (u,v) = uv
    Nu1 = (2*u-1)*(v-1)*v/4;  Nv1 = (u-1)*u*(2*v-1)/4;
    Nu2 = (2*u+1)*(v-1)*v/4;  Nv2 = (u+1)*u*(2*v-1)/4;
    Nu3 = (2*u+1)*(v+1)*v/4;  Nv3 = (u+1)*u*(2*v+1)/4;
    Nu4 = (2*u-1)*(v+1)*v/4;  Nv4 = (u-1)*u*(2*v+1)/4;
    Nu5 = -u*(v-1)*v;         Nv5 = -(u*u-1)*(2*v-1)/2;
    Nu6 = -(v*v-1)*(2*u+1)/2; Nv6 = -u*(u+1)*v;
    Nu7 = -u*(v+1)*v;         Nv7 = -(u*u-1)*(2*v+1)/2;
    Nu8 = -(v*v-1)*(2*u-1)/2; Nv8 = -u*(u-1)*v;
    Nu9 = 2*u*(v*v-1);        Nv9 = 2*v*(u*u-1);

    Nu = [Nu1 Nu2 Nu3 Nu4 Nu5 Nu6 Nu7 Nu8 Nu9]
    Nv = [Nv1 Nv2 Nv3 Nv4 Nv5 Nv6 Nv7 Nv8 Nv9]

    return [Nu; Nv]
end
uv = [.5 .5]
function Jacobian(uv,xs,ys)
    J = Nuv(uv)*[xs' ys']
    return J
end
function Nxy(uv,xs,ys)
    J = Jacobian(uv,xs,ys)
    nuv = Nuv(uv)
    nxy = J\nuv
    return nxy
end
function N(uv)
    (u,v) = uv
    N1 = u*v*(u-1)*(v-1)/4; N5 = v*(1-u*u)*(v-1)/2;
    N2 = u*v*(u+1)*(v-1)/4; N6 = u*(u+1)*(1-v*v)/2;
    N3 = u*v*(u+1)*(v+1)/4; N7 = v*(1-u*u)*(v+1)/2;
    N4 = u*v*(u-1)*(v+1)/4; N8 = u*(u-1)*(1-v*v)/2;
            N9 = (1-u*u)*(1-v*v);

    N = [N1 N2 N3 N4 N5 N6 N7 N8 N9 zeros(1,9);
              zeros(1,9)           N1 N2 N3 N4 N5 N6 N7 N8 N9]
    return N
end
function LN(uv,xs,ys)
    nxy = Nxy(uv,xs,ys)
    nx = nxy[1,:]'
    ny = nxy[2,:]'
    nx = [nx zeros(1,9);
          zeros(1,9) nx]
    ny = [ny zeros(1,9);
          zeros(1,9) ny]

    LN = [1 0; 0 0; 0 1]*nx + [0 0; 0 1; 1 0]*ny

    return LN
end
function ke(E,xs,ys)
    g = [-0.949108 -0.741531 -0.405845 0.405845 0.741531 0.949108 0.0]
    w = [0.129485 0.279705 0.38183 0.38183 0.279705 0.129485 0.417959]
    kre = zeros(18,18)
    for i=1:7,j=1:7
        uv = [g[i] g[j]]
        ln = LN(uv,xs,ys)
        J = Jacobian(uv,xs,ys)
        kre = kre + w[i]*w[j]*ln'*E*ln*abs(det(J))
    end
    return kre
end

xs = [0 1 1 0 0.5 1 .5 0 .5] .+ 1000
ys = [0 0 1 1 0 .5 1 .5 .5] .+ 1000
xs = 2 .*[0 1 1 0 0.5 1 .5 0 .5] .+1000
ys = 2 .*[0 0 1 1 0 .5 1 .5 .5] .+ 1000
xs = 1 .*[0 1 1 0 0.5 1 .5 0 .5].*1 .+ 0
ys = 1 .*[0 0 1 1 0 .5 1 .5 .5].*1 .+ 0
E = Ematrix(220e9,0.5,"stress")
kre = ke(E,xs,ys)

freedofs = ones(1,18)
fixedofs = [1 8 4 10]
freedofs[fixedofs] .= 0
freedofs = [i for i=1:18 if freedofs[i] .== 1]'

global f = [1e6 0]';
F = zeros(1,18)
global g = [-0.949108 -0.741531 -0.405845 0.405845 0.741531 0.949108 0.0]
global w = [0.129485 0.279705 0.38183 0.38183 0.279705 0.129485 0.417959]
global rhs = zeros(18,1)
for i=1:7
    global rhs,f,g,w
    uv = [1 g[i]]
    nuv = Nuv(uv)
    xdot = nuv[2,:]'*xs'
    ydot = nuv[2,:]'*ys'
    road = (xdot[1]^2+ydot[1]^2)^(0.5)
    rhs = rhs + w[i]*N(uv)'*f*road
end
F[[2 3 6]]   .= rhs[[2 3 6]]

U = zeros(18,1)

Ufree = kre[freedofs,freedofs][1,:,1,:]\F[freedofs]'
U[freedofs] .= Ufree'

σ = E*LN([.5 .5],xs,ys)*U

figure(11)
for j=0:0.1:.5
    for i=-1:0.01:1

    uv = [i j]
    (x,y)=quadra(uv,xs,ys)
    s = E*LN(uv,xs,ys)*U
    scatter(i,s[1])
end
end
for k=1:1
    figure(k+300)
    for i=-1:0.1:1
        for j= -1:0.1:1
            n = Jacobian([i j],xs,ys)
            #f = n[1,k]
            f = det(n)
            #f = Nxy([i j],xs,ys)[1,k]

            scatter3D([i],[j],[f])
        end
    end
end

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
"""
    serve.jl

This is a convenience script for docs  development that
"""

using LiveServer

# Make the documentation
include("make.jl")

# Serve the documentation for development
serve(dir="build")


text:
using SortedVectors, Test

@testset "basics" begin
    sv = SortedVector([2,3,1])
    @test collect(sv) == [1, 2, 3]
    @test length(sv) == 3
    @test size(sv) == (3, )
    @test sv[1] == 1
    @test sv[2] == 2
    @test sv[3] == 3
    @test_throws BoundsError sv[4]
    @test_throws BoundsError sv[0]
    @test_throws BoundsError sv[-1]

    let T = typeof(sv)
        @test eltype(T) ≡ Int
        @test Base.IndexStyle(T) ≡ Base.IndexLinear()
    end

    let sim =  similar(sv)
        @test sim isa Vector{Int}
        @test length(sim) == length(sv)
    end

    @test_throws ArgumentError sv[1] = 7
    @test_throws ArgumentError sv[3] = -1
    @test (sv[1] = 0) == 0
    @test (sv[2] = 1) == 1
    @test (sv[3] = 8) == 8
    @test sv[:] == [0, 1, 8]

    @test parent(sv) ≡ sv.sorted_contents
    @test copy(sv) == sv

    @test_throws ArgumentError SortedVector(SortedVectors.CheckSorted(), isless, [3, 1, 2])
    @test SortedVector(SortedVectors.CheckSorted(), isless, [1, 2, 3]) == [1, 2, 3]
end

@testset "search and cut" begin
    sv = SortedVector(1:5)
    @test SortedVectors.cut.([1, 1.5, 2, 5, 6], Ref(sv)) == [0, 1, 1, 4, 5]
end


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

import Combinatorics: combinations

"""
    HLLE{NN <: AbstractNearestNeighbors, T <: Real} <: AbstractDimensionalityReduction

The `HLLE` type represents a Hessian eigenmaps model constructed for `T` type data with a help of the `NN` nearest neighbor algorithm.
"""
struct HLLE{NN <: AbstractNearestNeighbors, T <: Real} <: AbstractDimensionalityReduction
    λ::AbstractVector{T}
    proj::Projection{T}
    nearestneighbors::NN
    component::AbstractVector{Int}
end

## properties
outdim(R::HLLE) = size(R.proj, 1)
eigvals(R::HLLE) = R.λ
neighbors(R::HLLE) = R.nearestneighbors.k
vertices(R::HLLE) = R.component

## show
summary(io::IO, R::HLLE) = print(io, "Hessian Eigenmaps(outdim = $(outdim(R)), neighbors = $(neighbors(R)))")

## interface functions
"""
    fit(HLLE, data; k=12, maxoutdim=2, nntype=BruteForce)

Fit a Hessian eigenmaps model to `data`.

# Arguments
* `data`: a matrix of observations. Each column of `data` is an observation.

# Keyword arguments
* `k`: a number of nearest neighbors for construction of local subspace representation
* `maxoutdim`: a dimension of the reduced space.
* `nntype`: a nearest neighbor construction class (derived from `AbstractNearestNeighbors`)

# Examples
```julia
M = fit(HLLE, rand(3,100)) # construct Hessian eigenmaps model
R = transform(M)          # perform dimensionality reduction
```
"""
function fit(::Type{HLLE}, X::AbstractMatrix{T};
             k::Int=12, maxoutdim::Int=2, nntype=BruteForce) where {T<:Real}
    # Construct NN graph
    NN = fit(nntype, X, k)
    D, E = knn(NN, X)
    G, C = largest_component(SimpleWeightedGraph(adjmat(D,E)))
    XX = @view X[:, C]
    n = length(C)

    # Obtain tangent coordinates and develop Hessian estimator
    hs = (maxoutdim*(maxoutdim+1)) >> 1
    W = spzeros(T, hs*n, n)
    for i=1:n
        II = @view E[:,C[i]]
        # re-center points in neighborhood
        VX = view(XX, :, II)
        μ = mean(VX, dims=2)
        N = VX .- μ
        # calculate tangent coordinates
        tc = svd(N).V[:,1:maxoutdim]

        # Develop Hessian estimator
        Yi = [ones(T, k) tc zeros(T, k, hs)]
        for ii=1:maxoutdim
            Yi[:,maxoutdim+ii+1] = tc[:,ii].^2
        end
        yi = 2*(1+maxoutdim)
        for (ii,jj) in combinations(1:maxoutdim, 2)
            Yi[:, yi] = tc[:, ii] .* tc[:, jj]
            yi += 1
        end
        F = qr(Yi)
        H = transpose(F.Q[:,(end-(hs-1)):end])
        W[(1:hs).+(i-1)*hs, II] = H
    end

    # decomposition
    λ, V = decompose(transpose(W)*W, maxoutdim)
    return HLLE{nntype, T}(λ, transpose(V) .* convert(T, sqrt(n)), NN, C)
end

"""
    transform(R::LLE)

Transforms the data fitted to the Hessian eigenmaps model `R` into a reduced space representation.
"""
transform(R::HLLE) = R.proj


text:
#!/usr/bin/env julia
include("gbd17riskpred.jl")
ihd17sex(2) |> pafrate_ols |> poplot
savefig("../../images/ihd_2017_5069_gbdpaf_kv.svg")
ihd17sex(1) |> pafrate_ols |> poplot
savefig("../../images/ihd_2017_5069_gbdpaf_m.svg")

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
# calculates unlabeled attachment score
function evaluate(parser::DepParser{Unlabeled}, ss::Vector{State{Unlabeled}})
    ignore = map(w -> parser.words[w], ["''", ",", ".", ":", "``"])
    num, den = 0, 0
    for s in ss
        pred = heads(s)
        gold = map(t -> t.head, s.tokens)
        for i in 1:length(pred)
            s.tokens[i].word in ignore && continue
            den += 1
            pred[i] == gold[i] && (num += 1)
        end
    end
    uas = num / den
    @printf "UAS: %1.4f\n" uas
    uas, 0.0
end

# calculates both unlabeled & labeled attachment score
function evaluate(parser::DepParser{Labeled}, ss::Vector{State{Labeled}})
    ignore = map(w -> parser.words[w], ["''", ",", ".", ":", "``"])
    unum, lnum, den = 0, 0, 0
    for s in ss
        predheads  = heads(s)
        predlabels = labels(s)
        goldheads  = map(t -> t.head, s.tokens)
        goldlabels = map(t -> t.label, s.tokens)
        for i in 1:length(s.tokens)
            s.tokens[i].word in ignore && continue
            den += 1
            if predheads[i] == goldheads[i]
                unum += 1
                if predlabels[i] == goldlabels[i]
                    lnum += 1
                end
            end
        end
    end
    uas = unum / den
    las = lnum / den
    @printf "UAS: %1.4f\n" uas
    @printf "LAS: %1.4f\n" las
    uas, las
end

function evaluate(parser::DepParser{Unlabeled}, ss::Vector{State{Unlabeled}}, labels::Vector{Vector{Int}})
    ignore = map(w -> parser.words[w], ["''", ",", ".", ":", "``"])
    unum, lnum, den = 0, 0, 0
    for j in 1:length(ss)
        s = ss[j]
        predheads  = heads(s)
        predlabels = labels[j]
        goldheads  = map(t -> t.head, s.tokens)
        goldlabels = map(t -> t.label, s.tokens)
        for i in 1:length(s.tokens)
            s.tokens[i].word in ignore && continue
            den += 1
            if predheads[i] == goldheads[i]
                unum += 1
                if predlabels[i] == goldlabels[i]
                    lnum += 1
                end
            end
        end
    end
    uas = unum / den
    las = lnum / den
    @printf "UAS: %1.4f\n" uas
    @printf "LAS: %1.4f\n" las
    uas, las
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
module Runtests
import Example

end


text:
module Examples

using Images
using Logging

using DeformableRegistration: ImageProcessing, Transformation, Interpolation, Distance, Regularizer, Optimization
using DeformableRegistration.regOptions
using Interpolations

export registerImagesParametric, registerImagesNonParametric

function registerImagesParametric(referenceImage,templateImage, options::regOptions;
                                  affineParameters = [1.0,0,0,0,1.0,0],
                                  measureDistance = ssdDistance)

  affineParametersInitial = affineParameters
  options.parametricOnly = true

  # start multilevel registration
  for level in options.levels

    # define image level
    R = restrictResolutionToLevel(referenceImage,level)
    centeredGrid = getCellCenteredGrid(R);
    T = restrictResolutionToLevel(templateImage,level)

    # output
    Logging.info("level ",level,": [",size(R.data)[1],"]x[",size(R.data)[2],"]")

    # define objective function
    Jfunc(param;doDerivative=false,doHessian=false) =
            measureDistance(R,T,transformGridAffine(centeredGrid,param).data,
                doDerivative=doDerivative,doHessian=doHessian,options=options, centeredGrid=centeredGrid.data)

    # gauss newton method
    affineParameters = optimizeGaussNewton(Jfunc,affineParameters,affineParametersInitial,options)
  end

  return affineParameters
end

function registerImagesNonParametric(referenceImage, templateImage, options::regOptions;
                                     affineParameters = [1.0,0,0,0,1.0,0],
                                     measureDistance = ssdDistance,
                                     regularizerOperator=createCurvatureOperatorCentered)

  affineParametersInitial = affineParameters
  options.parametricOnly = false
  referenceGrid = []; deformedGrid = []; imageSize = [];

  # start multilevel registration
  for level = options.levels

  	R = restrictResolutionToLevel(referenceImage,level)
  	centeredGrid = getCellCenteredGrid(R);
  	T = restrictResolutionToLevel(templateImage,level)

    # define initial Grid
  	affineParametersInitialGrid = transformGridAffine(centeredGrid,affineParametersInitial)
  	if level==options.levels[1]
  	      # initial deformedGrid
  	      referenceGrid = transformGridAffine(centeredGrid,affineParameters)
  	      deformedGrid = referenceGrid
  	else
  	      #deformation field interpolation
  	      displacementField = deformedGrid - referenceGrid
  	      referenceGrid = transformGridAffine(centeredGrid,affineParameters)
  	      deformedGrid = interpolateDeformationField(displacementField, referenceGrid, interpolationScheme=InterpLinearFast) + referenceGrid
  	end

  	imageSize = getSize(R)
  	Logging.info("level ",level,": [",size(R.data)[1],"]x[",size(R.data)[2],"]")

  	# define objective function
  	regularizerMatrix = regularizerOperator(R.voxelsize, imageSize)
  	Jfunc(grid;doDerivative=false,doHessian=false) =
        measureDistance(R, T, grid, doDerivative=doDerivative,
            doHessian=doHessian, options=options, centeredGrid=centeredGrid.data)[1:3] +
            options.regularizerWeight * regularizer(grid-referenceGrid.data,regularizerMatrix)

  	## gauss newton method
  	deformedGrid = optimizeGaussNewton(Jfunc, deformedGrid.data, affineParametersInitialGrid.data, options)
    deformedGrid = scaledArray(deformedGrid, size(R.data), R.voxelsize, R.shift)
  end

  # return deformation field
  displacementField = deformedGrid-referenceGrid
  if options.interpolateToReferenceImage
      referenceGrid = transformGridAffine(getCellCenteredGrid(referenceImage),affineParameters)
      return interpolateDeformationField(displacementField, referenceGrid, interpolationScheme=BSpline(Linear()))
  else
      return displacementField
  end

end

end


text:
__precompile__()

module NanoTimes

export TimeDate, TimeDateZone,
       TimeZone, timezonename, astimezone, @tz_str,
       at_time, on_date, in_zone, at_zone

import Base: promote_type, promote_rule, convert, string, show, showcompact,
             (==), (!=), (<=), (<), (>), (>=), isless, isequal,
             isempty, time, min, max, minmax

import Base: start, done, next, eltype, length, size

import Base.Math: (+), (-), (*), (/),
                  fld, cld, div, rem, mod, round

using Dates

import Dates: Time, Date, DateTime,
              Year, Month, Day, Hour, Minute, Second,
              Millisecond, Microsecond, Nanosecond,
              year, month, day, hour, minute, second,
              millisecond, microsecond, nanosecond

import Dates: Period, CompoundPeriod, AbstractTime

import Dates: yearmonthday, yearmonth, monthday, dayofmonth,
              dayofweek, isleapyear, daysinmonth, daysinyear,
              dayofyear, dayname, dayabbr,
              dayofweekofmonth, daysofweekinmonth, monthname, monthabbr,
              quarterofyear, dayofquarter,
              firstdayofweek, lastdayofweek,
              firstdayofmonth, lastdayofmonth,
              firstdayofyear, lastdayofyear,
              firstdayofquarter, lastdayofquarter,
              tonext, toprev

using CompoundPeriods
import CompoundPeriods: canonical

import TimeZones: @tz_str, ZonedDateTime, TimeZone,
    localzone, astimezone, UTCOffset, FixedTimeZone, VariableTimeZone,
    all_timezones, timezone_names

import TimeZones.value # value(offset::UTCOffset) = value(offset.std + offset.dst)

const AkoTimeZone = Union{FixedTimeZone, VariableTimeZone}

include("TimeDate.jl")
include("TimeDateZone.jl")
include("timeconsts.jl")

Date(x::Date) = x
Time(x::Time) = x
DateTime(x::DateTime) = x
ZonedDateTime(x::ZonedDateTime) = x

include("periods.jl")
include("select.jl")
include("arith.jl")
include("compare.jl")
include("showstring.jl")
include("passthru.jl")

end # NanoTimes


text:
# SPDX-License-Identifier: MIT

struct GSMPhonetic{T1,T2,T3} <: AbstractModel
    f::Function
    pL::T1
    pS::T2
    model::T3
end

function kldiv(gsm::GSMPhonetic, qθ, Tθ)
    qL, qS, qλ, qσ, qπ = qθ
    μL, μS, μλ, μσ, μπ = Tθ
    retval = (
        sum(kldiv.(qL, (gsm.pL,), μL))
        + sum(kldiv.(qS, (gsm.pS,), μS))
        + sum(kldiv.(qλ, (gsm.model.pθ.pλ,), μλ))
        + sum(kldiv.(qσ, (gsm.model.pθ.pσ,), μσ))
        + sum(kldiv.(qπ, (gsm.model.pθ.pπ,), μπ))
    )
end

function getstats(gsm::GSMPhonetic, qθ)
    qL, qS, qλ, qσ, qπ = qθ
    μλ = [μ(v) for v in qλ]
    μσ = [μ(v) for v in qσ]
    μπ = [μ(v) for v in qπ]
    μ.(qL), μ.(qS), μλ, μσ, μπ
end

function unpack(gsm::GSMPhonetic, qθ, Tθ)
    qL, qS, qλ, qσ, qπ = qθ
    μL, μS, μλ, μσ, μπ = Tθ

    L = getindex.(unpack.(qL, μL), 1)
    S = getindex.(unpack.(qS, μS), 1)
    λ̂ = [vcat(reparamtrick(unpack(q, μ)), 1) for (q, μ) in zip(qλ, μλ)]
    σ̂ = [vcat(reparamtrick(unpack(q, μ)), 1) for (q, μ) in zip(qσ, μσ)]
    π̂ = [vcat(reparamtrick(unpack(q, μ)), 1) for (q, μ) in zip(qπ, μπ)]

    L, S, λ̂, σ̂, π̂
end

function _llh(gsm, labels, x::AbstractVector, uTθ)
    L, S, λ̂, σ̂, π̂ = uTθ
    (l, s, p) = labels
    P = sum(L .* λ̂[l]) + sum(S .* σ̂[s])
    ψ = P' * π̂[p]
    η = gsm.f(ψ)
    loglikelihood(gsm.model, x, η)
end

function loglikelihood(gsm::GSMPhonetic, groups, uTθ)
    sum(groups) do labels_data
        (l, s, p), data = labels_data
        L, S, λ̂, σ̂, π̂ = uTθ
        P = sum(L .* λ̂[l]) + sum(S .* σ̂[s])
        ψ = P' * π̂[p]
        η = gsm.f(ψ)
        sum(loglikelihood(gsm.model, data, η))
    end
end

function newposterior(gsm::GSMPhonetic, qθ::Tuple, ∇μ; lrate=1, clip=true)
    q₀L, q₀S, q₀λ, q₀σ, q₀π = qθ
    ∇μL, ∇μS, ∇μλ, ∇μσ, ∇μπ = ∇μ

    qL = newposterior.((gsm,), q₀L, ∇μL; lrate=lrate, clip=clip)
    qS = newposterior.((gsm,), q₀S, ∇μS; lrate=lrate, clip=clip)
    qλ = newposterior.((gsm,), q₀λ, ∇μλ; lrate=lrate, clip=clip)
    qσ = newposterior.((gsm,), q₀σ, ∇μσ; lrate=lrate, clip=clip)
    qπ = newposterior.((gsm,), q₀π, ∇μπ; lrate=lrate, clip=clip)

    qL, qS, qλ, qσ, qπ
end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA