text:
using DPC, CairoMakie, ProgressMeter, Distributions, Makie.MakieLayout

#include(joinpath(@__DIR__, "utils.jl"))
include("utils.jl")
include(joinpath(@__DIR__, "place_cells", "utils.jl"))

## Initialize the network
(net,objects) = load_network(joinpath(@__DIR__, "place_cells", "network.yaml"), weight_type=BernoulliSynapseWeight{Float64})

# Set up the place-cell populations:
populations = [
    (
        # the receptive field functions are Gaussians centered at the gridcell_centers with radii determined by gridcell_radii
        rf = v-> ismissing(v) ? 0.0 : exp(-sum((v.-μ).^2)/(2*(grid_params.r/3)^2)),
        neurons = n
    ) for (μ,n) in zip(gridcell_centers,[[objects[Symbol("i$(i)$(lpad(j,2,"0"))")] for j in 1:20] for i in 1:3])
]

## Randomly sample stochastic paths
num_paths = 50
prog = Progress(num_paths+1, 1, "Sampling paths...")
i=0
tt = LinRange(path_trange..., 50)
ProgressMeter.update!(prog, i)
paths = Vector{NamedTuple{(:x,:y),Tuple{Vector{Float64},Vector{Float64}}}}(undef, num_paths)
special_path = nothing
while i < num_paths
    # generate a single path and the response
    path = generate_stochastic_path(path_trange, path_params, generate_u₀(path_params, domain))
    spikes, logger = run_path!(net, path_trange, path, populations, λ; t_jitter, background_rate=λ_background)

    # if this path resulted in a spike, save it
    if any(x->(x.object == :n && x.event == :spikes), eachrow(logger.data))
        global i+=1
        paths[i] = (x=path.(tt,idxs=1),y=path.(tt,idxs=2))

        if i == 1
            global special_path = (;path=paths[i], spikes, logger)
        end
    end
    ProgressMeter.update!(prog, i)
end

plateau_starts_1 = filter(x->(x.object==:seg1 && x.event==:plateau_starts), special_path.logger.data)
plateau_starts_2 = filter(x->(x.object==:seg2 && x.event==:plateau_starts), special_path.logger.data)
plateau_extended_1 = filter(x->(x.object==:seg1 && x.event==:plateau_extended), special_path.logger.data)
plateau_extended_2 = filter(x->(x.object==:seg2 && x.event==:plateau_extended), special_path.logger.data)
spike_times = filter(x->(x.object==:n && x.event==:spikes), special_path.logger.data)

## Vary the speed
α_opt = 2π/6
x₀_opt = [cos(α_opt+π)*r+0.5*grid_params.xscale,sin(α_opt+π)*r+0.5*grid_params.yscale]
vs = LinRange(v_opt/100,10*v_opt,100)

prob_speed = Vector{Float64}[]
prob_rotated = Vector{Float64}[]
prob_offset = Vector{Float64}[]
αs = LinRange(0,2π,51)[1:end-1]        # movement directions over which to sweep
offsets = LinRange(-2r,2r,31)
x₀s_offset = eachrow(x₀_opt' .+ [cos(α_opt+π/2) sin(α_opt+π/2)] .* offsets) # offsets over which to sweep

x₀s_rotated = [[cos(α+π)*r+0.5*grid_params.xscale,sin(α+π)*r+0.5*grid_params.yscale] for α ∈ αs] # path starting points
for θ ∈ [8,4]
    objects[:n].θ_syn = θ
    objects[:seg1].θ_syn = θ
    objects[:seg2].θ_syn = θ

    counts = zeros(Int, length(vs))
    @showprogress 1.0 "Sweeping run-speeds ..." for (i,v) ∈ enumerate(vs)
        path = generate_straight_path((path_trange[1],path_trange[1]+2r/v), α_opt, v, x₀_opt)
        for t ∈ 1:trials
            spikes, logger = run_path!(net, path_trange, path, populations, λ; t_jitter, background_rate=λ_background)

            if any(x->(x.object == :n && x.event == :spikes), eachrow(logger.data))
                counts[i] += 1
            end
        end
    end
    push!(prob_speed, counts ./ trials)

    ## Vary the orientation
    counts = zeros(Int, length(αs))
    @showprogress 1.0 "Sweeping directions ..." for (i,(α,x₀)) ∈ enumerate(zip(αs, x₀s_rotated))
        path = generate_straight_path(path_trange, α, v_opt, x₀)
        for t ∈ 1:trials
            spikes, logger = run_path!(net, path_trange, path, populations, λ; t_jitter, background_rate=λ_background)

            if any(x->(x.object == :n && x.event == :spikes), eachrow(logger.data))
                counts[i] += 1
            end
        end
    end
    push!(prob_rotated, counts ./ trials)

    ## Vary the shift
    counts = zeros(Int, length(offsets))
    @showprogress 1.0 "Sweeping offsets ..." for (i,x₀) ∈ enumerate(x₀s_offset)
        path = generate_straight_path(path_trange, α_opt, v_opt, x₀)
        for t ∈ 1:trials
            spikes, logger = run_path!(net, path_trange, path, populations, λ; t_jitter, background_rate=λ_background)

            if any(x->(x.object == :n && x.event == :spikes), eachrow(logger.data))
                counts[i] += 1
            end
        end
    end
    push!(prob_offset, counts ./ trials)
end

## Plotting
# Set up plots
show_spines = (
    :leftspinevisible => true,
    :rightspinevisible => true,
    :bottomspinevisible => true,
    :topspinevisible => true
)

yticks,ytickformat = make_manual_ticks(collect(0:-5:-60), vcat([""],["$(grp)$(sub)" for grp in ["A","B","C"] for sub in ["₅ ", "₁₀","₁₅","₂₀"]]))

fig = Figure(resolution = (0.75textwidth, 0.75textwidth))
ax11 = fig[1,1] = Axis(fig; title="a.    Effective paths", titlealign=:left, backgroundcolor=:transparent, xlabel="x coordinate [mm]", ylabel="y coordinate [mm]",
    show_spines...)
# hidedecorations!(ax11)
gl = fig[1,2:3] = GridLayout()
ax12 = gl[1,1] = Axis(fig;
    title="b.    Activity for the highlighted path", titlealign=:left, xlabel="time [ms]",
    yticks=([-50, -30, -10], ["C","B","A"]),
    yminorticks=Makie.IntervalsBetween(20, true),
    yminorticksvisible = true, yminorgridvisible = true,
)
ax13 = gl[1,2] = Axis(fig; aspect=DataAspect(), backgroundcolor=:transparent)
hidedecorations!(ax13)
ax21 = fig[2,1] = Axis(fig; title="c.    Optimal path", titlealign=:left, backgroundcolor=:transparent, ylabel="y coordinate", show_spines...)
hidexdecorations!(ax21)
ax22 = fig[2,2] = Axis(fig; title="d.    Rotated paths", titlealign=:left, backgroundcolor=:transparent, show_spines...)
hidedecorations!(ax22)
ax23 = fig[2,3] = Axis(fig; title="e.    Translated paths", titlealign=:left, backgroundcolor=:transparent, show_spines...)
hidedecorations!(ax23)
ax31 = fig[3,1] = Axis(fig, xlabel="run speed [mm/ms]", ylabel="spike probability", )
ax32 = fig[3,2] = Axis(fig, xlabel="angle [deg]")
ax33 = fig[3,3] = Axis(fig, xlabel="offset [mm]")

rowsize!(fig.layout, 1, Aspect(1, grid_params.yscale/grid_params.xscale))
rowsize!(fig.layout, 2, Aspect(1, grid_params.yscale/grid_params.xscale))
colsize!(gl, 2, Relative(0.1))
linkaxes!(ax11, ax21)
linkaxes!(ax11, ax23)
linkyaxes!(ax31, ax32)
linkyaxes!(ax31, ax33)

# Plot the stochastically sampled paths
for path in paths[2:end]
    lines!(ax11, 1000 .*(path.x.-domain[2][1]/2), 1000 .*(path.y.-domain[2][2]/2), linewidth=2, color=RGBAf0(0.2,0.2,0.2,0.5))
    arrows!(ax11, 1000 .* [Point2f0(path.x[end],path.y[end]).-domain[2]./2], [Point2f0(path.x[end] - path.x[end-1], path.y[end] - path.y[end-1])], linewidth=2, color=RGBAf0(0.2,0.2,0.2,0.5))
end
lines!(ax11, 1000 .*(paths[1].x.-domain[2][1]/2), 1000 .*(paths[1].y.-domain[2][2]/2), linewidth=4, color=color_4)
arrows!(ax11, 1000 .* [Point2f0(paths[1].x[end],paths[1].y[end]).-domain[2]./2], [Point2f0(paths[1].x[end] - paths[1].x[end-1], paths[1].y[end] - paths[1].y[end-1])], linewidth=4, color=color_4)

# Plot the spikes, plateaus and rectangles to highlight plateau triggers
seg1 = get_trace(:seg1, special_path.logger.data)
statetrace!(ax12, [0;1000 .* seg1.t;200], [DPC.voltage_low;seg1.state;DPC.voltage_low], Dict(DPC.voltage_low=>:transparent, DPC.voltage_elevated=>color_1_25, DPC.voltage_high=>color_1_50), -20.05, 19.9)
seg2 = get_trace(:seg2, special_path.logger.data)
statetrace!(ax12, [0;1000 .* seg2.t;200], [DPC.voltage_low;seg2.state;DPC.voltage_low], Dict(DPC.voltage_low=>:transparent, DPC.voltage_elevated=>color_2_25, DPC.voltage_high=>color_2_50), -40.05, 19.9)
for (i,color) in enumerate((color_1, color_2, color_3))
    for j in 1:20
        epsp = get_trace(Symbol("syn$(i)$(lpad(j,2,"0"))"), special_path.logger.data)
        statetrace!(ax12, [0;1000 .* epsp.t;200], [0;Int.(epsp.state);0], Dict(0=>:transparent, 1=>color), -(i-1)*20 .+ -j, 0.9)
    end
end
linesegments!(ax12, 1000 .* repeat(plateau_starts_1.t, inner=2), repeat([-20, 0], outer=length(plateau_starts_1.t)), color=color_1, linewidth=2)
for t in plateau_extended_1.t
    lines!(ax12, 1000 .* [t,t], [-20.25, -0.25], color=color_1, linewidth=2)
end
linesegments!(ax12, 1000 .* repeat(plateau_starts_2.t, inner=2), repeat([-40, -20], outer=length(plateau_starts_2.t)), color=color_2, linewidth=2)
for t in plateau_extended_2.t
    lines!(ax12, 1000 .* [t,t], [-40.25, -20.25], color=color_2, linewidth=2)
end
linesegments!(ax12, 1000 .* repeat(spike_times.t, inner=2), repeat([-60, -40], outer=length(spike_times.t)), color=color_3, linewidth=2)

xlims!(ax12, 1000 .* path_trange.+[-2,6])
# Plot the neuron
plot!(ax13, objects[:n],
    branch_width=1,
    branch_length=8.0,
    color=Dict(:n=>color_3, :seg1=>color_1, :seg2=>color_2)
)

# Plot the optimal path
path_opt = generate_straight_path(path_trange, α_opt, v_opt, x₀_opt)
(path_start,path_end) = path_opt.(path_trange)
arrows!(ax21, 1000 .* Point2f0[path_start .- domain[2]./2], 1000 .* Point2f0[path_end .- path_start], linewidth=2, color=:black)

# Plot the speed-dependent spike probability
vlines!(ax31, [v_opt], linestyle=:dash, color=:gray, linewidth=2)
for (prob, linestyle) in zip(prob_speed, (:solid, :dash))
    lines!(ax31, vs, prob; color=color_4, linewidth=2, linestyle)
end

# Plot the rotated paths
pos = Point2f0[]
dir = Point2f0[]
for (α,x₀) in zip(αs, x₀s_rotated)
    path = generate_straight_path(path_trange, α, v_opt, x₀)
    (path_start,path_end) = path.(path_trange)
    push!(pos, path_start .- domain[2]./2)
    push!(dir, path_end .- path_start)
end
arrows!(ax22, 1000 .* pos, 1000 .* dir, linewidth=2, color=RGBAf0.(0.2,0.2,0.2,prob_rotated[1]))

# Plot the orientation-dependent spike probability
vlines!(ax32, [α_opt].*180/π, linestyle=:dash, color=:gray, linewidth=2)
for (prob, linestyle) in zip(prob_rotated, (:solid, :dash))
    lines!(ax32, αs.*180/π, prob; color=color_4, linewidth=2, linestyle)
end

# Plot the offset paths
pos = Point2f0[]
dir = Point2f0[]
for (p,x₀) in zip(prob_offset[1], x₀s_offset)
    path = generate_straight_path(path_trange, α_opt, v_opt, x₀)
    (path_start,path_end) = path.(path_trange)
    push!(pos, path_start .- domain[2]./2)
    push!(dir, path_end .- path_start)
end
arrows!(ax23, 1000 .* pos, 1000 .* dir, linewidth=2, color=RGBAf0.(0.2,0.2,0.2,prob_offset[1]))

# Plot the offset-dependent spike probability
vlines!(ax33, [0], linestyle=:dash, color=:gray, linewidth=2)
for (prob, linestyle) in zip(prob_offset, (:solid, :dash))
    lines!(ax33, 1000 .*offsets, prob; color=color_4, linewidth=2, linestyle)
end

xlims!(ax11, -500 * domain[2][1], 500 * domain[2][1])
ylims!(ax11, -500 * domain[2][2], 500 * domain[2][2])
ylims!(ax12, -60.5, 0.5)

save(joinpath("figures","place_cells.pdf"), fig)
save(joinpath("figures","place_cells.svg"), fig)
save(joinpath("figures","place_cells.png"), fig)
fig


text:
abstract type Word2Vec{LANG} <: EmbeddingSystem{LANG} end

function init(::Type{Word2Vec})
    register(DataDep("word2vec 300d",
    """
    Pretrained Word2Vec Word emeddings
    Website: https://code.google.com/archive/p/word2vec/
    Author: Mikolov et al.
    Year: 2013

    Pre-trained vectors trained on part of Google News dataset (about 100 billion words). The model contains 300-dimensional vectors for 3 million words and phrases.

    Paper:
        Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg Corrado, and Jeffrey Dean. Distributed Representations of Words and Phrases and their Compositionality. In Proceedings of NIPS, 2013.
    """,
    "https://s3.amazonaws.com/dl4j-distribution/GoogleNews-vectors-negative300.bin.gz",
    "21c05ae916a67a4da59b1d006903355cced7de7da1e42bff9f0504198c748da8";
    post_fetch_method=DataDeps.unpack
    ))

    push!(language_files(Word2Vec{:en}), "word2vec 300d/GoogleNews-vectors-negative300.bin")
end

function _load_embeddings(::Type{<:Word2Vec}, embedding_file, max_vocab_size, keep_words)
    local LL, indexed_words, index
    open(embedding_file,"r") do fh
        vocab_size, vector_size = parse.(Int64, split(readline(fh)))
        max_stored_vocab_size = min(max_vocab_size, vocab_size)

        indexed_words = Array{String}(max_stored_vocab_size)
        LL = Array{Float32}(vector_size, max_stored_vocab_size)

        index = 1
        @inbounds for _ in 1:vocab_size
            word = readuntil(fh, ' ')[1:end-1] #remove last space
            vector = read(fh, Float32,vector_size )

            if !contains(word, "_") && (length(keep_words)==0 || word in keep_words ) #If it isn't a phrase
                LL[:,index]=vector./norm(vector)
                indexed_words[index] = word

                index+=1
                if index>max_stored_vocab_size
                    break
                end
            end

        end
    end

    LL = LL[:,1:index-1] #throw away unused columns
    indexed_words = indexed_words[1:index-1] #throw away unused columns
    LL, indexed_words
end

ext: jl


text:
function surf_ebal(ebm::EBM, cn::Constants, Ta, Qa, Ua, Ps, SW, LW)

    Qs = qsat(false, Ps, ebm.Tsurf, cn)

    psi = ebm.gs / (ebm.gs + ebm.CH * Ua)
    if (Qs < Qa || ebm.Sice[1] > 0.0)
        psi = 1.0
    end
    Lh = cn.Ls
    if (ebm.Tsurf > cn.Tm)
        Lh = cn.Lc
    end
    rho = Ps / (cn.Rgas * Ta)
    rKH = rho * ebm.CH * Ua

    # Surface energy balance without melt
    D = Lh * Qs / (cn.Rwat * ebm.Tsurf^2.0)
    Esurf = psi * rKH * (Qs - Qa)
    ebm.Gsurf = 2.0 * ebm.ksurf * (ebm.Tsurf - ebm.Ts1) / ebm.Dz1

    ebm.Hsurf = cn.cp * rKH * (ebm.Tsurf - Ta)
    ebm.LEsrf = Lh * Esurf
    ebm.Melt = 0.0
    ebm.Rnet = (1.0 - ebm.alb) * SW + LW - cn.sb * ebm.Tsurf^4.0
    dTs = (ebm.Rnet - ebm.Hsurf - ebm.LEsrf - ebm.Gsurf) / ((cn.cp + Lh * psi * D) * rKH + 2 * ebm.ksurf / ebm.Dz1 + 4.0 * cn.sb * ebm.Tsurf^3.0)
    dE = psi * rKH * D * dTs
    dG = 2.0 * ebm.ksurf * dTs / ebm.Dz1
    dH = cn.cp * rKH * dTs
    dR = -cn.sb * ebm.Tsurf^3.0 * dTs

    # Surface melting
    if (ebm.Tsurf + dTs > cn.Tm && ebm.Sice[1] > 0.0)
        ebm.Melt = sum(ebm.Sice) / ebm.dt
        dTs = (ebm.Rnet - ebm.Hsurf - ebm.LEsrf - ebm.Gsurf - cn.Lf * ebm.Melt) / ((cn.cp + cn.Ls * psi * D) * rKH + 2.0 * ebm.ksurf / ebm.Dz1 + 4.0 * cn.sb * ebm.Tsurf^3.0)
        dE = rKH * D * dTs
        dG = 2.0 * ebm.ksurf * dTs / ebm.Dz1
        dH = cn.cp * rKH * dTs
        if (ebm.Tsurf + dTs < cn.Tm)
            Qs = qsat(false, Ps, cn.Tm, cn)
            Esurf = rKH * (Qs - Qa)
            ebm.Gsurf = 2.0 * ebm.ksurf * (cn.Tm - ebm.Ts1) / ebm.Dz1
            ebm.Hsurf = cn.cp * rKH * (cn.Tm - Ta)
            ebm.LEsrf = cn.Ls * Esurf
            ebm.Rnet = (1.0 - ebm.alb) * SW + LW - cn.sb * cn.Tm^4.0
            ebm.Melt = (ebm.Rnet - ebm.Hsurf - ebm.LEsrf - ebm.Gsurf) / cn.Lf
            ebm.Melt = max(ebm.Melt, 0.0)
            dE = 0.0
            dG = 0.0
            dH = 0.0
            dR = 0.0
            dTs = cn.Tm - ebm.Tsurf
        end
    end

    # Update surface temperature and fluxes
    ebm.Tsurf = ebm.Tsurf + dTs
    Esurf = Esurf + dE
    ebm.Gsurf = ebm.Gsurf + dG
    ebm.Hsurf = ebm.Hsurf + dH
    ebm.Rnet  = ebm.Rnet + dR
    ebm.Esnow = 0.0
    Esoil = 0.0
    if (ebm.Sice[1] > 0.0 || ebm.Tsurf < cn.Tm)
        ebm.Esnow = Esurf
        ebm.LEsrf = cn.Ls * Esurf
    else
        Esoil = Esurf
        ebm.LEsrf = cn.Lc * Esurf
    end

    return nothing

end


text:
__precompile__(false)

module ChunkManager

using MathManager
using LogManager
using FrustumManager

using StaticArrays
using JLD2
using Images

export Block
export Chunk
export ChunkNode
export reset
export update
export showAll
export getData
export createSingle
export createExample
export createLandscape
export checkInFrustum
export getActiveChilds
export getVisibleChilds
export getValidChilds
export getFilteredChilds
export inFrustum

GPU_CHUNKS = false

mutable struct Neighbours6{T} <: FieldVector{6, T}
  left::Union{T,Nothing}
  right::Union{T,Nothing}
  front::Union{T,Nothing}
  back::Union{T,Nothing}
  top::Union{T,Nothing}
  bottom::Union{T,Nothing}

  Neighbours6{T}() where T = new(nothing,nothing,nothing,nothing,nothing,nothing)
end

abstract type INode end
abstract type IBlock <: INode end
abstract type IChunk <: INode end
#RChunk = Union{Nothing,IChunk}

const NodeNeighbours6 = Neighbours6{INode}
const BlockNeighbours = Neighbours6{IBlock}

mutable struct Node6{T} <: INode
  value::Union{T,Nothing}
  next::NodeNeighbours6

  Node6{T}() where T = new(nothing,NodeNeighbours6())
  Node6{T}(value::T) where T = new(value,NodeNeighbours6())
end

mutable struct Block <: IBlock
  active::Bool
  visible::Bool
  surrounded::Bool
  typ::Int32
  flags::UInt32
  pos::Vec3f
  index::Vec3f
  sides::Array{UInt32,1}
  parent::Union{Nothing,IChunk}
  next::BlockNeighbours
  nexts::Array{Float32,1}
  id::Float32

  Block(parent::IChunk, pos=Vec3f(), nexts=Float32[], typ=0; id=0) = new(true,true,false,typ,0,pos,Vec3f(pos.x,pos.z,pos.y),resetSides(),parent,BlockNeighbours(),nexts,id) #,zeros(Mat4x4f),zeros(Mat4x4f))
end

#Base.isless{T}(a::Ref{T}, b::Ref{T}) = a.x < b.x
Base.isless(a::Block, b::Block) = a.pos < b.pos

mutable struct ChunkNode <: IChunk
  active::Bool
  visible::Bool
  pos::Vec3f
  index::Vec3f
  nodes::Array{INode,1}
  parent::Union{Nothing,IChunk}

  function ChunkNode(parent::IChunk, len::Integer)
    this = new(true,true,Vec3f(),Vec3f(),Array{INode}(undef,len),parent)
  end
end

mutable struct Chunk <: IChunk
  id::UInt32
  active::Bool
  visible::Bool
  pos::Vec3f
  index::Vec3f
  count::UInt32
  fileredCount::UInt32
  childs::Array{Block,3}
  filtered::Array{Block,1}
  nodes::Array{ChunkNode,1}
  parent::Union{Nothing,IChunk}

  Chunk(;id=0,pos=Vec3f()) = new(id,true,true,pos,Vec3f(),0,0,Array{Block}(undef,0,0,0),Block[],Array{ChunkNode}(undef,0),nothing)
end

function init(this::Chunk, size::Integer)
    this.count = size^3
    this.childs = Array{Block}(undef,size,size,size)
    createBlocks(this)

    #if GPU_CHUNKS linkBlockNeighbours(this); return end
    #createChunkNodes(this)
    #linkBlockNeighbours(this)
end

function update(this::Chunk;unseen=true)
  #if GPU_CHUNKS && unseen hideUnseen(this); return end
  #updateChunkNodePos(this)
  #if unseen hideUnseen(this) end
  #r=create_spiral3D_searchlist(this)
end

function reset(this::Chunk; size=64)
  clean(this)
  init(this,size)
end

function clean(this::Union{Nothing,Chunk})
  if this == nothing return end
  this.childs = Array{Block}(undef,0,0,0)
  this.nodes = Array{ChunkNode}(undef,0)
  this.filtered = Block[]
  GC.gc() # force garbage collection, free memory
end

getPosition(this::INode) = this.pos

function isActive(this::INode)
  if !this.active return false end
  node=this
  #while (parent = parent.parent) != nothing if !parent.active return false end end
  true
end

function isVisible(this::INode)
  if !this.visible return false end
  node=this.parent
  #while node != nothing if !node.visible return false else node=node.parent; end end
  true
end

setActive(this::INode, active::Bool) = this.active = active
setVisible(this::INode, visible::Bool) = this.visible = visible

isType(this::Block, typ) = this.typ == typ

isValid(this::Chunk) = this.count > 0 && this.fileredCount > 0

isSeen(this::Block) = isActive(this) && isVisible(this) && this.typ > 0

i=0
i+=1; const LEFT_SIDE = i #0x1
i+=1; const RIGHT_SIDE = i #0x2
i+=1; const TOP_SIDE = i #0x4
i+=1; const BOTTOM_SIDE = i #0x8
i+=1; const FRONT_SIDE = i #0x10
i+=1; const BACK_SIDE = i #0x20
const MAX_SIDES = i

resetSides() = fill(UInt32(1),MAX_SIDES)
resetSides(this::Block) = this.sides=resetSides()

#=
function createChunkNodes_recurive(this::Chunk, depth=size(this.childs)[1], node::IChunk=this, index=0, pos=ones(UInt,3); len=8)
  if index > 0 node = ChunkNode(node, len) end
  depth = round(UInt,depth/2)
  if depth > 1 for i=1:len node.nodes[i] = createChunkNodes(this, depth, node, i, getIndex(i, depth, pos)) end
  else
    for i=1:len
      x,y,z = getIndex(i, depth, pos)
      node.nodes[i] = b = this.childs[x,z,y]
      this.parent=node
    end
  end
  node
end
=#

function getIndex(i, depth, vindex)
  x=(i == 1 || i == 3 || i == 5 || i == 7) ? 0 : 1
  z=(i == 1 || i == 2 || i == 5 || i == 6) ? 0 : 1
  y=(i >= 1 && i <= 4) ? 0 : 1
  vindex + [x,z,y] * depth
end

function createChunkNodes(this::Chunk)
  SIZE=size(this.childs)[1]
  this.nodes=Array{ChunkNode}(undef,8)
  depths=zeros(UInt,round(UInt,log(2,SIZE)))
  max_depth = length(depths)
  scale=[2^(max_depth-i) for i=1:max_depth] #[32,16,8,4,2,1]
  depth=1
  count = 0
  len=length(this.childs)
  default_index=ones(UInt,3)
  node = this

  while depth > 0 && count <= len
    depths[depth]+=1
    index = depths[depth]
    if index > 8
      depths[depth]=0
      depth-=1
      node=node.parent
      continue
    end

    x,z,y = cindex = getIndex(index, scale[depth], depth == 1 ? default_index : node.index)
    child = depth < max_depth ? ChunkNode(node, 8) : this.childs[UInt(x),UInt(z),UInt(y)]
    child.index = cindex

    node.nodes[index] = child
    child.parent=node

    if depth < max_depth depth+=1; node=child; else count+=1;end
  end
end

function itranslate(this::INode, SIZE::Integer)
  x,z,y = this.index
  X, Y, Z = SIZE, SIZE, SIZE
  DIST = Vec3f(2,2,2); #Vec3f(3,5,3) #r = 1f0/30 #0.005f0
  START = Vec3f(-(X*DIST.x) / 2.0f0, -(Y*DIST.y) / 2.0f0, (Z*DIST.z) / 2.0f0)
  Vec3f(START.x+(x-1)*DIST.x, START.y+(y-1)*DIST.y, START.z-(z-1)*DIST.z)
end

function itranslate(index::Vec3f)::Vec3f
  DIST = 0.25
  STARTDIST = (128*DIST) / 2f0
  Vec3f(-1,-1,1)*STARTDIST+index*Vec3f(1,1,-1)*DIST-Vec3f(0,STARTDIST,0)
end

function updateChunkNodePos(this::Chunk)
  SIZE=size(this.childs)[1]
  depths=zeros(UInt,round(UInt,log(2,SIZE)))
  max_depth = length(depths)
  depth=1
  count = 0
  len=length(this.childs)
  node = this

  while depth > 0 && count <= len
    depths[depth]+=1
    index = depths[depth]
    if index > 8
      node.index = Vec3f();
      for child in node.nodes node.index += child.index end
      node.index /= 8
      node.pos = itranslate(node, SIZE)
      depths[depth]=0
      depth-=1
      node=node.parent
      continue
    end
    child = node.nodes[index]
    if depth < max_depth depth+=1; node=child; else count+=1; child.pos = itranslate(child, SIZE); end
  end
end

inFrustum(this::Chunk, fstm::Frustum) = checkCube(fstm, itranslate(this.pos), Vec3f(128,128,128))[1] != :FRUSTUM_OUTSIDE #checkSphere(fstm, itranslate(this.pos), 128.0)

function checkInFrustum(this::Chunk, fstm::Frustum)
  depths=[0,0,0,0,0,0]
  scale=[32,16,8,4,2,1]
  max_depth = length(depths)
  depth=1
  count = 0
  len=length(this.childs)
  node = this
  list = Block[]

  #setFilteredChilds(this,getValidChilds(this))
  #return true

  #for b in this.childs setVisible(b, false) end

  #println("done.")
  #open("nodes.txt", "w") do f
  while depth > 0 && count <= len
    depths[depth]+=1
    index = depths[depth]

    if index <= 1
      #write(f, string(repeat(" ",depth-1),"{ ",scale[depth]," : ",  depth>1 ? depths[depth-1] : 0, "\n"))
    elseif index > 8
      #write(f, string(repeat(" ",depth-1),"} ", scale[depth]," : ", depth>1 ? depths[depth-1] : 0, "\n"))
      depths[depth]=0;
      depth-=1;
      node=node.parent
      continue
    end
    child = node.nodes[index]

    #visible = true
    #if typeof(child) == Block
    p = getPosition(child)
    #result = checkCube(fstm, p, Vec3f([1,1,1] * scale[depth]))[1]
    result = checkSphere(fstm, p, scale[depth] * 1.5)[1]
    visible = result != :FRUSTUM_OUTSIDE
    #write(f, string(repeat(" ",depth),depth >= max_depth ? "L" : "C",index,": ",p[1],", ",p[2],", ",p[3],visible ? "" : " ---> R ", "\n"))
    setVisible(child, visible)
    #end

     #result == :FRUSTUM_INTERSECT
    if depth >= max_depth
      count+=1;
      if isValid(child)
        push!(list, child)
        #=write(f, string(count,"\n"));=#
      end
    elseif visible
      #if depth == 5 && !visible continue end
      depth+=1; node=child
    end
  end
  #end

  setFilteredChilds(this,list)
  #setFilteredChilds(this,getValidChilds(this))
end

function createBlocks(this::Chunk)
  X, Y, Z = size(this.childs); XZ = X*Z; i=0;
  for y=1:Y; for z=1:Z; for x=1:X; i+=1
    #i = (y-1)*XZ+(z-1)*X+(x-1)
    this.childs[x,z,y] = Block(this,Vec3f(x,y,z), Float32[
    x>1 ? i-1 : 0,
    x<X ? i+1 : 0,
    y>1 ? i+XZ : 0,
    y<Y ? i-XZ : 0,
    z>1 ? i+X : 0,
    z<Z ? i-X : 0];id=i)
  end; end; end

  open("chunks.txt", "w") do f
    for b in this.childs
      write(f, string(b.id," : ",b.index.x,", ",b.index.y,", ",b.index.z,"\n"))
    end
  end
end

function linkBlockNeighbours(this::Chunk)
  a = this.childs
  X, Y, Z = size(this.childs)
  for y=1:Y; for z=1:Z; for x=1:X;
    pos = (x,y,z)
    b = getBlock(this,pos)
    b.next.left = getBlock(this,pos,LEFT_SIDE)
    b.next.right = getBlock(this,pos,RIGHT_SIDE)
    b.next.bottom = getBlock(this,pos,BOTTOM_SIDE)
    b.next.top = getBlock(this,pos,TOP_SIDE)
    b.next.back = getBlock(this,pos,BACK_SIDE)
    b.next.front = getBlock(this,pos,FRONT_SIDE)
  end; end; end
end

function getBlock(this::Chunk, pos::Tuple{Integer,Integer,Integer}, side=0)
  a = this.childs
  X, Y, Z = size(this.childs)
  (x,y,z) = pos
  b = nothing
  if side == LEFT_SIDE
    b=x>1 ? a[x-1,z,y] : nothing
  elseif side == RIGHT_SIDE
    b=x<X ? a[x+1,z,y] : nothing
  elseif side == BOTTOM_SIDE
    b=y>1 ? a[x,z,y-1] : nothing
  elseif side == TOP_SIDE
    b=y<Y ? a[x,z,y+1] : nothing
  elseif side == BACK_SIDE
    b=z>1 ? a[x,z-1,y] : nothing
  elseif side == FRONT_SIDE
    b=z<Z ? a[x,z+1,y] : nothing
  else
    b=x>=1&&x<=X&&y>=1&&y<=Y&&z>=1&&z<=Z ? a[x,z,y] : nothing
  end
  b
end

function hideUnseen(this::Chunk)
  for b in this.childs
    if !isSeen(b) continue end # skip invisible

    blocked=0

    if (next=b.next.left) != nothing
      r=b.sides[LEFT_SIDE]=!isSeen(next)
      if !r blocked+=1 end
    end

    if (next=b.next.right) != nothing
      r=b.sides[RIGHT_SIDE]=!isSeen(next)
      if !r blocked+=1 end
    end

    if (next=b.next.bottom) != nothing
      r=b.sides[BOTTOM_SIDE]=!isSeen(next)
      if !r blocked+=1 end
    end

    if (next=b.next.top) != nothing
      r=b.sides[TOP_SIDE]=!isSeen(next)
      if !r blocked+=1 end
    end

    if (next=b.next.back) != nothing
      r=b.sides[BACK_SIDE]=!isSeen(next)
      if !r blocked+=1 end
    end

    if (next=b.next.front) != nothing
      r=b.sides[FRONT_SIDE]=!isSeen(next)
      if !r blocked+=1 end
    end

    setSurrounded(b,blocked >= 6)
  end
end

function setFlag(this::Block, flag::Unsigned, add::Bool)
  if add this.flags |= flag
  else this.flags ⊻= flag #⊻ = xor and is not shown on windows...
  end
end

isSurrounded(this::Block) = this.surrounded
setSurrounded(this::Block, surrounded::Bool) = this.surrounded = surrounded

isValid(this::Block) = isActive(this) && isVisible(this) && !isSurrounded(this) && this.typ > 0

hideType(this::Chunk, typ::Integer) = for b in this.childs; if b.typ == typ; setVisible(b,false); end; end
removeType(this::Chunk, typ::Integer) = for b in this.childs; if b.typ == typ; setActive(b,false); end; end

function showAll(this::Chunk)
  #if GPU_CHUNKS
  setFilteredChilds(this,filter(b->true,this.childs))
  #=
  return end
  for b in this.childs
    setVisible(b,true)
    #setSurrounded(b,false)
    #resetSides(b)
  end
  setFilteredChilds(this,getValidChilds(this))
  =#
end

setFilteredChilds(this::Chunk, r::Array{Block,1}) = begin this.filtered = r; this.fileredCount = length(r); r end
getFilteredChilds(this::Chunk) = this.filtered

getActiveChilds(this::Chunk) = filter(b->isActive(b),this.childs)
getVisibleChilds(this::Chunk) = filter(b->isVisible(b),this.childs)
getValidChilds(this::Chunk) = filter(b->isValid(b),this.childs)

function getData(this::Block)
  # get visible sides
  i=0; sides=0
  for side in this.sides
    if side > 0 sides |= (0x1 << i) end
    i+=1
  end

  #if this.sides[LEFT_SIDE] > 0 sides |= 0x1 end
  #if this.sides[RIGHT_SIDE] > 0 sides |= 0x2 end
  #if this.sides[TOP_SIDE] > 0 sides |= 0x4 end
  #if this.sides[BOTTOM_SIDE] > 0 sides |= 0x8 end
  #if this.sides[FRONT_SIDE] > 0 sides |= 0x10 end
  #if this.sides[BACK_SIDE] > 0 sides |= 0x20 end

  #[n-1 for n in this.nexts]...,
  SVector(Float32[this.index.x-1, this.index.z-1, this.index.y-1, this.typ, sides]...)
end

getPos(this::Block) = SVector(Float32[this.pos...]...)

getData(this::Chunk) = isValid(this) ? vec((b->getData(b)).(getFilteredChilds(this))) : Float32[]
getPos(this::Chunk) = isValid(this) ? vec((b->getPos(b)).(getFilteredChilds(this))) : Float32[]

#fill copies references
#refblocks = Ptr{Nothing}[]
#for b in this.childs push!(refblocks, pointer_from_objref(b)) end

function createSingle(this::Chunk)
  b=this.childs[32,22,32]
  b.typ=1
end

function createExample(this::Chunk)
  X, Y, Z = sz = size(this.childs)
  #DIST = Vec3f(2,2,2); #Vec3f(3,5,3) #r = 1f0/30 #0.005f0
  #START = Vec3f(-(X*DIST.x) / 2.0f0, -(Y*DIST.y) / 2.0f0, (Z*DIST.z) / 2.0f0)
  #translate(x,y,z) = Vec3f(START.x+(x-1)*DIST.x, START.y+(y-1)*DIST.y, START.z-(z-1)*DIST.z)

  x=1; y=1; z=1;

  for i=1:length(sz)
    b = this.childs[i];
    b.typ=rand(1:16)
    #model = Mat4x4f(translation(Array(b.pos)))
    #push!(refblocks, pointer_from_objref(b))

    x += 1; if x > X
      z += 1; x=1;
      if z > Z y += 1; z=1;
        if y > Y break; end
      end
    end
  end

  h=round(UInt,X/2)
  b = this.childs[h,h,h]
  b.typ=7
end

function createLandscape(this::Chunk)
  if true return nothing end
  #Texture* heightTexture = m_pRenderer->GetTexture(m_pChunkManager->GetHeightMapTexture());

  rgbImage=Images.load("heightmap.jpg")
  (imgwidth, imgheight) = size(rgbImage)

  grayImage = Gray.(rgbImage)
  imgMatrix = reinterpret.(channelview(grayImage)) #convert(Array{Int32},raw(grayImage))

  X, Y, Z = size(this.childs)
  w = imgwidth/X
  h = imgheight/Y

  #DIST = Vec3f(2,2,2)
  #START = Vec3f(-(X*DIST.x) / 2.0f0, -(Y*DIST.y) / 2.0f0, (Z*DIST.z) / 2.0f0)
  #translate(x,y,z) = Vec3f(START.x+x*DIST.x, START.y+y*DIST.y, START.z-z*DIST.z)

  for z=1:Z
    for x=1:X
      # Use the height map texture to get the height value of x, z
      height = (imgMatrix[UInt32(trunc(z*h)), UInt32(trunc(x*w))] / 0xFF)* 1.25 * Y

      level_air = height * 0.95
      level_grass = height * 0.9
      level_dirt = height * 0.8
      level_stonebricks = height * 0.7
      level_stone = height * 0.4
      level_lava = height * 0

      for y=1:Y
        if y > height break end
        #id = y*X*Y+z*Z+x
        b = this.childs[x,z,y]
        #if y >= level_air b.typ = 0 # air or nothing
        if y >= level_grass  b.typ = 2 #grass
        elseif y >= level_dirt b.typ = 1 #dirt
        elseif y >= level_stonebricks b.typ = 5 #stonebricks
        elseif y >= level_stone b.typ = 4 #stone
        elseif y >= level_lava b.typ = 15 #lava
        end
      end
    end
  end

  h=round(UInt, X/2)
  b = this.childs[h,h,h]
  #b.pos=Vec3f(0,0,0)
  b.typ=7

  this.childs[1,1,1].typ=0
end

function saveChunk(this::Chunk, file::String)
  #childs=getValidChilds(this)
  #childs=filter(b->(b.pos, b.typ),childs)
  r=[(b.pos, b.typ) for b in this.childs]
  @save file*".jld2" r
end

function loadChunk(this::Chunk, file::String)
  @load file*".jld2" r
  i=1; for c in this.childs; b=r[i]; c.pos=b[1]; c.typ=b[2]; i+=1; end
end

"""
  Spiral 3D Search List

  clock wise neighbour search around a center point (shell by shell search)
  first and last layer (z-axis) full 2D spiral search (whole layer is shell)
  layers in between: search outside (next shell) and already searched area (core) will be skipped
"""
function create_spiral3D_searchlist(this::Chunk)
  println("Spiral 3D Search List")

  Z, Y, X = sz = size(this.childs)
  MAX = length(sz)
  r = Array{Int}(undef, MAX)

  file="spiral3D_searchlist_chunk_"*string(X)*"x"*string(Y)*"x"*string(Z)*".jld2"
  if isfile(file) @load file r; return r; end

  XH, YH, ZH = X/2, Y/2, Z/2

  lastElement = true
  endOflayer = false
  endOflBox = false
  turn=false
  log=false

  xyz=[0,0,0]
  count=[0,0]
  element=1; j=0;
  xt=0; yt=0;
  layerElement = 0
  layer = 0
  edge=0
  layerCount = 0
  maxlength = 0
  layerElementCount = 0
  layerElementCountPrev =0

  r[element]=ZH*YH*XH;

  println("create search coordinates...")
  open(file*".txt", "w") do f
  while true
    if lastElement
      j += 1
      if j >= ZH break end
      edge=0
      layer=1
      endOflayer = false
      endOflBox = false
      layerElement = 0

      dimMAX = 1 + j*2
      layerCount = dimMAX
      layerElementCountPrev = layerElementCount
      layerElementCount = dimMAX^2 # X*Y

      xyz=[0,0,j]
      xt=1; yt=-1;
      count=[1,1]
      turn=true
      lastElement = false
      #print(".")
    end

    if !lastElement
      layerElement += 1
      element+=1
      if element > MAX break end

      notFirst = !(layer == 1 && layerElement == 1)
      firstLayer = layer <= 1
      lastLayer = layer >= layerCount
      max = layerElementCount-(firstLayer || lastLayer ? (lastLayer ? 1 : 0) : layerElementCountPrev)
      endOflayer = (layerElement > max)
      border = [abs(xyz[1]) >= count[1],abs(xyz[2]) >= count[2]]

      if notFirst
        switch = [turn && border[1],!turn && border[2]]
        if switch[1] xt=-xt; elseif switch[2] yt=-yt; end
        if switch[1] || switch[2] turn=!turn; end
      end

      if border[1] && border[2] edge+=1 end

      if endOflayer
        layerElement=1
        layer+=1
        edge=0
        firstLayer = layer <= 1
        lastLayer = layer >= layerCount
        lastElement = layer > layerCount
        count=[j,j]

        if lastLayer
          xt=1
          xyz=[0,0,xyz[3]]
          turn=true
          count=[1,1]
        end

        if !lastElement xyz[3]-=1; end
        if log print("\$"); if !lastElement print("\n////// |") end end
      end

      if edge > 3
        if log print("#") end
        edge=0
        xt=-xt
        count+=[1,1]
        turn = true
      end

      endOfCube = lastLayer && endOflayer
      firstLayerElement = layerElement == 1
      lastLayerElement = endOflayer || layerElement == max
      lastElement = layer > layerCount

      if notFirst; if turn xyz[1]+=xt; else xyz[2]+=yt; end; end
      if lastElement xyz=[0,0,xyz[3]]; end

      x,y,z = (xyz[2]+XH), (xyz[1]+YH), (xyz[3]+ZH)
      r[element] = index = z*Y*X + y*X + x
      write(f, string([z,y,x]," = ", index, "\n"));

      if log print(x>=0 ? " " : "",x, "" , y>=0 ? " " : "",y,"",z>=0 ? " " : "", z," |") end
    end
  end
  end
  println("done.")
  @save file r
  r
end

const FrustumNode = Node6{Tuple{Block,Any}}

printList(list,n=1) = begin
  n = round(n)
  if n<1 error("n < 1 is invalid!") end
  println(typeof(list)," with ",length(list)," entries shows each ",n,n==1 ? "st" : (n==2 ? "nd" : (n==3 ? "rd" : "th"))," element:")
  i=0; for (k,v) in list; i+=1; if i == n i=0; println(" ",k) end; end
  println()

  #a = collect(keys(list))[1:X:end]
  #Base.showarray(STDOUT,a,false)
  #show(STDOUT, "text/plain", a)
  #display(a)
  #whos()
end

function checkInFrustum2(this::Chunk, fstm::Frustum)
  pos = getNearPosition(fstm)
  list = DataStructures.SortedDict()

  #open("frustum.txt", "w") do f
    for b in this.childs
      result = checkSphere(fstm, b.pos, 1.5)
      clist = result[2]
      visible = result[1] != :FRUSTUM_OUTSIDE
      setVisible(b,visible)

      # sort
      if visible
        (x,y,z) = (clist[:X],clist[:Y],clist[:Z])
        write(f, string("[",x,", ",y,", ",z,"]")*"\n")

        rx=round(x)
        ry=round(y)
        rz=round(z)

        if rx % 2 != 0
          rx+=1;
        end

        k = (rx,ry)

        if (abs(x) <= 1 && abs(y) <= 1 && abs(z) <= 1) || z < 0
          setVisible(b,false)
        else
          kv=z=>b
          if haskey(list,k)

            p = first(list[k])
            list[k][z] = b
            q = first(list[k])

            if p[1] == q[1] kv=p[1]=>p[2] else b=p[2] end
            setVisible(b,false)
          end
          list[k] = DataStructures.SortedDict(kv)
        end
      end
    end
  #end
end

function checkInFrustum3(this::Chunk, fstm::Frustum)
  println("checkInFrustum")
  pos = getNearPosition(fstm)
  list = DataStructures.SortedDict()

  proceed1 = function(b::Block)
    result = checkSphere(fstm, b.pos, 1.5)
    visible = result[1] != :FRUSTUM_OUTSIDE
    setVisible(b,visible)

    #=
    clist = result[2]
    if visible && false
      current = FrustumNode((b,clist))
      parent = nothing
      if parent == nothing parent = current
      else
        nlist = parent.value[2]
        while true
          lt = clist[:FRUSTUM_LEFT] > clist[:FRUSTUM_TOP]
          l = nlist[:FRUSTUM_LEFT] > clist[:FRUSTUM_LEFT]
          r = nlist[:FRUSTUM_RIGHT] < clist[:FRUSTUM_RIGHT]
          t = nlist[:FRUSTUM_TOP] > clist[:FRUSTUM_TOP]
          b = nlist[:FRUSTUM_BOTTOM] < clist[:FRUSTUM_BOTTOM]
          n = nlist[:FRUSTUM_NEAR] < clist[:FRUSTUM_NEAR]
          f = nlist[:FRUSTUM_FAR] > clist[:FRUSTUM_FAR]

          #lr = l>r ? l : r
          #tb = t>b ? t : b
          #nf = n<f ? n : f

          if lt
            if l
              next = parent.next.left
              if next == nothing parent.next.left = current
              else parent = next
              end
            elseif r
              next = parent.next.right
              if next == nothing parent.next.right = current
              else parent = next
              end
            end
          else
            if t
              next = parent.next.top
              if next == nothing parent.next.top = current
              else parent = next
              end
            elseif b
              next = parent.next.bottom
              if next == nothing parent.next.bottom = current
              else parent = next
              end
            end
          end
        end

      end
    end
    =#
    return visible
  end

  proceed2 = function(x,y,z)
      if !isassigned(this.childs, x,y,z) return false end
      proceed(this.childs[x,y,z])
  end
  for b in this.childs proceed1(b) end

  #printList(list,length(list)/10)

  #for (k,z) in list
    #p = k1 + Vec3f(0,0,0)
    #l = k1 + Vec3f(-1,-1,-1)
    #r = k1 + Vec3f(1,1,1)

    #skip = true
    #for (k2,b) in z
      #if skip skip=false; continue end
      #if k2 == nothing continue end
      #if p.z - 10 >= k2.z continue end
      #if p.y + 1 >= k2.y continue end
      #if p.y - 1 <= k2.y continue end
      #if l.x <= k2.x continue end
      #if r.x >= k2.x continue end
      #setVisible(b,false)
    #end
  #end

end

end #ChunkManager

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
using NLPModels
using LinearAlgebra

abstract type AbstractGNEP{T, S} end

mutable struct GNEP{T, S, NLP <: AbstractNLPModel{T, S}} <: AbstractGNEP{T, S}
  N::Int64 #number of player
  n::Array{Int64, 1} #array of integer (size of the player strategy)
  player::Array{NLP, 1} #array of nlpmodels

  x0::S
  sol::Function #fonction returning a boolean if x is a solution
  class::String
  gnsc::Bool #true if it is a GNSC
  proj::Bool #true if the projection function is furnished
  fnproj::Function #projection over the feasible set

  function GNEP(
    N::Int64,
    n::Array{Int64, 1},
    player::Array{NLP, 1};
    x0::S = ones(sum(n)),
    sol::Function = x -> true,
    class::String = "GNEP",
    gnsc::Bool = false,
    proj::Bool = false,
    fnproj::Function = x -> x,
  ) where {NLP, S}
    if length(n) != N
      throw(error("wrong size"))
    end
    if length(x0) != sum(n)
      throw(error("wrong size"))
    end
    if length(player) != N
      throw(error("wrong size"))
    end

    return new{eltype(x0), S, NLP}(N, n, player, x0, sol, class, gnsc, proj, fnproj)
  end
end

function pseudo_grad(gnep::AbstractGNEP{T, S}, x::S) where {T, S}
  g = zeros(T, sum(gnep.n))
  for i = 1:(gnep.N)
    nI = (1 + sum(gnep.n[1:(i - 1)])):sum(gnep.n[1:i])
    g[nI] = grad(gnep.player[i], x)[nI]
  end
  return g
end

function pseudo_cons(gnep::AbstractGNEP{T, S}, x::S) where {T, S}
  cx = zeros(T, nb_cons(gnep))
  c = 0
  for i = 1:(gnep.N)
    ncon = gnep.player[i].meta.ncon
    nI = (c + 1):(c + ncon)
    cx[nI] = cons(gnep.player[i], x)
    c += ncon
  end
  return cx
end

function pseudo_res(gnep::AbstractGNEP{T, S}, x::S) where {T, S}
  cx = zeros(T, nb_cons(gnep))
  c = 0
  for i = 1:(gnep.N)
    ncon = gnep.player[i].meta.ncon
    nI = (c + 1):(c + ncon)
    cx[nI] = cons(gnep.player[i], x) - gnep.player[i].meta.ucon
    c += ncon
  end
  return cx
end

function pseudo_vio(gnep::AbstractGNEP{T, S}, x::S) where {T, S}
  cx = zeros(T, nb_cons(gnep))
  c = 0
  for i = 1:(gnep.N)
    ncon = gnep.player[i].meta.ncon
    nI = (c + 1):(c + ncon)
    cx[nI] =
      max.(cons(gnep.player[i], x) - gnep.player[i].meta.ucon, 0) +
      max.(-cons(gnep.player[i], x) + gnep.player[i].meta.lcon, 0)
    c += ncon
  end
  return cx
end

function pseudo_hess(gnep::AbstractGNEP{T, S}, x::S) where {T, S}
  H = zeros(T, sum(gnep.n), sum(gnep.n))
  for i = 1:(gnep.N)
    nI = (1 + sum(gnep.n[1:(i - 1)])):sum(gnep.n[1:i])
    H[nI, nI] = hess(gnep.player[i], x)[nI, nI]
  end
  return H
end

function nb_cons(gnep::AbstractGNEP)
  count = 0
  for i = 1:(gnep.N)
    count += gnep.player[i].meta.ncon
  end
  return count
end

function pseudo_jac(gnep::AbstractGNEP{T, S}, x::S) where {T, S}
  J = zeros(T, nb_cons(gnep), sum(gnep.n))
  c = 0
  for i = 1:(gnep.N)
    nI = (1 + sum(gnep.n[1:(i - 1)])):sum(gnep.n[1:i])
    ncon = gnep.player[i].meta.ncon
    nC = (c + 1):(c + ncon)
    J[nC, nI] = jac(gnep.player[i], x)[:, nI]
    c += ncon
  end
  return J
end

# Lagrange Multiplier of the KKT-GNEP
# (without the bounds constraints)
# (without the signs constraints)
#
function lagrange_multiplier(gnep::AbstractGNEP{T, S}, x::S; prec::Float64 = 1e-6) where {T, S}
  lambda = zeros(T, nb_cons(gnep) + sum(gnep.n))
  c = 0

  for i = 1:(gnep.N)
    nI = (1 + sum(gnep.n[1:(i - 1)])):sum(gnep.n[1:i])
    nn = length(nI)
    ncon = gnep.player[i].meta.ncon

    cx = cons(gnep.player[i], x)
    ncP = findall(
      min.(abs.(cx - gnep.player[i].meta.ucon), abs.(gnep.player[i].meta.lcon - cx)) .<= prec,
    )
    ncB = findall(
      min.(
        abs.(x[nI] - gnep.player[i].meta.uvar[nI]),
        abs.(gnep.player[i].meta.lvar[nI] - x[nI]),
      ) .<= prec,
    )
    J = vcat(jac(gnep.player[i], x)[ncP, nI], Matrix(1.0I, nn, nn)[ncB, :])
    lambda[c .+ union(ncP, ncon .+ ncB)] =
      -LinearAlgebra.pinv(convert(Array, J')) * grad(gnep.player[i], x)[nI]
    c += ncon + nn
  end

  return lambda
end

#bounds problem
#function lagrange_multiplier_bool(gnep :: GNEP, x :: Vector; prec :: Float64 = 1e-6)
# lambda = lagrange_multiplier(gnep, x, prec = prec)
# return (norm(pseudo_jac(gnep,x)'*lambda+pseudo_grad(gnep,x)) <= prec)
#end

function lagrange_multiplier_bool(gnep::AbstractGNEP{T, S}, x::S; prec::Float64 = 1e-6) where {T, S}
  lambda = lagrange_multiplier(gnep, x, prec = prec)
  cl = 0
  res = Inf
  for i = 1:(gnep.N)
    nI = (1 + sum(gnep.n[1:(i - 1)])):sum(gnep.n[1:i])
    nn = length(nI)
    ncon = gnep.player[i].meta.ncon
    nL = (cl + 1):(cl + ncon + nn)
    J = vcat(jac(gnep.player[i], x)[:, nI], Matrix(1.0I, nn, nn))
    res = min(res, norm(J' * lambda[nL] + grad(gnep.player[i], x)[nI]))
    cl += ncon + nn
  end
  return res <= prec
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct CheckpointableObjectGraph_CheckpointableObject_ObjectReference <: ProtoType
    node_id::Int32
    local_name::AbstractString
    CheckpointableObjectGraph_CheckpointableObject_ObjectReference(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct CheckpointableObjectGraph_CheckpointableObject_ObjectReference

mutable struct CheckpointableObjectGraph_CheckpointableObject_SerializedTensor <: ProtoType
    name::AbstractString
    full_name::AbstractString
    checkpoint_key::AbstractString
    optional_restore::Bool
    CheckpointableObjectGraph_CheckpointableObject_SerializedTensor(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct CheckpointableObjectGraph_CheckpointableObject_SerializedTensor

mutable struct CheckpointableObjectGraph_CheckpointableObject_SlotVariableReference <: ProtoType
    original_variable_node_id::Int32
    slot_name::AbstractString
    slot_variable_node_id::Int32
    CheckpointableObjectGraph_CheckpointableObject_SlotVariableReference(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct CheckpointableObjectGraph_CheckpointableObject_SlotVariableReference

mutable struct CheckpointableObjectGraph_CheckpointableObject <: ProtoType
    children::Base.Vector{CheckpointableObjectGraph_CheckpointableObject_ObjectReference}
    attributes::Base.Vector{CheckpointableObjectGraph_CheckpointableObject_SerializedTensor}
    slot_variables::Base.Vector{CheckpointableObjectGraph_CheckpointableObject_SlotVariableReference}
    CheckpointableObjectGraph_CheckpointableObject(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct CheckpointableObjectGraph_CheckpointableObject

mutable struct CheckpointableObjectGraph <: ProtoType
    nodes::Base.Vector{CheckpointableObjectGraph_CheckpointableObject}
    CheckpointableObjectGraph(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct CheckpointableObjectGraph

export CheckpointableObjectGraph_CheckpointableObject_ObjectReference, CheckpointableObjectGraph_CheckpointableObject_SerializedTensor, CheckpointableObjectGraph_CheckpointableObject_SlotVariableReference, CheckpointableObjectGraph_CheckpointableObject, CheckpointableObjectGraph

max_stars_repo_name: ajnavarro/language-dataset

max_issues_repo_name: ajnavarro/language-dataset

max_forks_repo_name: ajnavarro/language-dataset


text:
reorder(M::Matrix, V::Vector{Int}) = M[V,V]

function upper_sum(array::Matrix)
    S, (H,W) = 0, size(array)
    @assert W == H
    for i=1:H, j=i+1:W  S += array[i,j] end
    return S
end

function score_mat_cum_zeros(score_mat::Matrix, permutation::Vector,
                             to_shift = Array{Pair{Int,Int},1}())
    N = size(score_mat)[1]
    @assert size(score_mat)[2] == N
    @assert size(permutation)[1] == N
    @assert isperm(permutation) == true

    for u = 1:N
        perm_u = permutation[u]
        cumsum = score_mat[perm_u,perm_u]
        @assert cumsum == 0
        for v = u-1:-1:1
            cumsum += score_mat[perm_u,permutation[v]]
            cumsum == 0 && push!(to_shift, u=>v)
        end
        cumsum = score_mat[perm_u,perm_u]
        for v = u+1:N
            cumsum += score_mat[perm_u,permutation[v]]
            cumsum == 0 && push!(to_shift, u=>v)
        end
    end
    return to_shift
end

function rotate_part!(transformed::Vector, pair::Pair{Int,Int})
    (u, v) = pair; to_move = transformed[u]
    if u > v
        for k = u:-1:v+1 transformed[k] = transformed[k-1] end
    else
        for k = u:v-1 transformed[k] = transformed[k+1] end
    end
    transformed[v] = to_move
end

function slide_on_level(score_mat::Matrix, some_perm::Vector)
    N = size(score_mat)[1]
    @assert size(score_mat)[2] == N
    @assert size(some_perm)[1] == N
    @assert isperm(some_perm)
    visited, perms = 0, Vector{Vector{Int}}()

    @info "Collecting equivalent permutations..."
    push!(perms, some_perm)

    amount = length(perms)
    while visited < amount
        for i = visited+1:amount
            next_perm = perms[i]
            pairs = score_mat_cum_zeros(score_mat, next_perm)
            for pair in pairs
                permutation = copy(next_perm)
                rotate_part!(permutation, pair)
                in(permutation, perms) && continue
                push!(perms, permutation)
            end
        end
        visited = amount
        amount = length(perms)
        @debug "$amount\r"
        if amount > 20000
            @warn "More than 20k permutations!"
            break
        end
    end
    return sort(perms)
end

function is_refinable(array::Matrix)
    (N,M) = size(array)
    @assert N == M

    for i=1:N
        cumsum = 0
        for j=i-1:-1:1
            cumsum += array[i,j]
            if cumsum < 0
                array[i,j] == 0 && continue
                @info "$i=>$j : $cumsum"
            end
        end
        cumsum = 0
        for j=i+1:N
            cumsum += array[i,j]
            if cumsum > 0
                array[i,j] == 0 && continue
                @info "$i=>$j : $cumsum"
            end
        end
    end
    nothing
end

function dive_into(score_mat::Matrix, some_perm::Vector)
    (N,M) = size(score_mat)
    @assert N == M
    @assert isperm(some_perm)
    @assert N == size(some_perm)[1]
    steps = Dict{Pair{Int,Int},Int}()

    for i=1:N
        cumsum = 0
        for j=i-1:-1:1
            elem = reorder(score_mat,some_perm)[i,j]
            cumsum += elem
            if cumsum < 0
                elem == 0 && continue
                steps[i=>j] = cumsum
		@debug "$i=>$j : $cumsum"
            end
        end
        cumsum = 0
        for j=i+1:N
            elem = reorder(score_mat,some_perm)[i,j]
            cumsum += elem
            if cumsum > 0
                elem == 0 && continue
                steps[i=>j] = cumsum
		@debug "$i=>$j : $cumsum"
            end
        end
    end
    return steps
end

function find_key_of_maxvalue(divedict::Dict{Pair{Int,Int},Int})
    maxvalue = 0
    maxvaluekey = nothing
    for ((from,to),val) in divedict
        abs(val) < maxvalue && continue
        maxvalue = abs(val)
        maxvaluekey = from=>to
    end
    return maxvaluekey
end

function col_optimize(score_mat::Matrix)
    permutation = sortperm(col_negs_sums(score_mat))
    @debug "$permutation : $(upper_sum(reorder(score_mat, permutation)))"
    step = 1
    while true
        divedict = dive_into(score_mat, permutation)
        way = find_key_of_maxvalue(divedict)
        way == nothing && break
        rotate_part!(permutation, way)
        @debug "($step) $way -> $permutation : $(upper_sum(reorder(score_mat, permutation)))"
        step += 1
    end
    @info "col_optimize() stopped at $permutation."
    @info "Upper-triangle sum is $(upper_sum(reorder(score_mat, permutation)))."
    return slide_on_level(score_mat, permutation)
end

function row_optimize(score_mat::Matrix)
    permutation = reverse(sortperm(row_negs_sums(score_mat)))
    @debug "$permutation : $(upper_sum(reorder(score_mat, permutation)))"
    step = 1
    while true
        divedict = dive_into(score_mat, permutation)
        way = find_key_of_maxvalue(divedict)
        way == nothing && break
        rotate_part!(permutation, way)
        @debug "($step) $way -> $permutation : $(upper_sum(reorder(score_mat, permutation)))"
        step += 1
    end
    @info "row_optimize() stopped at $permutation."
    @info "Upper-triangle sum is $(upper_sum(reorder(score_mat, permutation)))."
    return slide_on_level(score_mat, permutation)
end

function consensus(score_mat::Matrix)
    @assert score_mat' == -score_mat
    solution = Vector{Vector{Int}}()
    @info "Sum of negative values is $(sum(score_mat[score_mat.<0]))."
    if is_straight(score_mat)
        @info "Unique permutation"
        push!(solution, sortperm(col_negs(score_mat)))
    else
        if is_crossing(score_mat)
            @warn "Unavoidable excess!"
            @show Δ⁺(score_mat)
        else
            @show Δ⁺(score_mat)
            @assert Δ⁺(score_mat) == 0
        end
        col_pack = col_optimize(score_mat)
        col_reached = upper_sum(reorder(score_mat, col_pack[1]))
        row_pack = row_optimize(score_mat)
        row_reached = upper_sum(reorder(score_mat, row_pack[1]))
        if col_reached > row_reached
            @info "Upper-triangle sum reached: $row_reached"
            solution = row_pack
        else
            @info "Upper-triangle sum reached: $col_reached"
            solution = col_pack
        end
    end
    return solution
    nothing
end


text:
@testset "id3" begin
    attributes = filter(a -> a != "class", header)
    tree = id3(cars, 7, attributes)
    println(tree)
end


text:
module mmodal1
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexStructures.FEMMCorotBeamModule
using FinEtoolsFlexStructures.FEMMCorotBeamModule: FEMMCorotBeam
stiffness = FEMMCorotBeamModule.stiffness
mass = FEMMCorotBeamModule.mass
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using Test
function test()
        # Parameters:
        E=71240.0;#MPa
        nu=0.31;# Poisson ratio
        rho=5e-9;
        b=3.0; h=30.0; L=240.0; # cross-sectional dimensions and length of each leg in millimeters
        # Choose the mass formulation:
        mass_type=2;
        scale = 0.4

        # Reference frequencies
        reffs = [11.2732, 30.5269]
        neigvs = 2;

        # Cross-sectional properties
        cs = CrossSectionRectangle(s -> b, s -> h, s -> [0.0, 1.0, 0.0])

        # Select the number of elements per leg.
        xyz =
        n = 4;
        members = []
        push!(members, frame_member([0 0 L; L 0 L], n, cs))
        push!(members, frame_member([L 0 L; L 0 0], n, cs))
        fens, fes = merge_members(members; tolerance = L / 10000);

        # Material properties
        material = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

        # Construct the requisite fields, geometry and displacement
        # Initialize configuration variables
        geom0 = NodalField(fens.xyz)
        u0 = NodalField(zeros(size(fens.xyz,1), 3))
        Rfield0 = initial_Rfield(fens)
        dchi = NodalField(zeros(size(fens.xyz,1), 6))

        # Apply EBC's
        l1 = selectnode(fens; box = [0 0 0 0 L L], tolerance = L/10000)
        for i in [1,2,3,4,5,6]
            setebc!(dchi, l1, true, i)
        end
        applyebc!(dchi)
        numberdofs!(dchi);

        # Assemble the global discrete system
        femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material)
        K = stiffness(femm, geom0, u0, Rfield0, dchi);
        M = mass(femm, geom0, u0, Rfield0, dchi);

        # Solve the eigenvalue problem
        d,v,nev,nconv = eigs(K, M; nev=2*neigvs, which=:SM)
        fs = real(sqrt.(complex(d)))/(2*pi)
        # println("Natural frequencies: $fs [Hz]")
        # println("Reference: $reffs [Hz]")

        @test norm(reffs - fs[1:length(reffs)]) ./ norm(reffs) <= 3.0e-5

        # Visualize vibration modes
        # mode = 2
        # tbox = plot_space_box([[-L/2 -L/2 0]; [L/2 L/2 1.1*L]])
        # tenv0 = plot_solid(fens, fes; x = geom0.values, u = 0.0.*dchi.values[:, 1:3], R = Rfield0.values, facecolor = "rgb(125, 155, 125)", opacity = 0.3);
        # plots = cat(tbox, tenv0; dims = 1)
        # pl = render(plots; title = "Mode $(mode)")
        # Rfield1 = deepcopy(Rfield0)
        # for xscale in scale.*sin.(collect(0:1:72).*(2*pi/21))
        #     scattersysvec!(dchi, xscale.*v[:, mode])
        #     Rfield1 = deepcopy(Rfield0)
        #     update_rotation_field!(Rfield1, dchi)
        #     tenv1 = plot_solid(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield1.values, facecolor = "rgb(25, 255, 25)");
        #     plots = cat(tbox, tenv0, tenv1; dims = 1)
        #     react!(pl, plots, pl.plot.layout)
        #     sleep(0.115)
        # end

        return true
    end # argyr_l_frame_modal_anim

end
using .mmodal1
mmodal1.test()

module mmassform11
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexStructures.FEMMCorotBeamModule
using FinEtoolsFlexStructures.FEMMCorotBeamModule: FEMMCorotBeam
stiffness = FEMMCorotBeamModule.stiffness
mass = FEMMCorotBeamModule.mass
geostiffness = FEMMCorotBeamModule.geostiffness
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using FinEtoolsFlexStructures.FESetCorotBeamModule: FESetL2CorotBeam, local_frame_and_def!, local_mass_original!, local_stiffness!, natural_forces!, local_geometric_stiffness!, local_forces!
using FinEtoolsFlexStructures.FESetCorotBeamModule: local_mass_CONSISTENT_WITH_ROTATION_INERTIA!, MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA

using Test

function test()
    # Parameters:
    E = 1000000.0;
    nu = 0.3;
    rho = 2700
    L1 = L0 = L =   30.0; # Length of the beam
    b =  0.5; # width
    h =  4.0; # height
    ecoords0 = fill(0.0, 2, 3);
    ecoords1 = fill(0.0, 2, 3)
    edisp1 = fill(0.0, 2, 3);
    evel1 = fill(0.0, 2, 6);
    evel1f = fill(0.0, 2, 6)
    dofnums = zeros(FInt, 1, 12);
    F0 = fill(0.0, 3, 3);
    Ft = fill(0.0, 3, 3);
    FtI = fill(0.0, 3, 3);
    FtJ = fill(0.0, 3, 3)
    Te = fill(0.0, 12, 12)
    tempelmat1 = fill(0.0, 12, 12);
    tempelmat2 = fill(0.0, 12, 12);
    tempelmat3 = fill(0.0, 12, 12)
    elmat = fill(0.0, 12, 12);
    elmatTe = fill(0.0, 12, 12);
    elmato = fill(0.0, 12, 12)
    elvec = fill(0.0, 12);
    elvecf = fill(0.0, 12)
    aN = fill(0.0, 6, 12)
    dN = fill(0.0, 6)
    DN = fill(0.0, 6, 6)
    PN = fill(0.0, 6)
    LF = fill(0.0, 12)
    R1I = fill(0.0, 3, 3);
    R1J = fill(0.0, 3, 3);
    OS = fill(0.0, 3, 3)

    # Reference frequencies
    reffs = [48.5475, 124.1839]

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> b, s -> h, s -> [-1.0, 0.0, 0.0])

    # Select the number of elements per leg.
    n = 32;
    members = []
    push!(members, frame_member([0 0 0; 0 0 L], n, cs))
    fens, fes = merge_members(members; tolerance = L / 10000);

    A, I1, I2, I3, J, x1x2_vector = fes.A, fes.I1, fes.I2, fes.I3, fes.J, fes.x1x2_vector
    G = E / 2 / (1 + nu)
    i = 1

    A, I1, I2, I3, rho, L = 1.3, 35.1, 32.3, 16.5, 10000.0, 3.333
    MM = fill(0.0, 12, 12)
    local_mass_CONSISTENT_WITH_ROTATION_INERTIA!(MM, A, I1, I2, I3, rho, L)
    MMp = fill(0.0, 12, 12)
    local_mass_original!(MMp, A, I1, I2, I3, rho, L, MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA)
    @test norm(MM - MMp) / norm(MM) <= 1.0e-6
    true
end
end
using .mmassform11
mmassform11.test()

module mmassform10
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexStructures.FEMMCorotBeamModule
using FinEtoolsFlexStructures.FEMMCorotBeamModule: FEMMCorotBeam
stiffness = FEMMCorotBeamModule.stiffness
mass = FEMMCorotBeamModule.mass
geostiffness = FEMMCorotBeamModule.geostiffness
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using FinEtoolsFlexStructures.FESetCorotBeamModule: FESetL2CorotBeam, local_frame_and_def!, local_mass_original!, local_stiffness!, natural_forces!, local_geometric_stiffness!, local_forces!
using FinEtoolsFlexStructures.FESetCorotBeamModule: local_mass_CONSISTENT_NO_ROTATION_INERTIA!, MASS_TYPE_CONSISTENT_NO_ROTATION_INERTIA

using Test
function test()
    # Parameters:
    E = 1000000.0;
    nu = 0.3;
    rho = 2700
    L1 = L0 = L =   30.0; # Length of the beam
    b =  0.5; # width
    h =  4.0; # height
    ecoords0 = fill(0.0, 2, 3);
    ecoords1 = fill(0.0, 2, 3)
    edisp1 = fill(0.0, 2, 3);
    evel1 = fill(0.0, 2, 6);
    evel1f = fill(0.0, 2, 6)
    dofnums = zeros(FInt, 1, 12);
    F0 = fill(0.0, 3, 3);
    Ft = fill(0.0, 3, 3);
    FtI = fill(0.0, 3, 3);
    FtJ = fill(0.0, 3, 3)
    Te = fill(0.0, 12, 12)
    tempelmat1 = fill(0.0, 12, 12);
    tempelmat2 = fill(0.0, 12, 12);
    tempelmat3 = fill(0.0, 12, 12)
    elmat = fill(0.0, 12, 12);
    elmatTe = fill(0.0, 12, 12);
    elmato = fill(0.0, 12, 12)
    elvec = fill(0.0, 12);
    elvecf = fill(0.0, 12)
    aN = fill(0.0, 6, 12)
    dN = fill(0.0, 6)
    DN = fill(0.0, 6, 6)
    PN = fill(0.0, 6)
    LF = fill(0.0, 12)
    R1I = fill(0.0, 3, 3);
    R1J = fill(0.0, 3, 3);
    OS = fill(0.0, 3, 3)

    # Reference frequencies
    reffs = [48.5475, 124.1839]

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> b, s -> h, s -> [-1.0, 0.0, 0.0])

    # Select the number of elements per leg.
    n = 32;
    members = []
    push!(members, frame_member([0 0 0; 0 0 L], n, cs))
    fens, fes = merge_members(members; tolerance = L / 10000);

    A, I1, I2, I3, J, x1x2_vector = fes.A, fes.I1, fes.I2, fes.I3, fes.J, fes.x1x2_vector
    G = E / 2 / (1 + nu)
    i = 1

    A, I1, I2, I3, rho, L = 1.3, 35.1, 32.3, 16.5, 10000.0, 3.333
    MM = fill(0.0, 12, 12)
    local_mass_CONSISTENT_NO_ROTATION_INERTIA!(MM, A, I1, I2, I3, rho, L)
    MMp = fill(0.0, 12, 12)
    local_mass_original!(MMp, A, I1, I2, I3, rho, L, MASS_TYPE_CONSISTENT_NO_ROTATION_INERTIA)
    @test norm(MM - MMp) / norm(MM) <= 1.0e-6
    true
end
end
using .mmassform10
mmassform10.test()

module mmassform12
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexStructures.FEMMCorotBeamModule
using FinEtoolsFlexStructures.FEMMCorotBeamModule: FEMMCorotBeam
stiffness = FEMMCorotBeamModule.stiffness
mass = FEMMCorotBeamModule.mass
geostiffness = FEMMCorotBeamModule.geostiffness
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using FinEtoolsFlexStructures.FESetCorotBeamModule: FESetL2CorotBeam, local_frame_and_def!, local_mass_original!, local_stiffness!, natural_forces!, local_geometric_stiffness!, local_forces!
using FinEtoolsFlexStructures.FESetCorotBeamModule: local_mass_LUMPED_DIAGONAL_NO_ROTATION_INERTIA!, MASS_TYPE_LUMPED_DIAGONAL_NO_ROTATION_INERTIA

using Test
function test()
    # Parameters:
    E = 1000000.0;
    nu = 0.3;
    rho = 2700
    L1 = L0 = L =   30.0; # Length of the beam
    b =  0.5; # width
    h =  4.0; # height
    ecoords0 = fill(0.0, 2, 3);
    ecoords1 = fill(0.0, 2, 3)
    edisp1 = fill(0.0, 2, 3);
    evel1 = fill(0.0, 2, 6);
    evel1f = fill(0.0, 2, 6)
    dofnums = zeros(FInt, 1, 12);
    F0 = fill(0.0, 3, 3);
    Ft = fill(0.0, 3, 3);
    FtI = fill(0.0, 3, 3);
    FtJ = fill(0.0, 3, 3)
    Te = fill(0.0, 12, 12)
    tempelmat1 = fill(0.0, 12, 12);
    tempelmat2 = fill(0.0, 12, 12);
    tempelmat3 = fill(0.0, 12, 12)
    elmat = fill(0.0, 12, 12);
    elmatTe = fill(0.0, 12, 12);
    elmato = fill(0.0, 12, 12)
    elvec = fill(0.0, 12);
    elvecf = fill(0.0, 12)
    aN = fill(0.0, 6, 12)
    dN = fill(0.0, 6)
    DN = fill(0.0, 6, 6)
    PN = fill(0.0, 6)
    LF = fill(0.0, 12)
    R1I = fill(0.0, 3, 3);
    R1J = fill(0.0, 3, 3);
    OS = fill(0.0, 3, 3)

    # Reference frequencies
    reffs = [48.5475, 124.1839]

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> b, s -> h, s -> [-1.0, 0.0, 0.0])

    # Select the number of elements per leg.
    n = 32;
    members = []
    push!(members, frame_member([0 0 0; 0 0 L], n, cs))
    fens, fes = merge_members(members; tolerance = L / 10000);

    A, I1, I2, I3, J, x1x2_vector = fes.A, fes.I1, fes.I2, fes.I3, fes.J, fes.x1x2_vector
    G = E / 2 / (1 + nu)
    i = 1

    A, I1, I2, I3, rho, L = 1.3, 35.1, 32.3, 16.5, 10000.0, 3.333
    MM = fill(0.0, 12, 12)
    local_mass_LUMPED_DIAGONAL_NO_ROTATION_INERTIA!(MM, A, I1, I2, I3, rho, L)
    MMp = fill(0.0, 12, 12)
    local_mass_original!(MMp, A, I1, I2, I3, rho, L, MASS_TYPE_LUMPED_DIAGONAL_NO_ROTATION_INERTIA)
    @test norm(MM - MMp) / norm(MM) <= 1.0e-6
    true
end
end
using .mmassform12
mmassform12.test()

module mmassform13
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.CrossSectionModule: CrossSectionRectangle
using FinEtoolsFlexStructures.MeshFrameMemberModule: frame_member, merge_members
using FinEtoolsFlexStructures.FEMMCorotBeamModule
using FinEtoolsFlexStructures.FEMMCorotBeamModule: FEMMCorotBeam
stiffness = FEMMCorotBeamModule.stiffness
mass = FEMMCorotBeamModule.mass
geostiffness = FEMMCorotBeamModule.geostiffness
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using FinEtoolsFlexStructures.FESetCorotBeamModule: FESetL2CorotBeam, local_frame_and_def!, local_mass_original!, local_stiffness!, natural_forces!, local_geometric_stiffness!, local_forces!
using FinEtoolsFlexStructures.FESetCorotBeamModule: local_mass_LUMPED_DIAGONAL_WITH_ROTATION_INERTIA!, MASS_TYPE_LUMPED_DIAGONAL_WITH_ROTATION_INERTIA

using Test
function test()
    # Parameters:
    E = 1000000.0;
    nu = 0.3;
    rho = 2700
    L1 = L0 = L =   30.0; # Length of the beam
    b =  0.5; # width
    h =  4.0; # height
    ecoords0 = fill(0.0, 2, 3);
    ecoords1 = fill(0.0, 2, 3)
    edisp1 = fill(0.0, 2, 3);
    evel1 = fill(0.0, 2, 6);
    evel1f = fill(0.0, 2, 6)
    dofnums = zeros(FInt, 1, 12);
    F0 = fill(0.0, 3, 3);
    Ft = fill(0.0, 3, 3);
    FtI = fill(0.0, 3, 3);
    FtJ = fill(0.0, 3, 3)
    Te = fill(0.0, 12, 12)
    tempelmat1 = fill(0.0, 12, 12);
    tempelmat2 = fill(0.0, 12, 12);
    tempelmat3 = fill(0.0, 12, 12)
    elmat = fill(0.0, 12, 12);
    elmatTe = fill(0.0, 12, 12);
    elmato = fill(0.0, 12, 12)
    elvec = fill(0.0, 12);
    elvecf = fill(0.0, 12)
    aN = fill(0.0, 6, 12)
    dN = fill(0.0, 6)
    DN = fill(0.0, 6, 6)
    PN = fill(0.0, 6)
    LF = fill(0.0, 12)
    R1I = fill(0.0, 3, 3);
    R1J = fill(0.0, 3, 3);
    OS = fill(0.0, 3, 3)

    # Reference frequencies
    reffs = [48.5475, 124.1839]

    # Cross-sectional properties
    cs = CrossSectionRectangle(s -> b, s -> h, s -> [-1.0, 0.0, 0.0])

    # Select the number of elements per leg.
    n = 32;
    members = []
    push!(members, frame_member([0 0 0; 0 0 L], n, cs))
    fens, fes = merge_members(members; tolerance = L / 10000);

    A, I1, I2, I3, J, x1x2_vector = fes.A, fes.I1, fes.I2, fes.I3, fes.J, fes.x1x2_vector
    G = E / 2 / (1 + nu)
    i = 1

    A, I1, I2, I3, rho, L = 1.3, 35.1, 32.3, 16.5, 10000.0, 3.333
    MM = fill(0.0, 12, 12)
    local_mass_LUMPED_DIAGONAL_WITH_ROTATION_INERTIA!(MM, A, I1, I2, I3, rho, L)
    MMp = fill(0.0, 12, 12)
    local_mass_original!(MMp, A, I1, I2, I3, rho, L, MASS_TYPE_LUMPED_DIAGONAL_WITH_ROTATION_INERTIA)
    @test norm(MM - MMp) / norm(MM) <= 1.0e-6
    true
end
end
using .mmassform13
mmassform13.test()


text:
using NOVAS
import NOVAS.LibNOVAS

function fund_args(t)
    a = zeros(Cdouble, 5)
    LibNOVAS.fund_args(t, a)
    return a
end

function norm_ang(angle)
    return LibNOVAS.norm_ang(angle)
end

function ee_ct(jd_high, jd_low; accuracy)
    if accuracy == :full
        acc = 0
    elseif accuracy == :reduced
        acc = 1
    end
    return LibNOVAS.ee_ct(jd_high, jd_low, acc)
end

function iau2000a(jd_high, jd_low)
    dpsi = Ref{Cdouble}(0.0)
    deps = Ref{Cdouble}(0.0)
    LibNOVAS.iau2000a(jd_high, jd_low, dpsi, deps)
    return dpsi[], deps[]
end

function nu2000k(jd_high, jd_low)
    dpsi = Ref{Cdouble}(0.0)
    deps = Ref{Cdouble}(0.0)
    LibNOVAS.nu2000k(jd_high, jd_low, dpsi, deps)
    return dpsi[], deps[]
end

function nutation_angles(t; accuracy)
    if accuracy == :full
        acc = 0
    elseif accuracy == :reduced
        acc = 1
    end
    dpsi = Ref{Cdouble}(0.0)
    deps = Ref{Cdouble}(0.0)
    LibNOVAS.nutation_angles(t, acc, dpsi, deps)
    return dpsi[], deps[]
end

function e_tilt(jd_tdb; accuracy)
    if accuracy == :full
        acc = 0
    elseif accuracy == :reduced
        acc = 1
    end
    mobl = Ref{Cdouble}(0.0)
    tobl = Ref{Cdouble}(0.0)
    ee = Ref{Cdouble}(0.0)
    dpsi = Ref{Cdouble}(0.0)
    deps = Ref{Cdouble}(0.0)
    LibNOVAS.e_tilt(jd_tdb, acc, mobl, tobl, ee, dpsi, deps)
    return mobl[], tobl[], ee[], dpsi[], deps[]
end

function mean_obliq(jd_tdb)
    return LibNOVAS.mean_obliq(jd_tdb)
end

function frame_tie(pos, direction)
    if direction == :dynamic2icrs
        dir = -1
    else
        dir = 0
    end
    pos2 = zeros(Cdouble, 3)
    LibNOVAS.frame_tie(pos, dir, pos2)
    return pos2
end

function nutation(jd_tdb, pos; direction, accuracy)
    if accuracy == :full
        acc = 0
    elseif accuracy == :reduced
        acc = 1
    end

    if direction == :mean2true
        dir = 0
    elseif direction == :true2mean
        dir = 1
    end

    pos2 = zeros(Cdouble, 3)
    LibNOVAS.nutation(jd_tdb, dir, acc, pos, pos2)
    return pos2
end

function precession(jd_tdb1, pos, jd_tdb2)
    pos2 = zeros(Cdouble, 3)
    LibNOVAS.precession(jd_tdb1, pos, jd_tdb2, pos2)
    return pos2
end

function ira_equinox(jd_tdb; equinox, accuracy)
    if accuracy == :full
        acc = 0
    elseif accuracy == :reduced
        acc = 1
    end
    if equinox == :mean
        equ = 0
    elseif equinox == :true
        equ = 1
    end
    return LibNOVAS.ira_equinox(jd_tdb, equ, acc)
end

function cio_location(jd_tdb; accuracy)
    if accuracy == :full
        acc = 0
    elseif accuracy == :reduced
        acc = 1
    end
    ra_cio = Ref{Cdouble}(0.0)
    ref_sys = Ref{Cshort}(0)
    LibNOVAS.cio_location(jd_tdb, acc, ra_cio, ref_sys)
    return ra_cio[]
end

function cio_basis(jd_tdb, ra_cio; accuracy)
    if accuracy == :full
        acc = 0
    elseif accuracy == :reduced
        acc = 1
    end
    x = zeros(Cdouble, 3)
    y = zeros(Cdouble, 3)
    z = zeros(Cdouble, 3)
    LibNOVAS.cio_basis(jd_tdb, ra_cio, 2, acc, x, y, z)
    return x, y, z
end

function era(jd_high, jd_low)
    return LibNOVAS.era(jd_high, jd_low)
end

function tdb2tt(jdb_jd)
    tt_jd = Ref{Cdouble}(0.0)
    secdiff = Ref{Cdouble}(0.0)
    LibNOVAS.tdb2tt(jdb_jd, tt_jd, secdiff)
    return tt_jd[], secdiff[]
end

function sidereal_time(jd_high, jd_low, delta_t; gst_type, method, accuracy)
    if accuracy == :full
        acc = 0
    elseif accuracy == :reduced
        acc = 1
    end

    if gst_type == :mean
        gstt = 0
    elseif gst_type == :apparent
        gstt = 1
    end

    if method == :CIO
        meth = 0
    elseif method == :equinox
        meth = 1
    end

    gst = Ref{Cdouble}(0.0)
    LibNOVAS.sidereal_time(jd_high, jd_low, delta_t, gstt, meth, acc, gst)
    return gst[]
end

function wobble(tjd, xp, yp, pos; direction)
    if direction == :itrs2terr
        dir = 0
    elseif direction == :terr2itrs
        dir = 1
    end
    pos2 = zeros(Cdouble, 3)
    LibNOVAS.wobble(tjd, dir, xp, yp, pos, pos2)
    return pos2
end

function spin(angle, pos)
    pos2 = zeros(Cdouble, 3)
    LibNOVAS.spin(angle, pos, pos2)
    return pos2
end

function ter2cel(jd_ut_high, jd_ut_low, delta_t, pos; method, accuracy, option, xp, yp)
    if accuracy == :full
        acc = 0
    elseif accuracy == :reduced
        acc = 1
    end
    if method == :CIO
        meth = 0
    elseif method == :equinox
        meth = 1
    end
    if option == :GCRS
        opt = 0
    elseif option == :equinox
        opt = 1
    end
    vec2 = zeros(Cdouble, 3)
    LibNOVAS.ter2cel(jd_ut_high, jd_ut_low, delta_t, meth, acc, opt, xp, yp, pos, vec2)
    return vec2
end

function refract(location, zd_obs; ref_option)
    if ref_option == :standard
        ro = 1
    elseif ref_option == :location
        ro = 2
    end
    loc = convert(NOVAS.LibNOVAS.on_surface, location)
    return LibNOVAS.refract(Ref{NOVAS.LibNOVAS.on_surface}(loc), ro, zd_obs)
end

function equ2hor(jd_ut1, delta_t, ra, dec, location; accuracy, xp, yp, ref_option)
    zd = Ref{Cdouble}(0.0)
    az = Ref{Cdouble}(0.0)
    rar = Ref{Cdouble}(0.0)
    decr = Ref{Cdouble}(0.0)
    if accuracy == :full
        acc = 0
    elseif accuracy == :reduced
        acc = 1
    end
    if ref_option == :none
        ro = 0
    elseif ref_option == :standard
        ro = 1
    elseif ref_option == :location
        ro = 2
    end
    loc = Ref{NOVAS.LibNOVAS.on_surface}(convert(NOVAS.LibNOVAS.on_surface, location))
    LibNOVAS.equ2hor(jd_ut1, delta_t, acc, xp, yp, loc, ra, dec, ro, zd, az, rar, decr)
    return zd[], az[], rar[], decr[]
end

import Base.convert

function Base.convert(::Type{NOVAS.LibNOVAS.on_surface}, obj::OnSurface{Float64})
    # There's gotta be a better way to do this
    return NOVAS.LibNOVAS.on_surface(obj.latitude, obj.longitude, obj.height,
                                     obj.temperature, obj.pressure)
end


text:
function get_sunits_from_subfield_data(OK, N; recursive::Bool = false, compact::Bool = true)
  K = nf(OK)
  fbbound = Hecke.factor_base_bound_grh(OK)
  S = prime_ideals_up_to(OK, fbbound)
  @vprint :NormRelation 1 "Factor base bound: $fbbound\n"
  c, UZK = _setup_for_norm_relation_fun(K, S)

  docompact = false
  onlyp = 0

  if Hecke.isprime_power(index(N)) && compact
    docompact = true
    _, onlyp = ispower(index(N))
  end

  @vprint :NormRelation 1 "Doing something in the subfields\n"
  if !docompact
    _add_sunits_from_brauer_relation!(c, UZK, N)
  else
    @vprint :NormRelation 1 "Using the compact presentation\n"
    _add_sunits_from_brauer_relation!(c, UZK, N, compact = onlyp)
  end
  return c, UZK
end

function class_group_via_brauer(O::NfOrd, N::NormRelation; compact::Bool = true)
  K = N.K
  OK = lll(maximal_order(nf(O)))

  c, UZK = get_sunits_from_subfield_data(OK, N, compact = compact)
  for (p, e) in factor(index(N))
    @vprint :NormRelation 1 "Saturating at $p \n"
    b = Hecke.saturate!(c, UZK, Int(p), 3.5, easy_root = true)
    while b
      idx = Hecke._validate_class_unit_group(c, UZK)[1]
      @vprint :NormRelation 1 "Index bound from analysis $idx\n"
      b = Hecke.saturate!(c, UZK, Int(p), 3.5, easy_root = true)
    end
  end
  idx = Hecke._validate_class_unit_group(c, UZK)[1]
  @vprint :NormRelation 1 "Index is $idx (should be 1)!\n"
  if idx != 1
    @vprint :NormRelation 1 "Index is $idx (should be 1)!\n"
  end
  @assert idx == 1
  @vprint :NormRelation 1 "\n"
  c.finished = true
  UZK.finished = true
  Hecke._set_ClassGrpCtx_of_order(OK, c)
  Hecke._set_UnitGrpCtx_of_order(OK, UZK)
  return class_group(c, O)
end


text:
using LinearAlgebra
using Random, Distributions
using Plots
Plots.scalefontsizes(1.5)

# Global Terms
Random.seed!(1234)

# Logistic Regression

function f(β)
    #=
    Objective Function
    x: n x d
    y: d x 1
    β: d x 1
    λ: penalty
    =#
    result =  sum((X * β - y) .^ 2) / 2 / n + λ *norm(β,1)
    return result
end

function ∇f(β)
    #=
    This is already noisy gradient.
    β: 100 x 1
    σ: noise
    =#
    return X' * (X * β - y) ./n .+ λ .* sign.(β) .+ rand(Normal(0,σ^2), (d,1))
end

# MASG

function MASG(f, ∇f, nk, α, μ)
    #=
        f:  objective functin
        ∇f: gradient of f
        X:  design matrix
        y:  {-1,1} class
        a_n:an array containing each step each iteration
        K:  total number of iterations
        α:  an array containing step sizes
        μ:  μ-strongly convex
        λ:  penalty
    =#
    x_0 = x_1 = zeros(d)

    # Record Histories
    f_history = zeros(sum(nk))
    γ = 1

    for k in 1:K
        x_0 = x_1
        for m in 1:nk[k]
            w = (1 - sqrt(μ * α[k])) / (1 + sqrt(μ * α[k]))
            ymk = (1 + w) .* x_1 - w .* x_0
            x_0 = x_1
            x_1 = ymk .- α[k] .* ∇f(ymk)

            # Record History
            f_history[γ] = f(x_1)
            γ = γ + 1
        end
    end
    return f_history
end

function compute_nk(κ, p, Δ, label)
    C = floor(sqrt(κ) * log(2^(p+2))) + 1
    grid = [1:K;]
    nk = 2 .^ grid
    if label == 3.7
        nk[1] = floor(sqrt(κ) * log(κ)) + 1
    elseif label == 4
        nk[1] = floor(sqrt(κ) * log(2 * L * Δ / (σ^2 * sqrt(κ)))) + 1
    end
    return nk
end

function compute_MASG_α(L)
    grid = [1:K;]
    α = 1 ./ (2 .^ (2 .* [1:K;]) .* L)
    α[1] = 1 / L
    return α
end

# Vanilla Gradient Descent
function GD(f, ∇f)
    x = zeros(d)
    f_history = zeros(η)
    for i in 1:η
        x = x - 1 / L * ∇f(x)
        f_history[i] = f(x)
    end
    return f_history
end

# Accelerated Gradient Descent
function AGD(f, ∇f)
    x_0 = x_1 = zeros(d)
    f_history = zeros(η)
    for i in 1:η
        y = x_1 + (i-1)/(i+2) *(x_1-x_0)
        x_0 = x_1
        x_1 = y - 1/L * ∇f(y)

        f_history[i] = f(x_1)
    end
    return f_history
end

# AC-SA
## Bregman Distance
function D(β1, β2)
    #=
    Prox function, Bregman Distance
    =#
    return f(β1) - f(β2) - (∇f(β2)' * (β1 .- β2))[1]
end

function ORDA(τ)
    #=
    ξ is the turning param
    M is constant
    τ is quadratic constant
    =#
    M = 0
    c = (sqrt(τ) * (σ + M))/(2 * sqrt(D(βstart, zeros(d))) )
    θ(t) = 2/(t+2)
    ν(t) = 2/(t+1)
    γ(t) = c * (t+1)^(3/2)+τ*L
    G = 0
    zt = z0 = xt = x0 = zeros(d)
    f_history = zeros(η+1)
    for t in 0:η
        f_history[t+1] = f(xt)
        p1 = (1 - θ(t)) * (μ + θ(t)^2 * γ(t)) / (θ(t)^2 * γ(t) + (1-θ(t)^2) * μ) .* xt
        p2 = ((1 - θ(t)) * θ(t) * μ + θ(t)^3 * γ(t)) / (θ(t)^2 * γ(t) + (1-θ(t)^2) * μ) * zt
        yt = p1 + p2
        G = G .+ ∇f(yt) ./ ν(t)
        gt = θ(t) * ν(t) .* G

        zp1 = ∇f(x0) .- λ .- (gt .+ λ) ./ (θ(t) * ν(t) * γ(t+1))
        zt = inv(X' * X) * (zp1 .* n + X' * y)

        xp1 = ∇f(yt)  .- λ .- (∇f(yt) .+ λ) / (μ / (τ * θ(t)^2) + γ(t) / τ)
        xt = inv(X' * X) * (xp1 .* n + X' * y)
     end
     return f_history
end

function ACSA(ν)
    α(t) = 2/(t+1)
    γ(t) = 4 * 2 * L / (ν * t * (t+1))

    xt = xtag = zeros(d)
    f_history = zeros(η+1)
    f_history[1] = f(xt)

    for t in 1:η
        xp1 = (1 - α(t)) * (μ + γ(t)) / (γ(t) + (1-α(t)^2) * μ) * xtag
        xp2 = α(t)*((1-α(t)) * μ +γ(t))/(γ(t) + (1-α(t)^2) * μ) * xt
        xtmd = xp1 + xp2

        p1 = (∇f(xtmd) * (μ-1) * α(t) .- α(t) * λ .+ ((1-α(t))*μ + γ(t))*∇f(xt)) / (μ + γ(t))
        p2 = (p1 .- λ) .*n + X' * y
        xt = inv(X' * X) * p2

        xtag = α(t) * xt + (1- α(t)) * xtag
        f_history[t+1] = f(xt)
    end

    return f_history
end

## Data setup
global n = 10000   # number of datapoints
global d = 100    # number of dimensions
global λ = 0.01   # Penalty for beta term
global ρ = 0
global X = rand(Normal(0,1), (n,d))
β = rand(Normal(0, 1), (d, 1))
global βstart = β
global y = X * β  .+ rand(Normal(0,1),(n,1))

# Experiments
global K = 10     # number of iterations for MASG
global η = 2000   # Number of iterations for other methods
global L = maximum(eigvals(X' * X))
# global L = norm(X,2)^2/n
global μ = 0

global σ = sqrt(1e-1)
# MASG
α_MASG = compute_MASG_α(L)
nk_7 = compute_nk(κ, 1, norm(β,2),3.7)
f_MASG_7 = MASG(f, ∇f, nk_7, α_MASG, μ)

## MASG Thm 4
nk_4 = compute_nk(κ, 1, norm(β,2), 4.0)
f_MASG_4 = MASG(f, ∇f, nk_4, α_MASG, μ)

# Vanilla Gradient Descent
f_GD = GD(f, ∇f)

# Vanilla Accelerated Gradient Descent
f_AGD = AGD(f, ∇f)

# ORDA
f_ORDA = ORDA(L)

ν = norm(X,2) / n
f_ACSA = ACSA(ν)

plot(f_GD, label="GD",legend=:bottomleft,w=2,linestyle=:dash,color=:red, xscale = :log, yscale=:log)
plot!(f_AGD, label="AGD", w=2, linestyle=:dash,color=:purple)
plot!(f_ORDA, label = "ORDA", w=2, color=:blue)
plot!(f_ACSA, label="ACSA",w=2, color=:brown)
plot!(f_MASG_7, label="MASG",w=2, linestyle=:solid, color=:green)
plot!(f_MASG_4, label="MASG*",w=2, linestyle=:solid, color=:black)
xlabel!("Iteration")

savefig("./figures/non_sc_10000_4.png")

max_stars_repo_name: chenyangzhu/noisy-gradients

max_issues_repo_name: chenyangzhu/noisy-gradients

max_forks_repo_name: chenyangzhu/noisy-gradients


text:
@testset "Vax Ints" begin
    i2_vax = [ 0x0001, 0xFFFF, 0x0100, 0xFF00, 0x3039, 0xCFC7 ]
    i2_ieee = Array{Int16}([ 1, -1, 256, -256, 12345, -12345 ])

    i4_vax = [  0x00000001,
                0xFFFFFFFF,
                0x00000100,
                0xFFFFFF00,
                0x00010000,
                0xFFFF0000,
                0x01000000,
                0xFF000000,
                0x075BCD15,
                0xF8A432EB ]

    i4_ieee = Array{Int32}([          1,
                                     -1,
                                    256,
                                   -256,
                                  65536,
                                 -65536,
                               16777216,
                              -16777216,
                              123456789,
                             -123456789 ])

    @testset "VaxInt16" begin
        @testset "Conversion..." begin
            for (vax, ieee) in zip(i2_vax, i2_ieee)
                @test VaxInt16(vax) == VaxInt16(ieee)
                @test convert(Int16, VaxInt16(vax)) == ieee
            end
        end

        @testset "Promotion..." begin
            @test isa(one(Int8)*VaxInt16(1), Int16)
            @test isa(VaxInt16(1)*VaxInt16(1), Int16)
            for t in [Float16, Float32, Float64, Int16, Int32, Int64, Int128]
                @test isa(one(t)*VaxInt16(1), t)
            end
        end
    end

    @testset "VaxInt32" begin
        @testset "Conversion..." begin
            for (vax, ieee) in zip(i4_vax, i4_ieee)
                @test VaxInt32(vax) == VaxInt32(ieee)
                @test convert(Int32, VaxInt32(vax)) == ieee
            end
        end

        @testset "Promotion..." begin
            for t in [subtypes(VaxInt); Int8; Int16]
                @test isa(one(t)*VaxInt32(1), Int32)
            end

            for t in [Float16, Float32, Float64, Int32, Int64, Int128]
                @test isa(one(t)*VaxInt32(1), t)
            end
        end
    end
end


text:
module OhMyArtifacts

using Pkg
using Pkg.Artifacts
using Pkg.Types: parse_toml, write_env_usage
using SHA
using Dates
using Printf
using TOML
using Scratch
using Pidfile

const SCRATCH_DIR = Ref{String}()
const ARTIFACTS_DIR = Ref{String}()

const ARTIFACTS_TOML_VAR_SYM = Ref{Symbol}(:my_artifacts)

function __init__()
    global SCRATCH_DIR, ARTIFACTS_DIR
    # create scratch space
    ARTIFACTS_DIR[] = @get_scratch!("artifacts")
    SCRATCH_DIR[] = dirname(ARTIFACTS_DIR[])

    orphan_file = orphanages_toml_path()
    last_gc_time = if isfile(orphan_file)
        modified_time(orphan_file)
    else
        orphanages_toml()
        now()
    end
    if now() - last_gc_time >= Day(7)
        find_orphanages()
    end

    return
end

export my_artifacts_toml!, @my_artifacts_toml!, create_my_artifact, bind_my_artifact!,
    my_artifact_hash, my_artifact_path, my_artifact_exists,
    @my_artifact, download_my_artifact!

"""
    my_artifacts_toml!(pkg::Union{Module,Base.UUID,Nothing})

Return the path to (or creates) "Artifacts.toml" for the given `pkg`.

See also: [`@my_artifacts_toml!`](@ref)
"""
function my_artifacts_toml!(pkg::Union{Module,Base.UUID,Nothing})
    uuid = Scratch.find_uuid(pkg)
    path = get_scratch!(@__MODULE__, string(uuid), pkg)
    return touch(joinpath(path, "Artifacts.toml"))
end

## utilities

function get_scratch_dir()
    return mkpath(scratch_dir(string(Base.PkgId(@__MODULE__).uuid)))
end

function get_artifacts_dir()
    return mkpath(joinpath(get_scratch_dir(), "artifacts"))
end

function get_artifacts_toml_sym()
    global ARTIFACTS_TOML_VAR_SYM
    return ARTIFACTS_TOML_VAR_SYM[]
end

orphanages_toml_path() = joinpath(get_scratch_dir(), "my_artifact_orphanages.toml")

function orphanages_toml()
    path = orphanages_toml_path()
    mkpath(dirname(path))
    touch(path)
end

function usages_toml()
    path = joinpath(get_scratch_dir(), "my_artifact_usage.toml")
    mkpath(dirname(path))
    touch(path)
end

function modified_time(file)
    return Dates.unix2datetime(mtime(file)) + round(now() - now(Dates.UTC), Hour)
end

## macros

"""
    @my_artifacts_toml!()

Convenience macro that gets/creates a "Artifacts.toml" and parented to the package the calling module belongs to.

See also: [`my_artifacts_toml!`](@ref)
"""
macro my_artifacts_toml!()
    uuid = Base.PkgId(__module__).uuid
    return quote
        my_artifacts_toml!($(esc(uuid)))
    end
end

include("./macro.jl")

## SHA256 ##

struct SHA256
    bytes::NTuple{32, UInt8}
end

function SHA256(bytes::Vector{UInt8})
    length(bytes) == 32 ||
        throw(ArgumentError("wrong number of bytes for SHA256 hash: $(length(bytes))"))
    return SHA256(ntuple(i->bytes[i], Val(32)))
end
SHA256(s::AbstractString) = SHA256(hex2bytes(s))
Base.parse(::Type{SHA256}, s::AbstractString) = SHA256(s)
function Base.tryparse(::Type{SHA256}, s::AbstractString)
    try
        return parse(SHA256, s)
    catch e
        if isa(e, ArgumentError)
            return nothing
        end
        rethrow(e)
    end
end

Base.string(hash::SHA256) = bytes2hex(hash.bytes)
Base.print(io::IO, hash::SHA256) = bytes2hex(io, hash.bytes)
Base.show(io::IO, hash::SHA256) = print(io, "SHA256(\"", hash, "\")")

Base.isless(a::SHA256, b::SHA256) = isless(a.bytes, b.bytes)
Base.hash(a::SHA256, h::UInt) = hash((SHA256, a.bytes), h)
Base.:(==)(a::SHA256, b::SHA256) = a.bytes == b.bytes

## artifacts ##

function load_my_artifacts_toml(artifacts_toml::String)
    artifact_dict = mkpidlock(joinpath(dirname(artifacts_toml), "artifact_lock")) do
        parse_toml(artifacts_toml)
    end
    return artifact_dict
end

"""
    my_artifact_hash(name::AbstractString, artifacts_toml::String)

Return the hash found in `artifacts_toml` with given `name`, or `nothing` if not found.
"""
function my_artifact_hash(name::AbstractString, artifacts_toml::String)
    artifact_dict = load_my_artifacts_toml(artifacts_toml)
    if haskey(artifact_dict, name)
        return SHA256(artifact_dict[name]["sha256"])
    else
        return nothing
    end
end

"""
    my_artifact_path(hash::SHA256)

Given an artifact (identified by SHA256 content hash), return its installation path. If the artifact does not exist,
 returns the location it would be installed to.

See also: [`my_artifact_exists`](@ref)
"""
function my_artifact_path(hash::SHA256)
    artifacts_dir = get_artifacts_dir()
    return joinpath(artifacts_dir, string(hash))
end

"""
    my_artifact_exists(hash::SHA256)

Returns whether or not the given artifact (identified by its SHA256 content hash) exists on-disk.
"""
function my_artifact_exists(hash::SHA256)
    return isfile(my_artifact_path(hash))
end

"""
    create_my_artifact(f::Function)

Creates a new artifact by doing `path = f(working_dir)`, hashing the returned `path`, and moving it to
 the artifact store. Returns the identifying hash of this artifact.

`f(working_dir)` should return an absolute path to a single file at the top level of `working_dir`.
"""
function create_my_artifact(f::Function)
    artifacts_dir = get_artifacts_dir()
    temp_dir = mktempdir(artifacts_dir)

    try
        # create file in temp dir
        filepath = f(temp_dir)

        # make sure the returned file path is correct
        !isabspath(filepath) && error("returned file path $filepath is not a abspath")
        !startswith(filepath, temp_dir) && error("returned file path $filepath result outside of the given dir")

        file = basename(filepath)
        filepath = joinpath(temp_dir, file)
        !isfile(filepath) && error("returned file is not at the top level of given dir")

        # calculate the file hash
        artifact_hash = SHA256(Pkg.GitTools.blob_hash(SHA.SHA256_CTX, filepath))
        new_path = joinpath(artifacts_dir, string(artifact_hash))

        # skip if file already exist
        filelock = mkpidlock(joinpath(dirname(new_path), "$(string(artifact_hash)).lock"))
        try
            if !isfile(new_path)
                mv(filepath, new_path)
                fmode = filemode(new_path)
                # read-only
                chmod(new_path, fmode & (typemax(fmode) ⊻ 0o222))
            end
        finally
            close(filelock)
        end

        return artifact_hash
    finally
        rm(temp_dir; recursive=true, force=true)
    end
end

"""
    bind_my_artifact!(artifacts_toml::String, name::AbstractString, hash::SHA256; force::Bool = false)

Writes a mapping of `name` -> `hash` within the given "Artifacts.toml" file. If `force` is set to `true`,
 this will overwrite a pre-existant mapping, otherwise an error is raised.
"""
function bind_my_artifact!(artifacts_toml::String, name::AbstractString, hash::SHA256; force::Bool = false)
    artifact_lock = mkpidlock(joinpath(dirname(artifacts_toml), "artifact_lock"))
    try
        artifact_dict = parse_toml(artifacts_toml)

        if !force && haskey(artifact_dict, name)
            error("Mapping for $name within $(artifacts_toml) already exists!")
        end

        meta =  Dict{String,Any}(
            "sha256" => string(hash),
        )

        artifact_dict[name] = meta
        let artifact_dict = artifact_dict
            open(artifacts_toml, "w") do io
                TOML.print(io, artifact_dict, sorted=true)
            end
        end
    finally
        close(artifact_lock)
    end

    track_my_artifacts(artifacts_toml, name, hash)
    return
end

function track_my_artifacts(artifacts_toml::String, name::AbstractString, hash::SHA256)
    artifacts_dir = get_artifacts_dir()

    usage_file = usages_toml()
    mkpath(dirname(usage_file))

    artifact_path = joinpath(artifacts_dir, string(hash))

    usage_lock = mkpidlock(joinpath(dirname(usage_file), ".my_artifact_usage_lock"))
    try
        usage_toml = if isfile(usage_file)
            parse_toml(usage_file)
        else
            Dict{String, Any}()
        end

        if haskey(usage_toml, artifact_path)
            artifact_usage = usage_toml[artifact_path]
            if haskey(artifact_usage, artifacts_toml)
                artifact_usage[artifacts_toml][name] = now()
            else
                artifact_usage[artifacts_toml] = Dict(name => now())
            end
        else
            usage_dict = Dict{String, Any}(
                artifacts_toml => Dict{String, Any}(
                    name => now()
                )
            )
            usage_toml[artifact_path] = usage_dict
        end

        let usage_toml = usage_toml
            open(usage_file, "w") do io
                TOML.print(io, usage_toml, sorted=true)
            end
        end
    finally
        close(usage_lock)
    end
    return
end

"""
    download_my_artifact!([downloadf::Function = Base.download], url, name::AbstractString, artifacts_toml::String;
                         force_bind::Bool = false, downloadf_kwarg...)

Convenient function that do download-create-bind together and return the content hash.
 Download function `downloadf` should take two position arguments
 (i.e. `downloadf(url, dest; downloadf_kwarg...)`). if `force_bind` is `true`,
 it will overwrite the pre-existant binding.

See also: [create_my_artifact](@ref), [bind_my_artifact!](@ref)
"""
function download_my_artifact!(downloadf::Function, url, name::AbstractString, artifacts_toml::String;
                              force_bind::Bool = false, downloadf_kwarg...)
    lockname = bytes2hex(sha1(name))
    hash = mkpidlock(joinpath(dirname(artifacts_toml), "$(lockname).lock")) do
        create_my_artifact() do artifact_dir
            downloadf(url, tempname(artifact_dir; cleanup=false); downloadf_kwarg...)
        end
    end
    bind_my_artifact!(artifacts_toml, name, hash; force = force_bind)

    return hash
end
download_my_artifact!(
    url, name::AbstractString, artifacts_toml::String;
    force_bind::Bool = false, downloadf_kwarg...,
) =
    download_my_artifact!(Base.download, url, name, artifacts_toml; force_bind, downloadf_kwarg...)

"""
    unbind_my_artifact!(artifacts_toml::String, name::AbstractString)

Unbind the given `name` from the "Artifacts.toml" file. Silently fails if no such binding exists within the file.
"""
function unbind_my_artifact!(artifacts_toml::String, name::AbstractString)
    artifact_lock = mkpidlock(joinpath(dirname(artifacts_toml), "artifact_lock"))
    try
        artifact_dict = parse_toml(artifacts_toml)

        !haskey(artifact_dict, name) && return

        hash = artifact_dict[name]["sha256"]

        delete!(artifact_dict, name)

        let artifact_dict = artifact_dict
            open(artifacts_toml, "w") do io
                TOML.print(io, artifact_dict, sorted=true)
            end
        end
    finally
        close(artifact_lock)
    end
    return
end

function find_orphanages(; collect_delay::Period=Day(7))
    artifacts_dir = get_artifacts_dir()

    usage_file = usages_toml()
    orphanage_file = orphanages_toml()

    !isfile(usage_file) && return

    usage_lock = mkpidlock(joinpath(dirname(usage_file), ".my_artifact_usage_lock"))
    try
        usage_toml = if isfile(usage_file)
            parse_toml(usage_file)
        else
            Dict{String, Any}()
        end

        # find orphan artifacts
        orphanage = Pair{String, DateTime}[]

        # artifact without binding
        for artifact_path in readdir(get_artifacts_dir(), join=true)
            if !haskey(usage_toml, artifact_path)
                push!(orphanage, artifact_path=>modified_time(artifact_path))
            end
        end

        # artifact with binding
        curr_gc_time = now()
        for (artifact_path, usage_dict) in usage_toml
            # check if all bindings are removed

            for (artifacts_toml, entrys) in usage_dict
                if !isfile(artifacts_toml)
                    # the toml is already removed, so no such usage exist
                    delete!(usage_dict, artifacts_toml)
                else
                    # else we read the toml
                    artifact_dict = load_my_artifacts_toml(artifacts_toml)

                    # check that all usage entry is still exist, or remove them
                    for entry in keys(entrys)
                        if !haskey(artifact_dict, entry)
                            delete!(entrys, entry)
                        end
                    end

                    # if no entrys, remove this usage
                    isempty(entrys) && delete!(usage_dict, artifacts_toml)
                end
            end

            if isempty(usage_dict)
                # no usage of this artifact exist
                # mark it as orphan and record the current time
                push!(orphanage, artifact_path=>curr_gc_time)
                delete!(usage_toml, artifact_path)
            end
        end

        # update usage toml
        let usage_toml = usage_toml
            open(usage_file, "w") do io
                TOML.print(io, usage_toml, sorted=true)
            end
        end

        # update the orphanage file
        old_orphans = if isfile(orphanage_file)
            parse_toml(orphanage_file)
        else
            Dict{String, DateTime}()
        end
        # merge old and new orphan list
        # *notice*: if the entry already on the list, keep the old recorded time
        new_orphans = merge(Dict(orphanage), old_orphans)

        # mark orphanage for deletion
        gc_time = now()
        deletion_list = String[]
        for artifact_path in keys(new_orphans)
            last_gc_time = new_orphans[artifact_path]
            if gc_time - last_gc_time >= collect_delay
                push!(deletion_list, artifact_path)
                delete!(new_orphans, artifact_path)
            end
        end

        # update the resulted orphanage to file
        let new_orphans = new_orphans
            open(orphanage_file, "w") do io
                TOML.print(io, new_orphans, sorted=true)
            end
        end

        # delete my artifacts
        pretty_byte_str = (size) -> begin
            bytes, mb = Base.prettyprint_getunits(size, length(Base._mem_units), Int64(1024))
            return @sprintf("%.3f %s", bytes, Base._mem_units[mb])
        end

        function file_size(path)
            size = try
                lstat(path).size
            catch ex
                @error("Failed to calculate size of $path", exception=ex)
            end
            return size
        end

        function delete_path(path)
            path_size = file_size(path)
            try
                Base.Filesystem.prepare_for_deletion(path)
                Base.rm(path; recursive=true, force=true)
            catch e
                @warn("Failed to delete $path", exception=e)
                return 0
            end
            return path_size
        end

        space_freed = 0
        for artifact_path in deletion_list
            space_freed += delete_path(artifact_path)
        end

        ndel = length(deletion_list)
        if ndel > 0
            s = ndel == 1 ? "" : "s"
            @info "$ndel MyArtifact$(s) deleted ($(pretty_byte_str(space_freed)))"
        end

    finally
        close(usage_lock)
    end
    return
end

end


text:
a, b, c, d, w, x, y, z = map(RealVariable, [:a, :b, :c, :d, :w, :x, :y, :z])
ix = TypedVariable{Integer}(:x)

@test x == x
@test x != y
@test x != ix
@test 1x == x
@test 1x == 1.0x
@test -x == -1x
@test x+2x == 3x
@test x-x == 0
@test 1+x-x == 1
@test 2*(x-x) == 0
@test x+y-x == y
@test 2.2x+4.4y == 2*(1.1x+2.2y)

@test_throws MethodError (1+2x) * (3+4y) # no polynomials

@test conj(x) == x
@test conj(2x) == 2x

@test 2 * [x y] == [2x 2y]
@test 2.0 * [x 3y] == [2.0x 6.0y]
@test [x 3y] * 2.0 == [2.0x 6.0y]

A = rand(Int, 2, 2)
B = rand(Int, 2, 2)

@test A + A*x + B + B*x == (A+B) + (A+B)*x

le1 = y + 2.3 * (w + z) - 3.4 * (a + b) - c - d + 1.23
le2 = 1.23 - d - 3.4b - 3.4a + y + 2.3w + 2.3z - c

@test le1 == le2
@test differentiate(le1, x) == 0
@test differentiate(le1, w) == 2.3

function test_eq()
  a1 = LinExpr(-1.2, Dict(x=>-1., y=>2.))
  a2 = LinExpr(-1.2, Dict(x=>-1., y=>2.))
  b = LinExpr(-3.4, Dict(x=>-3., y=>4.))
  @test a1 == a2
  @test a1 != b
end

test_eq()

function test_sum()
  a = LinExpr(-1.2, Dict(x=>-1., y=>2.))
  b = LinExpr(-3.4, Dict(x=>-3., y=>4.))
  c = LinExpr(-3.4, Dict(x=>1., z=>4.))
  r_a_b = LinExpr(-4.6, Dict(x=>-4., y=>6.))
  r_a_c = LinExpr(-4.6, Dict(y=>2., z=>4.))
  @test a+b == r_a_b
  @test a+c == r_a_c
end

test_sum()

@test convert(RealVariable, LinExpr{Float64, RealVariable}(Dict(x=>1.0))) == x

@test_throws InexactError convert(RealVariable, LinExpr(0.0, Dict(x=>1.0, y=>1.0)))
@test_throws InexactError convert(RealVariable, LinExpr{Float64, RealVariable}(1.23))

@test convert(Float64, LinExpr{Float64, RealVariable}(1.23)) == 1.23
@test_throws InexactError convert(Float64, LinExpr(1.23, Dict(x=>4.56)))