text:
using LinearAlgebraicRepresentation, ViewerGL, SparseArrays
Lar = LinearAlgebraicRepresentation; GL = ViewerGL

n,m,p = 1,1,1
V,(VV,EV,FV,CV) = Lar.cuboidGrid([n,m,p],true)
cube = V,FV,EV

threecubes = Lar.Struct([ cube,
    Lar.t(.3,.4,.25), Lar.r(pi/5,0,0), Lar.r(0,0,pi/12), cube,
    Lar.t(-.2,.4,-.2), Lar.r(0,pi/5,0), Lar.r(0,pi/12,0), cube ])

V,FV,EV = Lar.struct2lar(threecubes)
GL.VIEW([ GL.GLGrid(V,FV), GL.GLFrame ])

cop_EV = convert(Lar.ChainOp, Lar.coboundary_0(EV::Lar.Cells));
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');

V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement( W, cop_EV, cop_FE)

#triangulated_faces = Lar.triangulate(V, [copEV, copFE])
W = convert(Lar.Points, V')
V,CVs,FVs,EVs = Lar.pols2tria(W, copEV, copFE, copCF)

GL.VIEW(GL.GLExplode(V,FVs,1.5,1.5,1.5,99,1));
GL.VIEW(GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1));
meshes = GL.GLExplode(V,CVs[1:end],8,4,6,99,0.5);
GL.VIEW( push!( meshes, GL.GLFrame) );


text:
# TODO: Use Preferences?
const DEBUG = true

if !isdefined(Base, Symbol("@var_str"))
    macro var_str(s)
        Symbol(s)
    end
end

macro noop(args...)
end

if DEBUG
    const var"@_debug" = var"@debug"
    const var"@_assert" = var"@assert"
else
    const var"@_debug" = var"@noop"
    const var"@_assert" = var"@noop"
end

const KEY_SENTINEL = Symbol("!!__TASKDAGRECORDERS_KEY_SENTINEL__!!")

"""
Like `get!(f, d, k)`, but safe to mutate `f` in `d` except for the slot `k`.
"""
function _get!(f, d, k)
    y = get(d, k, KEY_SENTINEL)
    if y === KEY_SENTINEL
        y = f()
        @_assert !haskey(d, k)
        d[k] = y
        return y
    else
        return y
    end
end

ext: jl


text:
# Bidiagonal matrices
type Bidiagonal{T} <: AbstractMatrix{T}
    dv::AbstractVector{T} # diagonal
    ev::AbstractVector{T} # sub/super diagonal
    isupper::Bool # is upper bidiagonal (true) or lower (false)
    function Bidiagonal{T}(dv::AbstractVector{T}, ev::AbstractVector{T}, isupper::Bool)
        length(ev)==length(dv)-1 ? new(dv, ev, isupper) : throw(DimensionMismatch(""))
    end
end

Bidiagonal{T}(dv::AbstractVector{T}, ev::AbstractVector{T}, isupper::Bool)=Bidiagonal{T}(copy(dv), copy(ev), isupper)
Bidiagonal{T}(dv::AbstractVector{T}, ev::AbstractVector{T}) = error("Did you want an upper or lower Bidiagonal? Try again with an additional true (upper) or false (lower) argument.")

#Convert from BLAS uplo flag to boolean internal
Bidiagonal(dv::AbstractVector, ev::AbstractVector, uplo::BlasChar) = Bidiagonal(copy(dv), copy(ev), uplo=='U' ? true : uplo=='L' ? false : error("Bidiagonal can only be upper 'U' or lower 'L' but you said '$uplo''"))

function Bidiagonal{Td,Te}(dv::AbstractVector{Td}, ev::AbstractVector{Te}, isupper::Bool)
    T = promote_type(Td,Te)
    Bidiagonal(convert(Vector{T}, dv), convert(Vector{T}, ev), isupper)
end

Bidiagonal(A::AbstractMatrix, isupper::Bool)=Bidiagonal(diag(A), diag(A, isupper?1:-1), isupper)

#Converting from Bidiagonal to dense Matrix
full{T}(M::Bidiagonal{T}) = convert(Matrix{T}, M)
convert{T}(::Type{Matrix{T}}, A::Bidiagonal{T})=diagm(A.dv) + diagm(A.ev, A.isupper?1:-1)
promote_rule{T,S}(::Type{Matrix{T}}, ::Type{Bidiagonal{S}})=Matrix{promote_type(T,S)}

#Converting from Bidiagonal to Tridiagonal
Tridiagonal{T}(M::Bidiagonal{T}) = convert(Tridiagonal{T}, M)
function convert{T}(::Type{Tridiagonal{T}}, A::Bidiagonal{T})
    z = zeros(T, size(A)[1]-1)
    A.isupper ? Tridiagonal(A.ev, A.dv, z) : Tridiagonal(z, A.dv, A.ev)
end
promote_rule{T,S}(::Type{Tridiagonal{T}}, ::Type{Bidiagonal{S}})=Tridiagonal{promote_type(T,S)}

#################
# BLAS routines #
#################

#Singular values
svdvals{T<:BlasReal}(M::Bidiagonal{T})=LAPACK.bdsdc!(M.isupper?'U':'L', 'N', copy(M.dv), copy(M.ev))
svd    {T<:BlasReal}(M::Bidiagonal{T})=LAPACK.bdsdc!(M.isupper?'U':'L', 'I', copy(M.dv), copy(M.ev))

####################
# Generic routines #
####################

function show(io::IO, M::Bidiagonal)
    println(io, summary(M), ":")
    print(io, " diag:")
    print_matrix(io, (M.dv)')
    print(io, M.isupper?"\nsuper:":"\n  sub:")
    print_matrix(io, (M.ev)')
end

size(M::Bidiagonal) = (length(M.dv), length(M.dv))
size(M::Bidiagonal, d::Integer) = d<1 ? error("dimension out of range") : (d<=2 ? length(M.dv) : 1)

#Elementary operations
for func in (conj, copy, round, iround)
    func(M::Bidiagonal) = Bidiagonal(func(M.dv), func(M.ev), M.isupper)
end

transpose(M::Bidiagonal) = Bidiagonal(M.dv, M.ev, !M.isupper)
ctranspose(M::Bidiagonal) = Bidiagonal(conj(M.dv), conj(M.ev), !M.isupper)

istriu(M::Bidiagonal) = M.isupper
istril(M::Bidiagonal) = !M.isupper

function diag{T}(M::Bidiagonal{T}, n::Integer=0)
    if n==0
        return M.dv
    elseif n==1
        return M.isupper ? M.ev : zeros(T, size(M,1)-1)
    elseif n==-1
        return M.isupper ? zeros(T, size(M,1)-1) : M.ev
    elseif -size(M,1)<n<size(M,1)
        return zeros(T, size(M,1)-abs(n))
    else
        throw(BoundsError())
    end
end

function +(A::Bidiagonal, B::Bidiagonal)
    if A.isupper==B.isupper
        Bidiagonal(A.dv+B.dv, A.ev+B.ev, A.isupper)
    else
        apply(Tridiagonal, A.isupper ? (B.ev,A.dv+B.dv,A.ev) : (A.ev,A.dv+B.dv,B.ev))
    end
end

function -(A::Bidiagonal, B::Bidiagonal)
    if A.isupper==B.isupper
        Bidiagonal(A.dv-B.dv, A.ev-B.ev, A.isupper)
    else
        apply(Tridiagonal, A.isupper ? (-B.ev,A.dv-B.dv,A.ev) : (A.ev,A.dv-B.dv,-B.ev))
    end
end

-(A::Bidiagonal)=Bidiagonal(-A.dv,-A.ev)
*(A::Bidiagonal, B::Number) = Bidiagonal(A.dv*B, A.ev*B, A.isupper)
*(B::Number, A::Bidiagonal) = A*B
/(A::Bidiagonal, B::Number) = Bidiagonal(A.dv/B, A.ev/B, A.isupper)
==(A::Bidiagonal, B::Bidiagonal) = (A.dv==B.dv) && (A.ev==B.ev) && (A.isupper==B.isupper)

SpecialMatrix = Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular)
*(A::SpecialMatrix, B::SpecialMatrix)=full(A)*full(B)

#Generic multiplication
for func in (:*, :Ac_mul_B, :A_mul_Bc, :/, :A_rdiv_Bc)
    @eval begin
        ($func){T}(A::Bidiagonal{T}, B::AbstractVector{T}) = ($func)(full(A), B)
        #($func){T}(A::AbstractArray{T}, B::Triangular{T}) = ($func)(full(A), B)
    end
end

#Linear solvers
function \{TA<:Number,Tb<:Number}(A::Union(Bidiagonal{TA},Triangular{TA}), b::AbstractVector{Tb})
    TAb = typeof(one(TA)/one(Tb))
    naivesub!(A, b, similar(b, TAb))
end
function \{TA<:Number,TB<:Number}(A::Union(Bidiagonal{TA},Triangular{TA}), B::AbstractMatrix{TB})
    TAB = typeof(one(TA)/one(TB))
    hcat([naivesub!(A, B[:,i], similar(B[:,i], TAB)) for i=1:size(B,2)]...)
end
A_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(A,b)
At_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(transpose(A),b)
Ac_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(ctranspose(A),b)
for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin
    function ($func)(A::Union(Bidiagonal,Triangular), B::AbstractMatrix)
        tmp = similar(B[:,1])
        n = size(B, 1)
        for i=1:size(B,2)
            copy!(tmp, 1, B, (i-1)*n+1, n)
            ($func)(A, tmp)
            copy!(B, (i-1)*n+1, tmp, 1, n) # Modify this when array view are implemented.
        end
        B
    end
end end
for func in (:A_ldiv_Bt!, :Ac_ldiv_Bt!, :At_ldiv_Bt!) @eval begin
    function ($func)(A::Union(Bidiagonal,Triangular), B::AbstractMatrix)
        tmp = similar(B[:,2])
        m, n = size(B)
        nm = n*m
        for i=1:size(B,1)
            copy!(tmp, 1, B, i:m:nm, n)
            ($func)(A, tmp)
            copy!(B, i:m:nm, tmp, 1, n)
        end
        B
    end
end end

#Generic solver using naive substitution
function naivesub!{T}(A::Bidiagonal{T}, b::AbstractVector, x::AbstractVector=b)
    N = size(A, 2)
    N==length(b)==length(x) || throw(DimensionMismatch(""))

    if !A.isupper #do forward substitution
        for j = 1:N
            x[j] = b[j]
            j>1 && (x[j] -= A.ev[j-1] * x[j-1])
            x[j]/= A.dv[j]==zero(T) ? throw(SingularException(j)) : A.dv[j]
        end
    else #do backward substitution
        for j = N:-1:1
            x[j] = b[j]
            j<N && (x[j] -= A.ev[j] * x[j+1])
            x[j]/= A.dv[j]==zero(T) ? throw(SingularException(j)) : A.dv[j]
        end
    end
    x
end

# Eigensystems
eigvals(M::Bidiagonal) = M.dv
function eigvecs{T}(M::Bidiagonal{T})
    n = length(M.dv)
    Q=Array(T, n, n)
    blks = [0; find(x->x==0, M.ev); n]
    if M.isupper
        for idx_block=1:length(blks)-1, i=blks[idx_block]+1:blks[idx_block+1] #index of eigenvector
            v=zeros(T, n)
            v[blks[idx_block]+1] = one(T)
            for j=blks[idx_block]+1:i-1 #Starting from j=i, eigenvector elements will be 0
                v[j+1] = (M.dv[i]-M.dv[j])/M.ev[j] * v[j]
            end
            Q[:,i] = v/norm(v)
        end
    else
        for idx_block=1:length(blks)-1, i=blks[idx_block+1]:-1:blks[idx_block]+1 #index of eigenvector
            v=zeros(T, n)
            v[blks[idx_block+1]] = one(T)
            for j=(blks[idx_block+1]-1):-1:max(1,(i-1)) #Starting from j=i, eigenvector elements will be 0
                v[j] = (M.dv[i]-M.dv[j+1])/M.ev[j] * v[j+1]
            end
            Q[:,i] = v/norm(v)
        end
    end
    Q #Actually Triangular
end
eigfact(M::Bidiagonal) = Eigen(eigvals(M), eigvecs(M))

#Singular values
function svdfact(M::Bidiagonal, thin::Bool=true)
    U, S, V = svd(M)
    SVD(U, S, V')
end

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
#= TODO: fix conventions (abstract, type)

abstract GObjectI
#const
const GObject GObjectI

# Alternative object construction style. This would let us share constructors
# by creating const aliases: `const Z = GObject{:Z}`

type  GObjectAny{Name} <: GObjectI
    handle::Ptr{GObject}
    GObjectAny(handle::Ptr{GObject}) = (handle != C_NULL ? gc_ref(new(handle)) : error("Cannot construct $gname with a NULL pointer"))
end

=#

#type GtkWidgetAny{T} <: GtkWidgetI
#    handle::Ptr{GObject}
#    GtkWidgetAny(handle::Ptr{GObject}) = gc_ref(new(handle))
#end
#type GtkContainerAny{T} <: GtkContainerI
#    handle::Ptr{GObject}
#    GtkContainerAny(handle::Ptr{GObject}) = gc_ref(new(handle))
#end
#type GtkBinAny{T} <: GtkBinI
#    handle::Ptr{GObject}
#    GtkBinAny(handle::Ptr{GObject}) = gc_ref(new(handle))
#end
#type GtkBoxAny{T} <: GtkBoxI
#    handle::Ptr{GObject}
#    GtkBoxAny(handle::Ptr{GObject}) = gc_ref(new(handle))
#end

macro GType(gname)
    if isa(gname,Expr)
        @assert(gname.head == :comparison && length(gname.args) == 3 && gname.args[2] == :<:, "invalid GType expr")
        super = gname.args[3]
        gname = gname.args[1]
    else
        super = :GObject
    end
    gname = gname::Symbol
    quote
        type $(esc(gname)) <: $(esc(symbol(string(super,'I'))))
            handle::Ptr{GObjectI}
            $(esc(gname))(handle::Ptr{GObjectI}) = (handle != C_NULL ? gc_ref(new(handle)) : error("Cannot construct $gname with a NULL pointer"))
  end
end

macro quark_str(q)
    :( ccall((:g_quark_from_string, libglib), Uint32, (Ptr{Uint8},), bytestring($q)) )
end
const jlref_quark = quark"julia_ref"

# All GtkWidgets are expected to have a 'handle' field
# of type Ptr{GObjectI} corresponding to the Gtk object
# and an 'all' field which has type GdkRectangle
# corresponding to the rectangle allocated to the object,
# or to override the size, width, and height methods
convert(::Type{Ptr{GObjectI}},w::GObjectI) = w.handle
convert{T<:GObjectI}(::Type{T},w::Ptr{T}) = convert(T,convert(Ptr{GObjectI},w))
function convert{T<:GObjectI}(::Type{T},w::Ptr{GObjectI})
    x = ccall((:g_object_get_qdata, libgobject), Ptr{GObjectI}, (Ptr{GObjectI},Uint32), w, jlref_quark)
    x == C_NULL && error("GObject didn't have a corresponding Julia object")
    unsafe_pointer_to_objref(x)::T
end
eltype{T<:GObjectI}(::GSList{T}) = T

show(io::IO, w::GObjectI) = print(io,typeof(w))

### Miscellaneous types
#const
#const  Enum  Int32
baremodule GConnectFlags
    const AFTER = 1
    const SWAPPED = 2
    get(s::Symbol) =
        if s === :after
            AFTER
        elseif s === :swapped
            SWAPPED
        else
            Main.Base.error(Main.Base.string("invalid GConnectFlag ",s))
        end
end

### Garbage collection [prevention]
const gc_preserve = ObjectIdDict() # reference counted closures
function gc_ref(x::ANY)
    global gc_preserve
    gc_preserve[x] = (get(gc_preserve, x, 0)::Int)+1
    x
end
function gc_unref(x::ANY)
    global gc_preserve
    count = get(gc_preserve, x, 0)::Int-1
    if count <= 0
        delete!(gc_preserve, x)
    end
    nothing
end
gc_ref_closure{T}(x::T) = (gc_ref(x);cfunction(gc_unref, Void, (T, Ptr{Void})))
gc_unref(x::Any, ::Ptr{Void}) = gc_unref(x)

#= const gc_preserve_gtk = WeakKeyDict{GObjectI,Union(Bool,GObjectI)}() # gtk objects =#
function gc_ref{T<:GObjectI}(x::T)
    global gc_preserve_gtk
    addref = function()
        ccall((:g_object_ref_sink,libgobject),Ptr{GObjectI},(Ptr{GObjectI},),x)
        finalizer(x,function(x)
                global gc_preserve_gtk
                gc_preserve_gtk[x] = x # convert to a strong-reference
                ccall((:g_object_unref,libgobject),Void,(Ptr{GObjectI},),x) # may clear the strong reference
            end)
        gc_preserve_gtk[x] = true # record the existence of the object, but allow the finalizer
    end
    ref = get(gc_preserve_gtk,x,nothing)
    if isa(ref,Nothing)
        ccall((:g_object_set_qdata_full, libgobject), Void,
            (Ptr{GObjectI}, Uint32, Any, Ptr{Void}), x, jlref_quark, x,
            cfunction(gc_unref, Void, (T,))) # add a circular reference to the Julia object in the GObjectI
        addref()
    elseif !isa(ref,WeakRef)
        # oops, we previously deleted the link, but now it's back
        addref()
    else
        # already gc-protected, nothing to do
    end
    x
end

function gc_unref_weak(x::GObjectI)
    # this strongly destroys and invalidates the object
    # it is intended to be called by Gtk, not in user code function
    # note: this may be called multiple times by Gtk
    x.handle = C_NULL
    global gc_preserve_gtk
    delete!(gc_preserve_gtk, x)
    nothing
end
function gc_unref(x::GObjectI)
    # this strongly destroys and invalidates the object
    # it is intended to be called by Gtk, not in user code function
    ref = ccall((:g_object_get_qdata,libgobject),Ptr{Void},(Ptr{GObjectI},Uint32),x,jlref_quark)
    if ref != C_NULL && x !== unsafe_pointer_to_objref(ref)
        # We got called because we are no longer the default object for this handle, but we are still alive
        warn("Duplicate Julia object creation detected for GObject")
        ccall((:g_object_weak_ref,libgobject),Void,(Ptr{GObjectI},Ptr{Void},Any),x,cfunction(gc_unref_weak,Void,(typeof(x),)),x)
    else
        ccall((:g_object_steal_qdata,libgobject),Any,(Ptr{GObjectI},Uint32),x,jlref_quark)
        gc_unref_weak(x)
    end
    nothing
end
gc_unref(::Ptr{GObjectI}, x::GObjectI) = gc_unref(x)
gc_ref_closure(x::GObjectI) = C_NULL

baremodule GtkJustification
    const LEFT   = 0
    const RIGHT  = 1
    const CENTER = 2
    const FILL   = 3
end

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
"""
Author: Eric Ford and collaborators (see header of each file)
Created: August 2020
Contact: https://github.com/eford/RvSpectML.jl
"""

"""
[**RvSpectML.jl**](https://github.com/eford/RvSpectML.jl) is a Julia package for measuring radial velocities (RVs) from stellar spectroscopic timeseries via machine learning (ML).

"""
module RvSpectML

# Packages from RvSpectML ecosystem
using RvSpectMLBase
export speed_of_light_mps
export calc_normalization, normalize_spectrum!, get_λ_range
export make_chunk_list, make_orders_into_chunks
export make_chunk_list_timeseries, make_order_list_timeseries
export make_grid_for_chunk, filter_bad_chunks

using EchelleInstruments
export EchelleInstruments

using EchelleCCFs
export EchelleCCFs
CCFs = EchelleCCFs

using Scalpels
export Scalpels

#using Experimental
#using RvSpectMLPlots

# Packages we know we'll use in many places
using LinearAlgebra, Statistics
using DataFrames, Query
using Dates
#using AstroLib  # For solar position
include("astrolib/astrolib.jl")
# Packages that are being used and likely can be shared
# using Distributions, Interpolations, MultivariateStats, PDMats
# Packages we might use soon
#using Optim
#using Stheno, TemporalGPs
# Let's try to keep Plots out the main package to reduce compilation times when not plotting.
#using  Plots

include("pipeline/pipeline.jl")
using .Pipeline
export Pipeline
export extract_orders, prepare_line_list, ccf_total, ccf_orders, calc_rvs_from_ccf_total

#=
export PipelinePlan
export make_plot,  save_plot, save_data, need_to, has_cache, read_cache, set_cache!  # Query pipeline
export need_to!, dont_need_to!, reset_all_needs!             # Write to pipeline
export make_plot!, dont_make_plot!, make_all_plots!,  make_no_plots!
=#

include("interp/interp.jl")
#=
export TemporalGPInterpolation
export construct_gp_posterior, gp_marginal, predict_gp
export predict_mean, predict_deriv, predict_deriv2, predict_mean_and_deriv, predict_mean_and_derivs
=#

include("line_finder/line_finder.jl")
export LineFinderPlan

include("line_shapes/get_line_shapes.jl")

# alg.jl  will be responsible for exporting its own types, modules & functions
#include("alg/alg.jl")
# For now treating parts individually
include("alg/project_flux_common_wavelengths.jl")
include("alg/make_template_spectrum.jl")
#export interp_chunk_to_shifted_grid_gp_temporal
include("alg/combine_obs.jl")
export bin_times_consecutive, bin_rvs_consecutive, bin_spectra_consecutive
export bin_times_nightly, bin_rvs_nightly, bin_times_and_rvs_nightly, bin_spectra_nightly, rms_rvs_within_night
export bin_times_max_Δt, bin_rvs_max_Δt,  bin_times_and_rvs_max_Δt, bin_spectra_max_Δt

include("alg/dcpca.jl")
#export DCPCA
using .DCPCA
export clean_rvs_dcpca, calc_sigma_pca_scores

#include("alg/ppcca.jl")
#include("alg/rvs_from_gp_pairs.jl")

include("util/spectra.jl")
export calc_depth_and_expected_rv_precission, calc_formal_rv_precission
include("util/earth_mover_distance.jl")
export earth_mover_distance, wasserstein_distance_presorted, wasserstein_distance_presorted_common_x
include("util/sun.jl")
export calc_solar_alt

end


text:
using Plots;
using Latexify;
using Random;
using LinearAlgebra;
pyplot()

# range(0,2π,50)

# 定义一些常量、工具函数（添加噪声、计算均方误差）
MAX_ORDER = 9                         # 阶数
MAX_DATA_NUM = 20                     #
K = 10                                # K-折交叉验证（k-fold cross-validation）
MAX_NOISE = 0.2                       #
α = 0.05                              # 学习率
lr_decay = 0                          # 学习率权重衰减
Momentum = 0.9                        # 动量，为了快速收敛
MAX_ITER = 10000                      # 最大迭代次数
BATCH_SIZE = 25                       # 批次
rng = MersenneTwister(1234);
noise(x) = x + MAX_NOISE * randn(rng, Float64, size(x))

f(x) = sin(2π * x)

# 损失函数 无正则/有正则
E(y_real,y_predict) = (sum([Δy^2 for Δy in y_predict .- y_real])) / length(y_real)

# 使用开发集，验证集，测试集求解
# GetW(X_train,Y_train,λ)
function Calc(X_r, Y_r, GetW::Function)
    λs = []
    # 训练数据集（train dataset）：用来构建机器学习模型
    # 验证数据集（validation dataset）：辅助构建模型，用于在构建过程中评估模型，为模型提供无偏估计，进而调整模型超参数
    # 测试数据集（test dataset）：用来评估训练好的最终模型的性能
    train = 1:round(Integer, 0.6MAX_DATA_NUM)
    valid = round(Integer, 0.6MAX_DATA_NUM) + 1:round(Integer, 0.8MAX_DATA_NUM)
    test  = round(Integer, 0.8MAX_DATA_NUM) + 1:round(Integer, MAX_DATA_NUM)

    X_train, Y_train = X_r[train,:], Y_r[train]
    X_valid, Y_valid = X_r[valid,:], Y_r[valid]
    X_test, Y_test  = X_r[test,:], Y_r[test]

    W_hat = zeros(MAX_ORDER + 1, 1)
    lnλ = 0 # 正则项系数的自然对数
    err = E(Y_valid, X_valid * W_hat)
    for t_lnλ in -30:0.01:10
        W = GetW(X_train, Y_train, exp(t_lnλ))
        new_error = E(Y_valid, X_valid * W)
        if new_error <= err
            err = new_error
            lnλ = t_lnλ
            W_hat = W
        end
        push!(λs, t_lnλ => new_error)
    end
    return W_hat, lnλ, λs, E(Y_test, X_test * W_hat)
end

# 用于绘图
function W_Plot(X, Y, W)
    x_Plot = 0:0.01:1
    X_Plot = reshape([xd^i for i in MAX_ORDER:-1:0 for xd in x_Plot], length(x_Plot), MAX_ORDER + 1)
    Y_Plot = X_Plot * W

    scatter(X, Y, label="data")
    plot!(x_Plot, f.(x_Plot), label="real(x)")
    plot!(x_Plot, Y_Plot, label="predict(x)")
end

# 生成标签数据并打乱
x_s = shuffle(range(0, stop=1, length=MAX_DATA_NUM))
X_s = reshape([xd^i for i in MAX_ORDER:-1:0 for xd in x_s], length(x_s), MAX_ORDER + 1)
Y_s = noise(f.(x_s))

train = 1:round(Integer, 0.8MAX_DATA_NUM)
test = round(Integer, 0.8MAX_DATA_NUM) + 1:round(Integer, MAX_DATA_NUM)

X_train, Y_train = X_s[train,:], Y_s[train]
X_test, Y_test = X_s[test,:], Y_s[test]

# 解析解 无正则项
# $E_{\hat{W}} = \frac{1}{length(X)} ∑ (y - x\hat{w}))^2$
# $E_{\hat{W}} = \frac{1}{length(X)} (y - X\hat{W})^T(y - X\hat{W})$
# $\frac{∂E_{\hat{W}}}{∂\hat{W}} = \frac{1}{length(X)} 2X^T(X\hat{W} - y)$
# $\frac{∂E_{\hat{W}}}{∂\hat{W}} = 0$
#      ==> $\hat{W} = (X^TX)^{-1}X^Ty$

W_hat = (X_train' * X_train)^-1 * X_train' * Y_train

W_Plot(x_s,Y_s,W_hat)
savefig("ResolveNoLambda:$(E(Y_test, X_test * W_hat)).png")

# 解析解 有正则项
# $E_{\hat{W}} = \frac{1}{length(X)} [ ∑(x\hat{w} -y)^2 + λ∑\hat{w}^2 ]$
# $E_{\hat{W}} = \frac{1}{length(X)} [ (X\hat{W} - y)^T(X\hat{W} - y) + λ\hat{W}^T * \hat{W}]$
# $\frac{∂E_{\hat{W}}}{∂\hat{W}} = \frac{1}{length(X)} [2X^T(X\hat{W} - y) + 2λ\hat{W}]$
# $\frac{∂E_{\hat{W}}}{∂\hat{W}} = 0$
#      ==> $\hat{W} = (X^T * X + λI )^{-1}X^Ty$

GetW(X_train,Y_train,λ) = (X_train' * X_train + Matrix(I, MAX_ORDER + 1, MAX_ORDER + 1) * λ)^-1 * X_train' * Y_train

W_hat, lnλ, λs, err = Calc(X_s, Y_s, GetW)

plot([x[1] for x in λs],[y[2] for y in λs])
savefig("ResolveWithLambda-lambda-err.png")
W_Plot(x_s,Y_s,W_hat)
savefig("ResolveWithLambda:$(E(Y_test, X_test * W_hat)).png")

# 梯度下降 无正则项
# $\hat{W} = \hat{W} - α\frac{∂E_{\hat{W}}}{∂\hat{W}}$
W_hat = ones(MAX_ORDER + 1)
v = zeros(MAX_ORDER + 1)

for i = 1:MAX_ITER
    global W_hat, v
    ∂E_W = (X_train' * (X_train * W_hat - Y_train)) * (2 / MAX_DATA_NUM)
    lr_i = α / (1 + lr_decay * i)
    v = Momentum * v - lr_i * ∂E_W
    W_hat += v
end

W_Plot(x_s,Y_s,W_hat)
savefig("SGDNoLambda:$(E(Y_test, X_test * W_hat)).png")

# 梯度下降 有正则项
# $E_{\hat{W}} = \frac{1}{length(X)} [ ∑(x\hat{w} -y)^2 + λ∑\hat{w}^2 ]$
# $E_{\hat{W}} = \frac{1}{length(X)} [ (X\hat{W} - y)^T(X\hat{W} - y) + λ\hat{W}^T * \hat{W}]$
# $\frac{∂E_{\hat{W}}}{∂\hat{W}} = \frac{1}{length(X)} [2X^T(X\hat{W} - y) + 2λ\hat{W}]$

function GetW(X_train, Y_train, λ)
    W_hat = zeros(MAX_ORDER + 1)
    v = zeros(MAX_ORDER + 1)

    for i = 1:MAX_ITER
        ∂E_W = (X_train' * (X_train * W_hat - Y_train) + λ * W_hat) * (2 / MAX_DATA_NUM)
        lr_i = α / (1 + lr_decay * i)
        v = Momentum * v - lr_i * ∂E_W
        W_hat += v
    end
    return W_hat
end

W_hat, lnλ, λs, err = Calc(X_s, Y_s, GetW)

plot([x[1] for x in λs],[y[2] for y in λs])
savefig("SGDWithLambda-lambda-err.png")
W_Plot(x_s,Y_s,W_hat)
savefig("SGDWithLambda:$(E(Y_test, X_test * W_hat)).png")

# 随机梯度下降 无正则项
# $\hat{W} = \hat{W} - α\frac{∂E_{\hat{W}}}{∂\hat{W}}$

W_hat = zeros(MAX_ORDER + 1)
v = zeros(MAX_ORDER + 1)

for i = 1:MAX_ITER
    global W_hat, v
    p = rand(1:length(Y_train), BATCH_SIZE)
    ∂E_W = (X_train[p,:]' * (X_train[p,:] * W_hat - Y_train[p])) * (2 / BATCH_SIZE)
    lr_i = α / (1 + lr_decay * i)
    v = Momentum * v - lr_i * ∂E_W
    W_hat += v
end

W_Plot(x_s,Y_s,W_hat)
E(Y_test,X_test * W_hat)

# 随机梯度下降 有正则项
# $E_{\hat{W}} = \frac{1}{length(X)} [ ∑(x\hat{w} -y)^2 + λ∑\hat{w}^2 ]$
# $E_{\hat{W}} = \frac{1}{length(X)} [ (X\hat{W} - y)^T(X\hat{W} - y) + λ\hat{W}^T * \hat{W}]$
# $\frac{∂E_{\hat{W}}}{∂\hat{W}} = \frac{1}{length(X)} [2X^T(X\hat{W} - y) + 2λ\hat{W}]$

function GetW(X_train, Y_train, λ)
    W_hat = zeros(MAX_ORDER + 1)
    v = zeros(MAX_ORDER + 1)

    for i = 1:MAX_ITER
        p = rand(1:length(Y_train), BATCH_SIZE)
        ∂E_W = (X_train[p,:]' * (X_train[p,:] * W_hat - Y_train[p]) + λ * W_hat) * (2 / BATCH_SIZE)
        lr_i = α / (1 + lr_decay * i)
        v = Momentum * v - lr_i * ∂E_W
        W_hat += v
    end
    return W_hat
end

W_hat, lnλ, λs, err = Calc(X_s, Y_s, GetW)

plot([x[1] for x in λs],[y[2] for y in λs])
W_Plot(x_s,Y_s,W_hat)

# 共轭梯度下降 无正则项
# $\frac{∂E_{\hat{W}}}{∂\hat{W}} = \frac{1}{length(X)} 2X^T(X\hat{W} - y)$
# $X^T X \hat{W} = X^T Y$
#    $A  X  = b$

A = X_train' * X_train
b = X_train' * Y_train
W_hat = ones(MAX_ORDER + 1)

rk = b - A * W_hat
pk = rk

while true
    global W_hat, rk, pk
    αk = rk' * rk / (pk' * A * pk)
    W_hat = W_hat + αk * pk
    rk2 = rk -  αk * A * pk

    rk_l2 = sum([x^2 for x in rk2])
    if rk_l2 < 1e-3 || isnan(rk_l2)
        break
    end

    Βk = rk2' * rk2 /  (rk' * rk)
    pk2 = rk2 + Βk * pk

    rk = rk2
    pk = pk2
end

W_Plot(x_s,Y_s,W_hat)
savefig("ConjugationNoLambda:$(E(Y_test, X_test * W_hat)).png")

# 共轭梯度下降 有正则项
# $E_{\hat{W}} = \frac{1}{length(X)} [ ∑(x\hat{w} -y)^2 + λ∑\hat{w}^2 ]$
# $E_{\hat{W}} = \frac{1}{length(X)} [ (X\hat{W} - y)^T(X\hat{W} - y) + λ\hat{W}^T * \hat{W}]$
# $\frac{∂E_{\hat{W}}}{∂\hat{W}} = \frac{1}{length(X)} [2X^T(X\hat{W} - y) + 2λ\hat{W}]$
# $\frac{∂E_{\hat{W}}}{∂\hat{W}} = 0$
#      ==> $\hat{W} = (X^T * X + λI )^{-1}X^Ty$
#      ==> $(X^T * X + λI ) \hat{W} = X^Ty$
#                    $A  X  = b$

function GetW(X_train, Y_train, λ)
    A = X_train' * X_train + Matrix(I, MAX_ORDER + 1, MAX_ORDER + 1) * λ
    b = X_train' * Y_train
    W_hat = ones(MAX_ORDER + 1)

    rk = b - A * W_hat
    pk = rk

    while true
        αk = rk' * rk / (pk' * A * pk)
        W_hat = W_hat + αk * pk
        rk2 = rk -  αk * A * pk

        rk_l2 = sum([x^2 for x in rk2])
        if rk_l2 < 0.01 || isnan(rk_l2)
            break
        end

        Βk = rk2' * rk2 /  (rk' * rk)
        pk2 = rk2 + Βk * pk

        rk = rk2
        pk = pk2
    end
    return W_hat
end

W_hat, lnλ, λs, err = Calc(X_s, Y_s, GetW)

plot([x[1] for x in λs],[y[2] for y in λs])
savefig("ConjugationWithLambda-Lambda-err.png")
W_Plot(x_s,Y_s,W_hat)
savefig("ConjugationWithLambda:$(E(Y_test, X_test * W_hat)).png")

max_stars_repo_name: AlongWY/HIT-WorkShop

max_issues_repo_name: AlongWY/HIT-WorkShop

max_forks_repo_name: AlongWY/HIT-WorkShop


text:
# This file holds the definition of the functions pertaining to the BaseNode
#   type.

export  setNodeName!,
        setNodeTarget!,
        setNodeAttritionScheme!,
        addNodeRequirement!,
        removeNodeRequirement!,
        clearNodeRequirements!,
        setNodeRequirements!

"""
```
setNodeName!(
    node::BaseNode,
    name::String )
```
This function sets the name of the node `node` to `name`.

This function returns `true`, indication that the node name is set successfully.
"""
function setNodeName!( node::BaseNode, name::String )::Bool

    node.name = name
    return true

end  # setNodeName!( node, name )

"""
```
setNodeTarget!(
    node::BaseNode,
    target::Integer )
```
This function sets the target number of personnel members in node `node` to
`target`. If the number is less than zero, it means there's no target.

This function returns `true`, indicating the node population target is successfully set.
"""
function setNodeTarget!( node::BaseNode, target::Integer )::Bool

    node.target = max( target, -1 )
    return true

end  # setNodeTarget!( node, target )

"""
```
setNodeAttritionScheme!(
    node::BaseNode,
    attrition::String )
```
This function sets the name of the attrition scheme of the node `node` to `attrition`. The user must make sure that an attrition scheme with this name exists in the simulation.

This function returns `true`, indicating the attrition scheme is successfully set.
"""
function setNodeAttritionScheme!( node::BaseNode, attrition::String )::Bool

    node.attrition = lowercase( attrition ) ∈ ["", "default"] ? "default" :
        attrition
    return true

end  # setNodeAttritionScheme!( node, attrition )

"""
```
addNodeRequirement!(
    node::BaseNode,
    attrVals::NTuple{2,String}... )
```
This function adds the attribute requirements listed in `attrVals` to the node `node`. If attributes in the list already have a requirement on them, it gets overwritten.

If there are multiple entries for the same attribute, this function issues a warning and doesn't make any changes.

This function returns `true` if the requirements have been added successfully, and `false` if the list contained duplicates.
"""
function addNodeRequirement!( node::BaseNode,
    attrVals::NTuple{2,String}... )::Bool

    newAttrs = map( attrVal -> attrVal[1], collect( attrVals ) )

    if length( newAttrs ) != length( unique( newAttrs ) )
        @warn "Duplicate attribute entries in the requirement list, not making any changes."
        return false
    end  # if length( newAttrs ) != length( unique( newAttrs ) )

    newVals = map( attrVal -> attrVal[2], collect( attrVals ) )

    for ii in eachindex( newAttrs )
        node.requirements[newAttrs[ii]] = newVals[ii]
    end  # for ii in eachindex( newAttrs )

    return true

end  # addNodeRequirement!( node, attriVals )

"""
```
addNodeRequirement!(
    node::BaseNode,
    attribute::String,
    value::String )
```
This function adds the requirement that the attribute `attribute` must have the value `value` to qualify for the node `node`. If the node already has a requirement on that attribute, it gets overwritten.

This function returns `true`, indicating the requirement has been successfully added.
"""
addNodeRequirement!( node::BaseNode, attribute::String, value::String )::Bool =
    addNodeRequirement!( node, (attribute, value) )

"""
```
removeNodeRequirement!(
    node::BaseNode,
    attributes::String... )
```
This function removes the requirements on the attributes listed in `attributes` from the node `node`.

This function returns `true` if any requirements have been successfully removed, and `false` if the nodes didn't have any requirements on the given attributes.
"""
function removeNodeRequirement!( node::BaseNode, attributes::String... )::Bool

    hasRequirement = collect( haskey.( Ref( node.requirements ), attributes ) )

    if !any( hasRequirement )
        return false
    end  # if !any( hasRequirement )

    delete!.( Ref( node.requirements ), attributes[hasRequirement] )
    return true

end  # removeRequirement!( node, attributes )

"""
```
clearNodeRequirements!( node::BaseNode )
```
This function clears all the requirements for node `node`.

This function returns `true`, indicating all requirements on the node have been successfully cleared.
"""
function clearNodeRequirements!( node::BaseNode )::Bool

    empty!( node.requirements )
    return true

end  # clearNodeRequirements!( node )

"""
```
setNodeRequirements!(
    node::BaseNode,
    attrVals::Dict{String,String} )
```
This function sets the requirements for the node `node` to the attribute/value combinations in `attrVals`.

This function returns `true`, indicating the requirements have been successfully set.
"""
function setNodeRequirements!( node::BaseNode,
    attrVals::Dict{String,String} )::Bool

    node.requirements = deepcopy( attrVals )
    return true

end  # setNodeRequirements!( node, attrVals )

"""
```
setNodeRequirements!(
    node::BaseNode,
    attrVals::NTuple{2,String}... )
```
This function sets the requirements for the node `node` to the attribute/value combinations in `attrVals`.

If there are multiple entries for the same attribute, this function isues a warning and doesn't make any changes.

This function returns `true` when the requirements have been set successfully, and `false` when the list contains duplicate entries.
"""
function setNodeRequirements!( node::BaseNode,
    attrVals::NTuple{2,String}... )::Bool

    newAttrs = map( attrVal -> attrVal[1], collect( attrVals ) )

    if length( newAttrs ) != length( unique( newAttrs ) )
        @warn "Duplicate attribute entries in the requirement list, not making any changes."
        return false
    end  # if length( newAttrs ) != length( unique( newAttrs ) )

    attrValDict = Dict{String,String}()

    for attrVal in attrVals
        attrValDict[attrVal[1]] = attrVal[2]
    end  # for attrVal in attrVals

    return setNodeRequirements!( node, attrValDict )

end  # setNodeRequirements!( node, attrVals )

"""
```
setNodeRequirements!(
    node::BaseNode,
    attributes::Vector{String},
    values::Vector{String} )
```
This function sets the requirements of the node `node` to the attributes and corresponding values in `attributes` resp. `values`.

If the two vectors differ in length or if there are duplicate entries in the list of attributes, the function will issue a warning and not make any changes.

This function returns `true` when the requirements have been set successfully, and `false` when the list contains duplicate entries or the attributes and values lists don't match up.
"""
function setNodeRequirements!( node::BaseNode, attributes::Vector{String}, values::Vector{String} )::Bool

    if length( attributes ) != length( values )
        @warn "Mismatched lengths of vector of attributes and values, not making any changes."
        return false
    end  # if length( values ) != length( weights )

    if length( attributes ) != length( unique( attributes ) )
        @warn "Duplicate entries in the attributes list, not making any changes."
        return false
    end  # if length( values ) != length( unique( values ) )

    attrValDict = Dict{String,String}()

    for ii in eachindex( attributes )
        attrValDict[attributes[ii]] = values[ii]
    end  # for ii in eachindex( attributes )

    return setNodeRequirements!( node, attrValDict )

end  # setNodeRequirements!( node, attributes, values )

function Base.show( io::IO, node::BaseNode )::Nothing

    print( io, "Node: ", node.name )
    print( io, "\n  Associated attrition scheme: ", node.attrition )

    if node.target >= 0
        print( io, "\n  Population target of ", node.target,
            " personnel members in node." )
    else
        print( io, "\n  Node has no population target." )
    end  # if node.target >= 0

    if isempty( node.requirements )
        print( io, "\n  Node '", node.name, "' has no attribute requirements." )
        return
    end  # if isempty( node.requirements )

    print( io, "\n  Requirements" )

    for attribute in keys( node.requirements )
        print( io, "\n    ", attribute, " = ", node.requirements[attribute] )
    end  # for attr in keys( node.requirements )

    return

end  # show( io, node )

include( joinpath( privPath, "basenode.jl" ) )


text:
module Mutabilities

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end Mutabilities

export freeze,
    freezeindex,
    freezeproperties,
    freezevalue,
    melt,
    meltindex,
    meltproperties,
    meltvalue,
    move!,
    readonly

import SplittablesBase
using ConstructionBase: constructorof, setproperties
using Requires: @require

include("utils.jl")
include("interface.jl")
include("baseimpl.jl")
include("arrays.jl")
include("dicts.jl")
include("sets.jl")
include("datatypes.jl")

const AbstractReadOnlyContainers =
    Union{AbstractReadOnlyArray,AbstractReadOnlyDict,AbstractReadOnlySet}

include("splittables.jl")

function __init__()
    @require StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" begin
        include("staticarrays.jl")
    end
    @require StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" begin
        include("structarrays.jl")
    end
end

end # module


text:
module Lookup3
import JLD2: jlunsafe_load

# Original source at http://www.burtleburtle.net/bob/c/lookup3.c

rot(x::UInt32, k) = xor(((x)<<(k)), ((x)>>(32-(k))))

# -------------------------------------------------------------------------------
# mix -- mix 3 32-bit values reversibly.

# This is reversible, so any information in (a,b,c) before mix() is
# still in (a,b,c) after mix().

# If four pairs of (a,b,c) inputs are run through mix(), or through
# mix() in reverse, there are at least 32 bits of the output that
# are sometimes the same for one pair and different for another pair.
# This was tested for:
# * pairs that differed by one bit, by two bits, in any combination
#   of top bits of (a,b,c), or in any combination of bottom bits of
#   (a,b,c).
# * "differ" is defined as +, -, ^, or ~^.  For + and -, I transformed
#   the output delta to a Gray code (a^(a>>1)) so a string of 1's (as
#   is commonly produced by subtraction) look like a single 1-bit
#   difference.
# * the base values were pseudorandom, all zero but one bit set, or
#   all zero plus a counter that starts at zero.

# Some k values for my "a-=c; a^=rot(c,k); c+=b;" arrangement that
# satisfy this are
#     4  6  8 16 19  4
#     9 15  3 18 27 15
#    14  9  3  7 17  3
# Well, "9 15 3 18 27 15" didn't quite get 32 bits diffing
# for "differ" defined as + with a one-bit base and a two-bit delta.  I
# used http://burtleburtle.net/bob/hash/avalanche.html to choose
# the operations, constants, and arrangements of the variables.

# This does not achieve avalanche.  There are input bits of (a,b,c)
# that fail to affect some output bits of (a,b,c), especially of a.  The
# most thoroughly mixed value is c, but it doesn't really even achieve
# avalanche in c.

# This allows some parallelism.  Read-after-writes are good at doubling
# the number of bits affected, so the goal of mixing pulls in the opposite
# direction as the goal of parallelism.  I did what I could.  Rotates
# seem to cost as much as shifts on every machine I could lay my hands
# on, and rotates are much kinder to the top and bottom bits, so I used
# rotates.
# -------------------------------------------------------------------------------
@inline function mix(a, b, c)
  a -= c;  a = xor(a, rot(c, 4));  c += b
  b -= a;  b = xor(b, rot(a, 6));  a += c
  c -= b;  c = xor(c, rot(b, 8));  b += a
  a -= c;  a = xor(a, rot(c,16));  c += b
  b -= a;  b = xor(b, rot(a,19));  a += c
  c -= b;  c = xor(c, rot(b, 4));  b += a
  (a, b, c)
end

# # -------------------------------------------------------------------------------
# final -- final mixing of 3 32-bit values (a,b,c) into c

# Pairs of (a,b,c) values differing in only a few bits will usually
# produce values of c that look totally different.  This was tested for
# * pairs that differed by one bit, by two bits, in any combination
#   of top bits of (a,b,c), or in any combination of bottom bits of
#   (a,b,c).
# * "differ" is defined as +, -, ^, or ~^.  For + and -, I transformed
#   the output delta to a Gray code (a^(a>>1)) so a string of 1's (as
#   is commonly produced by subtraction) look like a single 1-bit
#   difference.
# * the base values were pseudorandom, all zero but one bit set, or
#   all zero plus a counter that starts at zero.

# These constants passed:
#  14 11 25 16 4 14 24
#  12 14 25 16 4 14 24
# and these came close:
#   4  8 15 26 3 22 24
#  10  8 15 26 3 22 24
#  11  8 15 26 3 22 24
# -------------------------------------------------------------------------------
@inline function final(a, b, c)
  c = xor(c, b); c -= rot(b,14)
  a = xor(a, c); a -= rot(c,11)
  b = xor(b, a); b -= rot(a,25)
  c = xor(c, b); c -= rot(b,16)
  a = xor(a, c); a -= rot(c,4)
  b = xor(b, a); b -= rot(a,14)
  c = xor(c, b); c -= rot(b,24)
  c
end

# -------------------------------------------------------------------------------
# hashlittle() -- hash a variable-length key into a 32-bit value
#   k       : the key (the unaligned variable-length array of bytes)
#   length  : the length of the key, counting by bytes
#   initval : can be any 4-byte value
# Returns a 32-bit value.  Every bit of the key affects every bit of
# the return value.  Two keys differing by one or two bits will have
# totally different hash values.

# The best hash table sizes are powers of 2.  There is no need to do
# mod a prime (mod is sooo slow!).  If you need less than 32 bits,
# use a bitmask.  For example, if you need only 10 bits, do
#   h = (h & hashmask(10));
# In which case, the hash table should have hashsize(10) elements.

# If you are hashing n strings (uint8_t **)k, do it like this:
#   for (i=0, h=0; i<n; ++i) h = hashlittle( k[i], len[i], h);

# By Bob Jenkins, 2006.  bob_jenkins@burtleburtle.net.  You may use this
# code any way you wish, private, educational, or commercial.  It's free.

# Use for hash table lookup, or anything where one collision in 2^^32 is
# acceptable.  Do NOT use for cryptographic purposes.
# -------------------------------------------------------------------------------
function hash(k::AbstractVector{UInt8}, istart::Integer=1, n::Integer=length(k), initval::UInt32=UInt32(0))
    n <= length(k) || throw(BoundsError())
    hash(pointer(k, istart), n, initval)
end

function hash(k::Ptr{UInt8}, n::Integer=length(k), initval::UInt32=UInt32(0))
    # Set up the internal state
    a = b = c = 0xdeadbeef + convert(UInt32, n) + initval

    # --------------- all but the last block: affect some 32 bits of (a,b,c)
    ptr = k
    # @inbounds while n > 12
    #     a += k[offset]
    #     a += convert(UInt32, k[offset+1])<<8
    #     a += convert(UInt32, k[offset+2])<<16
    #     a += convert(UInt32, k[offset+3])<<24
    #     @show a
    #     b += k[offset+4]
    #     b += convert(UInt32, k[offset+5])<<8
    #     b += convert(UInt32, k[offset+6])<<16
    #     b += convert(UInt32, k[offset+7])<<24
    #     @show b
    #     c += k[offset+8]
    #     c += convert(UInt32, k[offset+9])<<8
    #     c += convert(UInt32, k[offset+10])<<16
    #     c += convert(UInt32, k[offset+11])<<24
    #     @show c
    #     (a, b, c) = mix(a, b, c)
    #     n -= 12
    #     offset += 12
    # end
    @inbounds while n > 12
        a += jlunsafe_load(convert(Ptr{UInt32}, ptr))
        ptr += 4
        b += jlunsafe_load(convert(Ptr{UInt32}, ptr))
        ptr += 4
        c += jlunsafe_load(convert(Ptr{UInt32}, ptr))
        (a, b, c) = mix(a, b, c)
        ptr += 4
        n -= 12
    end

    # -------------------------------- last block: affect all 32 bits of (c)
    # @inbounds if n > 0
    #     n >= 12 && (c += convert(UInt32, k[offset+11])<<24)
    #     n >= 11 && (c += convert(UInt32, k[offset+10])<<16)
    #     n >= 10 && (c += convert(UInt32, k[offset+9])<<8)
    #     n >= 9  && (c += k[offset+8])
    #     n >= 8  && (b += convert(UInt32, k[offset+7])<<24)
    #     n >= 7  && (b += convert(UInt32, k[offset+6])<<16)
    #     n >= 6  && (b += convert(UInt32, k[offset+5])<<8)
    #     n >= 5  && (b += k[offset+4])
    #     n >= 4  && (a += convert(UInt32, k[offset+3])<<24)
    #     n >= 3  && (a += convert(UInt32, k[offset+2])<<16)
    #     n >= 2  && (a += convert(UInt32, k[offset+1])<<8)
    #     n >= 1  && (a += k[offset])
    #     c = final(a, b, c)
    # end
    @inbounds if n > 0
        if n == 12
            c += jlunsafe_load(convert(Ptr{UInt32}, ptr+8))
            @goto n8
        elseif n == 11
            c += UInt32(jlunsafe_load(Ptr{UInt8}(ptr+10)))<<16
            @goto n10
        elseif n == 10
            @label n10
            c += UInt32(jlunsafe_load(Ptr{UInt8}(ptr+9)))<<8
            @goto n9
        elseif n == 9
            @label n9
            c += jlunsafe_load(ptr+8)
            @goto n8
        elseif n == 8
            @label n8
            b += jlunsafe_load(convert(Ptr{UInt32}, ptr+4))
            @goto n4
        elseif n == 7
            @label n7
            b += UInt32(jlunsafe_load(Ptr{UInt8}(ptr+6)))<<16
            @goto n6
        elseif n == 6
            @label n6
            b += UInt32(jlunsafe_load(Ptr{UInt8}(ptr+5)))<<8
            @goto n5
        elseif n == 5
            @label n5
            b += jlunsafe_load(ptr+4)
            @goto n4
        elseif n == 4
            @label n4
            a += jlunsafe_load(convert(Ptr{UInt32}, ptr))
        elseif n == 3
            @label n3
            a += UInt32(jlunsafe_load(Ptr{UInt8}(ptr+2)))<<16
            @goto n2
        elseif n == 2
            @label n2
            a += UInt32(jlunsafe_load(Ptr{UInt8}(ptr+1)))<<8
            @goto n1
        elseif n == 1
            @label n1
            a += jlunsafe_load(ptr)
        end
        c = final(a, b, c)
    end
    c
end
end


text:
using Plots
using DelimitedFiles
using FileIO
pgfplotsx()
thetas = readdlm("thetas.txt")
p = plot(0:45, thetas, color=:black, ylims = [(0,2.5) (0, 9)], layout = (2,1), label=["\$\\theta_1\$" "\$\\theta_2\$"]);
vspan!(p, [[0,15] [0,15]],
        color = :black, alpha = 0.15,
        labels = ["20 x 24" "20 x 24"])
        #; "41 x 49" "55 x 65"
vspan!(p, [[15,30] [15,30]], color = :white, alpha = 0.15, labels = ["41 x 49" "41 x 49"]);
vspan!(p, [[30,45] [30,45]], color =  :blue, alpha = 0.15, labels = "55 x 65");

p

save("thetas.png", p)

@show thetas


text:
using MultipleViewGeometry, MultipleViewGeometry.ModuleMove, LinearAlgebra, Test

# Intrinsic and extrinsic parameters for the first camera.
𝐊₁ = zeros(3,3)
𝐊₁[1,1] = 10
𝐊₁[2,2] = 10
𝐊₁[3,3] = 1
𝐑₁ = rotxyz(deg2rad(10), deg2rad(15), deg2rad(45))
𝐭₁ = [-250.0, 0.0, 2500.0]

# Intrinsic and extrinsic parameters for the second camera.
𝐊₂ = zeros(3,3)
𝐊₂[1,1] = 5
𝐊₂[2,2] = 5
𝐊₂[3,3] = 1
𝐑₂ =  rotxyz(deg2rad(10), deg2rad(15), deg2rad(45))
𝐭₂ = [250.0,   0.0, 2500.0]

𝐅 = construct(FundamentalMatrix(),𝐊₁,𝐑₁,𝐭₁,𝐊₂,𝐑₂,𝐭₂)
𝐄 = construct(EssentialMatrix(),𝐅, 𝐊₁, 𝐊₂)

# Result 9.17 of R. Hartley and A. Zisserman, “Two-View Geometry,” Multiple View Geometry in Computer Vision
# A 3 by 3 matrix is an essential matrix if and only if two of its singular values
# are equal, and the third is zero.
U, S , V = svd(𝐄)

@test isapprox.(S[1], S[2]; atol = 1e-14)
@test isapprox.(S[3], 0.0; atol = 1e-10)


text:
using Documenter, Example

makedocs(modules = [Example], sitename = "Example2.jl")

deploydocs(repo = "github.com/mortenpi/Example2.jl.git")


text:
using InteractiveUtils

versioninfo()
#=
Multiple dispatch is a key feature of Julia, that we will explore
It helps make software fast, It also makes software extensible,
programmable, and downright fun to play with.

It may just hearld a breakthrough for parallel computation.
=#

# Let's define a new type to represent a Roman numerical.
mutable struct Roman
    n::Int
end

#For coding simplicity, we'll just deal with numbers between 0 and 9
Base.show(io::IO, r::Roman) = print(io, 'ⅰ' + (r.n - 1) % 10 ) # 'ⅰ' is a Unicode Roman numeral

@show Roman(6) # creating object
@show typeof.([5 5.0 Roman(5) "Five" '5'  5//1])

# We would like to display it nicely, in Roman numerals:
romans = [1, 2, 3, 4, 5, 6, 7, 8, 9]
@show typeof(romans)
@show Roman.(romans) # equivalent to map(Roman, romans)

import Base: +, *

+(a::Roman, b::Roman) = Roman(a.n + b.n)
*(a::Roman, b::Roman) = Roman(a.n * b.n)

@show Roman(3) + Roman(7)
@show Roman(7) * Roman(3)
@show Roman(3) ^ 2

#=
@show Roman(3) * 2
ERROR: LoadError: MethodError: no method matching *(::Roman, ::Int64)
=#

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function tryincrement(counter :: Dict{T, Int}, key :: T, amount :: Int = 1) where T
	haskey(counter, key) ? counter[key] += amount : counter[key] = amount
end

function setup()
	lines = readlines("inputD14.txt")
	pattern = lines[1]
	pairs = Dict{String, String}()

	for line in lines[3:end]
		left, right = split(line, " -> ")
		right = left[1] * right * left[2]
		pairs[left] = right
	end

	return (pattern, pairs)
end

function part1()
	# A brute-force solution, for reference
	# Won't work for part 2 because the size of the string gets far too big
	pattern, pairs = setup()

	for step in 1:10
		broken = [pattern[index] * pattern[index + 1] for index in 1:length(pattern)-1]
		broken = [get(pairs, b, b) for b in broken]

		pattern = join(b[1:end-1] for b in broken)
		pattern *= broken[end][end]
	end

	counter = Dict{Char, Int}()
	for char in pattern
		tryincrement(counter, char)
	end

	return maximum(values(counter)) - minimum(values(counter))
end

function part2()
	pattern, pairs = setup()
	letterpairs = Dict{String, Int}()
	totalcounts = Dict{Char, Int}()

	for char in pattern
		tryincrement(totalcounts, char)
	end

	for index in 1:length(pattern)-1
		tryincrement(letterpairs, pattern[index] * pattern[index + 1])
	end

	for step in 1:40
		newletterpairs = Dict{String, Int}()
		for letters in keys(letterpairs)
			result = get(pairs, letters, letters)
			amount = letterpairs[letters]
			if length(result) == 2
				tryincrement(newletterpairs, letters, amount)
			else # implies length(result) == 3
				tryincrement(newletterpairs, result[1:2], amount)
				tryincrement(newletterpairs, result[2:3], amount)
				tryincrement(totalcounts, result[2], amount)
			end
		end

		letterpairs = copy(newletterpairs)
	end

	return maximum(values(totalcounts)) - minimum(values(totalcounts))
end

println(part1())
println(part2())


text:
using SIMDTypes
using Documenter

DocMeta.setdocmeta!(SIMDTypes, :DocTestSetup, :(using SIMDTypes); recursive=true)

makedocs(;
    modules=[SIMDTypes],
    authors="Julia Computing, Inc. and Contributors",
    repo="https://github.com/JuliaSIMD/SIMDTypes.jl/blob/{commit}{path}#{line}",
    sitename="SIMDTypes.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaSIMD.github.io/SIMDTypes.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaSIMD/SIMDTypes.jl",
)


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

# Run a local search algorithm to optimise a nested compliance table given some objective function (objFn).
# The local search is simple hill climbing, and the neighbourhood of a nested compliance table
# involves swapping or inserting items (station indices) within/into the table.
# The local search stops when a local optimum (as evaluated by simulating) is found.
# Ambulances are assumed to be identical.

using JEMSS
using Random
using Statistics
using StatsFuns
using Base.Iterators: Stateful

global ObjVal = Float

function nestedCompTableLocalSearch!(sim::Simulation, nestedCompTables::Vector{NestedCompTable}; outputFolder::String = "")
	isdir(outputFolder) || mkdir(outputFolder)
	@assert(isdir(outputFolder))
	@assert(!isempty(nestedCompTables))

	# parameters
	global solFilename = "$outputFolder/solutions.csv" # final solutions from each local search
	global nestedCompTablesOutputFilename = "$outputFolder/nested_comp_tables.csv" # save final solutions from each search (saved as priority lists)
	global logFilename = "$outputFolder/log.csv" # log search progress to file
	global sense = :max # :min or :max; direction of optimisation for objective function
	global conf = 0.95 # statistical confidence level
	global doPrint = true
	# reps:
	global minNumReps = 2 # min number of reps to simulate each nested compliance table
	global maxNumRepsList = [2,5] # = vcat(10 * 2 .^ [0:10;], Inf) # can start with a small maximum, and increase each time local optimum found
	global maxNumRepsIncrease = 1 # = 20 # max number of reps to increase by at a time when comparing two nested compliance tables
	# parallel (multi-threading):
	global parallel = false
	global numThreads = 1

	# some parameter checks
	@assert(sense == :min || sense == :max)
	@assert(0 <= conf < 1) # should be between 0 and 1, but not equal to 1 otherwise confidence interval is infinite
	@assert(minNumReps >= 2) # cannot calculate confidence interval for less than 2 reps
	@assert(minNumReps <= sim.numReps)
	@assert(all(maxNumRepsList .>= minNumReps))
	@assert(issorted(maxNumRepsList, lt=<=)) # maxNumRepsList should be strictly increasing
	@assert(maxNumRepsIncrease >= 1)
	@assert(isa(parallel, Bool))
	@assert(numThreads >= 1)

	global nullObjVal = sense == :max ? -Inf : Inf

	# keep track of nested compliance tables tried, and their objective values per replication
	global nestedCompTablesObjVals = Dict{NestedCompTable, Vector{ObjVal}}() # nestedCompTablesObjVals[nestedCompTable][repIndex] gives sim replication objective value

	global nestedCompTablesPeriodStatsList = Dict{NestedCompTable, Vector{SimPeriodStats}}() # nestedCompTablesPeriodStatsList[nestedCompTable][repIndex] gives sim replication period stats

	# save the locally optimal nested compliance tables found for each search
	global nestedCompTableSols = Vector{NestedCompTable}() # nestedCompTableSols[i] is solution for nestedCompTables[i]

	global stationCapacities = [station.capacity for station in sim.stations] # stationCapacities[i] gives ambulance holding capacity of station i

	repeatedLocalSearch!(sim, nestedCompTables)
end

function objValLookup(nestedCompTable::NestedCompTable, repIndex::Int)::ObjVal
	global nestedCompTablesObjVals
	if haskey(nestedCompTablesObjVals, nestedCompTable) && length(nestedCompTablesObjVals[nestedCompTable]) >= repIndex
		return nestedCompTablesObjVals[nestedCompTable][repIndex]
	end
	return nullObjVal
end
function objValsLookup(nestedCompTable::NestedCompTable, repIndices::Vector{Int})::Vector{ObjVal}
	return [objValLookup(nestedCompTable, repIndex) for repIndex in repIndices]
end
function objValsLookup(nestedCompTable::NestedCompTable)::Vector{ObjVal}
	global nestedCompTablesObjVals
	return get(nestedCompTablesObjVals, nestedCompTable, ObjVal[])
end

function setNestedCompTableObjVal!(nestedCompTable::NestedCompTable, repIndex::Int, objVal::ObjVal)
	global nestedCompTablesObjVals
	objVals = get!(nestedCompTablesObjVals, nestedCompTable, ObjVal[])
	@assert(length(objVals) == repIndex - 1) # otherwise will overwrite existing value / will have missing values in objVals
	push!(objVals, objVal)
end

function periodStatsListLookup(nestedCompTable::NestedCompTable)::Vector{SimPeriodStats}
	global nestedCompTablesPeriodStatsList
	return get(nestedCompTablesPeriodStatsList, nestedCompTable, SimPeriodStats[])
end

function setNestedCompTablePeriodStats!(nestedCompTable::NestedCompTable, repIndex::Int, period::SimPeriodStats)
	global nestedCompTablesPeriodStatsList
	periods = get!(nestedCompTablesPeriodStatsList, nestedCompTable, SimPeriodStats[])
	@assert(length(periods) == repIndex - 1) # otherwise will overwrite existing value / will have missing values in periods
	push!(periods, period)
end
function setNestedCompTablePeriodStatsList!(nestedCompTable::NestedCompTable, repIndices::UnitRange{Int}, periods::Vector{SimPeriodStats})
	for (repIndex, period) in zip(repIndices, periods)
		setNestedCompTablePeriodStats!(nestedCompTable, repIndex, period)
	end
end

function resetNestedCompTablesPeriodStatsList!()
	global nestedCompTablesPeriodStatsList
	empty!(nestedCompTablesPeriodStatsList)
end

# Return mean and half-width of mean of x.
# Assumes that values in x are from population with normal distribution with unknown standard deviation.
function calcMeanAndHalfWidth(x::Vector{T}; conf::Float = conf)::Tuple{Float,Float} where T <: Real
	x = convert(Vector{Float}, x)
	return mean(x), tDistrHalfWidth(x; conf = conf)
end

# Objective value for a single period
function objFn(period::SimPeriodStats)::Float
	return period.call.numResponsesInTime / period.call.numCalls # fraction of calls reached in time
end

# Return objective values from multiple replications (repIndices).
# Uses look-up if possible.
# Mutates: sim.reps, nestedCompTablesObjVals
function simObjVals!(sim::Simulation, nestedCompTable::NestedCompTable, repIndices::UnitRange{Int})::Tuple{Vector{ObjVal}, Vector{SimPeriodStats}, Int}
	# look-up any saved data
	objValsFound = objValsLookup(nestedCompTable)
	periodsFound = periodStatsListLookup(nestedCompTable)
	n = min(length(objValsFound), length(periodsFound))
	lookupIndices= intersect(1:n, repIndices)
	simRepIndices = setdiff(repIndices, lookupIndices)

	objVals = ObjVal[]
	periods = SimPeriodStats[]
	if !isempty(simRepIndices)
		# need to simulate for simRepIndices
		global parallel, numThreads
		function runRep!(rep::Simulation)
			reset!(rep)
			applyNestedCompTable!(rep, nestedCompTable)
			simulate!(rep)
		end
		reps = sim.reps[simRepIndices]
		if parallel == true
			runParallel!(runRep!, reps...; numThreads = numThreads)
		else
			for rep in reps runRep!(rep) end
		end
		periods = getRepsPeriodStatsList(reps)
		objVals = objFn.(periods)
		for (i,repIndex) in enumerate(simRepIndices)
			if length(objValsFound) < repIndex setNestedCompTableObjVal!(nestedCompTable, repIndex, objVals[i]) end
			if length(periodsFound) < repIndex setNestedCompTablePeriodStats!(nestedCompTable, repIndex, periods[i]) end
		end
	end

	objVals = vcat(objValsFound[lookupIndices], objVals)
	periods = vcat(periodsFound[lookupIndices], periods)
	@assert(length(objVals) == length(repIndices))
	@assert(length(periods) == length(repIndices))

	return objVals, periods, length(lookupIndices)
end

function applyNestedCompTable!(sim::Simulation, nestedCompTable::NestedCompTable)
	@assert(!sim.used)
	for i = 1:sim.numAmbs
		setAmbStation!(sim, sim.ambulances[i], sim.stations[nestedCompTable[i]])
	end
	initCompTable!(sim, nestedCompTable)
end

function printNestedCompTable(nestedCompTable::NestedCompTable)
	for i = 1:length(nestedCompTable)
		println("item ", i, ": station ", nestedCompTable[i])
	end
end

function getPeriodsStatsDict(periods::Vector{SimPeriodStats})::Dict{String,Any}
	global periodDuration, simStatsKeys, simStatsEmpty
	@assert(all(p -> isapprox(p.duration, periodDuration), periods))
	stats = flatten(statsDictFromPeriodStatsList(periods))
	return merge(Dict(["$(key)_mean" => stats[key].mean for key in simStatsKeys]), Dict(["$(key)_halfWidth" => stats[key].halfWidth for key in simStatsKeys]))
end

global simStatsEmpty = flatten(statsDictFromPeriodStatsList(SimPeriodStats[]))
global simStatsKeys = collect(keys(simStatsEmpty))
global simStatsEmpty = merge(Dict(["$(key)_mean" => NaN for key in simStatsKeys]), Dict(["$(key)_halfWidth" => NaN for key in simStatsKeys]))

global solFileHeader = vcat(["search", "numReps", "objVal", "objValHalfWidth"], sort(collect(keys(simStatsEmpty)))) # ... and nestedCompTable
global logFileHeader = vcat(["search", "maxNumRepsListIndex", "maxNumReps", "iter", "move", "i", "j", "usedMove", "numReps", "numLookups", "searchDurationSeconds", "objVal", "objValHalfWidth", "bestObjVal", "objValDiff", "objValDiffHalfWidth"], sort(collect(keys(simStatsEmpty)))) # ... and nestedCompTable
global logFileDict = Dict{String,Any}([s => "" for s in logFileHeader])

# mutates: file
function fileWriteDlmLine!(file::IOStream, fileHeader::Vector{String}, data::Dict{String,T}, nestedCompTable::NestedCompTable) where T <: Any
	line = []
	for header in fileHeader
		push!(line, data[header])
	end
	line = vcat(line, nestedCompTable)
	writeDlmLine!(file, line...)
	flush(file)
end

# perform local search, starting at each of the nested compliance tables provided
# mutates: sim
function repeatedLocalSearch!(sim::Simulation, nestedCompTables::Vector{NestedCompTable})
	global doPrint

	# check that sim is using nested compliance table for move-up
	@assert(sim.moveUpData.useMoveUp == true && sim.moveUpData.moveUpModule == compTableModule)

	# open files for writing solution
	global solFilename, logFilename
	solFile = open(solFilename, "w")
	logFile = open(logFilename, "w")

	# write misc data to files
	global warmUpDuration, periodDuration, conf
	numSearches = length(nestedCompTables)
	miscTable = Table("miscData", ["numAmbs", "numStations", "numSearches", "conf"];
		rows = [[sim.numAmbs, sim.numStations, numSearches, conf]]) # note that numCalls may vary between sim replications
	repsTable = Table("repsData", ["numReps", "warmUpDuration", "periodDuration"];
		rows = [[sim.numReps, warmUpDuration, periodDuration]])
	writeTablesToFile!(logFile, [miscTable, repsTable])
	writeTablesToFile!(solFile, [miscTable, repsTable])

	# write file headers
	writeDlmLine!(solFile, "solutions")
	writeDlmLine!(logFile, "log")
	writeDlmLine!(solFile, solFileHeader..., ["item_$i stationIndex" for i = 1:sim.numAmbs]...)
	writeDlmLine!(logFile, logFileHeader..., ["item_$i stationIndex" for i = 1:sim.numAmbs]...)
	flush(solFile)
	flush(logFile)

	for i = 1:numSearches
		logFileDict["search"] = i

		# initial nested compliance table
		nestedCompTable = nestedCompTables[i]
		doPrint && println("Initial nested compliance table $i (of $numSearches):")
		doPrint && printNestedCompTable(nestedCompTable)

		# perform local search until no improvements can be made, then increase max number of reps allowed
		doPrint && println()
		global minNumReps, maxNumRepsList
		numReps = minNumReps # init
		resetNestedCompTablesPeriodStatsList!()
		for j = 1:length(maxNumRepsList)
			maxNumReps = maxNumRepsList[j] < sim.numReps ? Int(maxNumRepsList[j]) : sim.numReps
			logFileDict["maxNumRepsListIndex"] = j
			logFileDict["maxNumReps"] = maxNumReps

			doPrint && println("Search $i (of $numSearches)")
			doPrint && println("Setting maximum number of reps to ", maxNumReps)
			nestedCompTable, numReps = localSearch!(sim, nestedCompTable, numReps, maxNumReps, logFile)
			if maxNumReps == sim.numReps break end # increasing max number of reps will not help
		end
		doPrint && println("Finished local search $i (of $numSearches); solution:")
		doPrint && printNestedCompTable(nestedCompTable)

		# save best nested compliance table found for this search
		global nestedCompTableSols
		push!(nestedCompTableSols, nestedCompTable)

		# write solution to file
		objVal, objValHalfWidth = calcMeanAndHalfWidth(objValsLookup(nestedCompTable, [1:numReps;]))
		solFileDict = Dict("search" => i, "numReps" => numReps, "objVal" => objVal, "objValHalfWidth" => objValHalfWidth)
		merge!(solFileDict, getPeriodsStatsDict(periodStatsListLookup(nestedCompTable)[1:numReps]))
		fileWriteDlmLine!(solFile, solFileHeader, solFileDict, nestedCompTable)

		# write solutions as priority list file
		global nestedCompTablesOutputFilename
		writePriorityListsFile(nestedCompTablesOutputFilename, nestedCompTableSols)
	end

	# print out all results from each finished local search
	doPrint && println()
	doPrint && println("Nested compliance tables from completed local searches:")
	for i = 1:numSearches
		doPrint && println()
		doPrint && println("Search: ", i)
		nestedCompTable = nestedCompTableSols[i]
		doPrint && printNestedCompTable(nestedCompTable)
		doPrint && println("objective value = ", mean(objValsLookup(nestedCompTable)))
	end

	close(solFile)
	close(logFile)
end

# Local search of nested compliance table, starting at nestedCompTable.
# Move items in nested compliance table, making changes that optimise the
# objective function (objFn) value for the given sense (:min or :max).
# Continue search until no improvement can be made.
# Mutates: sim, logFile
function localSearch!(sim::Simulation, nestedCompTable::NestedCompTable, nestedCompTableNumReps::Int, maxNumReps::Int, logFile::IOStream)::Tuple{NestedCompTable,Int}

	global minNumReps, maxNumRepsIncrease, logFileHeader, logFileDict, sense, doPrint

	# shorthand
	numAmbs = sim.numAmbs
	numStations = sim.numStations

	# track iterations
	iter = 1 # count iterations performed
	startTime = time()
	getSearchDuration() = round(time() - startTime, digits = 2)

	# track nested compliance tables tried
	nestedCompTablesTried = Set{NestedCompTable}()

	# keep track of how many reps were simulated for each nested compliance table
	nestedCompTablesRepsDone = Dict{NestedCompTable, Int}()
	getNestedCompTableRepsDone!(nestedCompTable::NestedCompTable) = get!(nestedCompTablesRepsDone, nestedCompTable, 0)
	updateNestedCompTablesRepsDone!(nestedCompTable::NestedCompTable, n::Int) =
		(nestedCompTablesRepsDone[nestedCompTable] = max(n, getNestedCompTableRepsDone!(nestedCompTable)))

	# calculate objective value for starting point
	numReps = nestedCompTableNumReps # do not just use minNumReps, as nestedCompTable may have already been simulated a number of reps
	objVals, periods, numLookups = simObjVals!(sim, nestedCompTable, 1:numReps)
	(objVal, objValHalfWidth) = calcMeanAndHalfWidth(objVals)
	updateNestedCompTablesRepsDone!(nestedCompTable, numReps)
	push!(nestedCompTablesTried, nestedCompTable)
	doPrint && println("starting objective value: ", objVal)

	# starting solution is current best
	bestNestedCompTable = copy(nestedCompTable)
	bestObjVals = objVals
	bestObjVal = objVal

	# write starting point
	for (key, value) in logFileDict
		if !in(key, ("search", "maxNumRepsListIndex", "maxNumReps"))
			logFileDict[key] = "" # reset value
		end
	end
	merge!(logFileDict, getPeriodsStatsDict(periods))
	merge!(logFileDict, Dict("iter" => iter, "numReps" => numReps, "numLookups" => numLookups, "searchDurationSeconds" => getSearchDuration(),
		"objVal" => objVal, "objValHalfWidth" => objValHalfWidth, "bestObjVal" => bestObjVal, "objValDiff" => 0, "objValDiffHalfWidth" => 0))
	fileWriteDlmLine!(logFile, logFileHeader, logFileDict, nestedCompTable)

	# local search
	doPrint && println("starting local search")
	endPoint = nothing
	iter += 1
	while true
		for i = (numAmbs+numStations):-1:1, j = 1:numAmbs, move = ["insert", "swap"]
			# starting with i > numAmbs will try insert/swap new items into nested compliance table first (instead of only moving existing items within list); should help with speed

			# change nested compliance table, using re-insertion / swapping
			# if i <= numAmbs, move item in nestedCompTable[i] to nestedCompTable[j]
			# if i > numAmbs, try move slot from station i-numAmbs (if num spare slots > 0)
			nestedCompTable = copy(bestNestedCompTable)
			if move == "insert"
				insert!(nestedCompTable, i, j)
			elseif move == "swap"
				swap!(nestedCompTable, i, j)
			end

			newBestFound = false
			if !in(nestedCompTable, nestedCompTablesTried)
				doPrint && print(move, ": i = ", i, ", j = ", j)

				# simulate until there is a significant difference in objective value between nestedCompTable and bestNestedCompTable,
				# or until running out of replications to simulate
				numReps = minNumReps
				meanDiff = halfWidth = 0.0 # init
				objVals = ObjVal[] # for nestedCompTable
				periods = SimPeriodStats[] # for nestedCompTable
				numLookups = 0
				while numReps <= maxNumReps
					newObjVals, newPeriods, newNumLookups = simObjVals!(sim, nestedCompTable, (length(objVals)+1):numReps)
					push!(nestedCompTablesTried, nestedCompTable)
					push!(objVals, newObjVals...)
					push!(periods, newPeriods...)
					numLookups += newNumLookups
					updateNestedCompTablesRepsDone!(nestedCompTable, numReps)
					if length(bestObjVals) < length(objVals)
						push!(bestObjVals, simObjVals!(sim, bestNestedCompTable, (length(bestObjVals)+1):numReps)[1]...)
						updateNestedCompTablesRepsDone!(bestNestedCompTable, numReps)
					end
					meanDiff, halfWidth = calcMeanAndHalfWidth(objVals[1:numReps] - bestObjVals[1:numReps])
					if numReps == maxNumReps
						break # ran out of reps, will have to guess which nested compliance table is better based on mean
					elseif (sense == :max && meanDiff < -halfWidth) || (sense == :min && meanDiff > halfWidth)
						break # bestNestedCompTable is better, stop
					elseif (sense == :max && meanDiff > halfWidth) || (sense == :min && meanDiff < -halfWidth)
						# nestedCompTable may be better
						# need to simulate nestedCompTable up to much as bestNestedCompTable, to avoid possibility of infinite loop
						newNumReps = getNestedCompTableRepsDone!(bestNestedCompTable)
						if numReps < newNumReps
							numReps = min(maxNumReps, newNumReps)
						elseif numReps == newNumReps
							break # done, no need to simulate nestedCompTable further
						else
							error("numReps > newNumReps, this should not happen.")
						end
					else # not sure which nested compliance table is better, simulate more reps
						newNumReps = numReps * (halfWidth / meanDiff) ^ 2 # since standard error is roughly proportional to 1/sqrt(n) for n samples
						if isnan(newNumReps) newNumReps = maxNumReps end # meanDiff was zero, just use all reps
						numReps = ceil(Int, max(numReps + 1, min(newNumReps, numReps + maxNumRepsIncrease, maxNumReps))) # bound newNumReps
					end
				end

				objVal, objValHalfWidth = calcMeanAndHalfWidth(objVals)
				bestObjVal = mean(bestObjVals)
				# note that length(objVals) <= length(bestObjVals), as the nested compliance tables may have been simulated different amounts (reps)
				doPrint && println("; objective value: ", objVal, " (best: ", bestObjVal, ")")

				usedMove = false
				if (sense == :max && meanDiff > 0) || (sense == :min && meanDiff < 0) # do not compare objVals and bestObjVals (or their means) as lengths may differ
					# keep change
					bestNestedCompTable = nestedCompTable
					bestObjVals = objVals
					bestObjVal = objVal
					usedMove = true
					newBestFound = true

					resetNestedCompTablesPeriodStatsList!() # remove periods to save on memory; new solution has almost completely new neighbourhood
					setNestedCompTablePeriodStatsList!(nestedCompTable, 1:length(periods), periods)

					doPrint && println("made ", move, ": i = ", i, ", j = ", j, "; new objective value: ", bestObjVal)
				end

				merge!(logFileDict, getPeriodsStatsDict(periods))
				merge!(logFileDict, Dict("iter" => iter, "move" => move, "i" => i, "j" => j, "numReps" => numReps, "numLookups" => numLookups, "usedMove" => Int(usedMove),
					"searchDurationSeconds" => getSearchDuration(), "objVal" => objVal, "objValHalfWidth" => objValHalfWidth, "bestObjVal" => bestObjVal,
					"objValDiff" => meanDiff, "objValDiffHalfWidth" => halfWidth))
				fileWriteDlmLine!(logFile, logFileHeader, logFileDict, nestedCompTable)

				iter += 1
			end

			if newBestFound || endPoint == nothing
				endPoint = (i,j,move)
			elseif endPoint == (i,j,move)
				doPrint && println("Found local optimum.")
				return bestNestedCompTable, getNestedCompTableRepsDone!(bestNestedCompTable)
			end
		end
	end
end

# for the given nested compliance table, count how many slots are still available for the station
function countStationSpareSlots(nestedCompTable::NestedCompTable, stationIndex::Int)
	global stationCapacities
	stationSpareSlots = stationCapacities[stationIndex]
	for i in nestedCompTable
		if i == stationIndex
			stationSpareSlots -= 1
		end
	end
	return stationSpareSlots
end

# Insert item at index i into index j in nested compliance table.
# If i > numAmbs, insertion is done with station index i-numAmbs.
# Return true if insertion completed, false otherwise.
function insert!(nestedCompTable::NestedCompTable, i::Int, j::Int)
	if i == j return false end

	# shorthand
	numAmbs = length(nestedCompTable)
	numStations = length(stationCapacities) # global stationCapacities

	if i <= numAmbs
		@assert(i >= 1)
		# re-insert item in nestedCompTable[i] to nestedCompTable[j]
		temp = nestedCompTable[i]
		if i < j
			for k = i:j-1
				nestedCompTable[k] = nestedCompTable[k+1]
			end
		elseif i > j
			for k = i:-1:j+1
				nestedCompTable[k] = nestedCompTable[k-1]
			end
		else
			error()
		end
		nestedCompTable[j] = temp
	else
		i -= numAmbs # so now i = station index
		@assert(1 <= i && i <= numStations)
		stationSpareSlots = countStationSpareSlots(nestedCompTable, i)
		if stationSpareSlots == 0
			return false
		end
		# insert slot from station i to nestedCompTable[j] (if stationSpareSlots > 0)
		temp = nestedCompTable[numAmbs] # will be bumped from end of nestedCompTable
		for k = numAmbs:-1:j+1
			nestedCompTable[k] = nestedCompTable[k-1]
		end
		nestedCompTable[j] = i
	end

	return true
end

# Swap items at index i and index j in nested compliance table.
# If i > numAmbs, swap is done with station index i-numAmbs.
# Return true if swap completed, false otherwise.
function swap!(nestedCompTable::NestedCompTable, i::Int, j::Int)
	if i == j return false end

	# shorthand
	numAmbs = length(nestedCompTable)
	numStations = length(stationCapacities) # global stationCapacities

	@assert(1 <= i && i <= numAmbs + numStations)
	@assert(1 <= j && j <= numAmbs)

	# perform swap if stations are different
	if i <= numAmbs
		if nestedCompTable[i] == nestedCompTable[j]
			return false
		end
		(nestedCompTable[i], nestedCompTable[j]) = (nestedCompTable[j], nestedCompTable[i])
	else # elseif i <= numAmbs + numStations
		i -= numAmbs # so now i = station index
		@assert(1 <= i && i <= numStations)
		stationSpareSlots = countStationSpareSlots(nestedCompTable, i)
		if i == nestedCompTable[j] || stationSpareSlots == 0
			return false
		end
		# swap i with nestedCompTable[j]
		nestedCompTable[j] = i
	end

	return true
end

# mutates: sim
function setSimStatsCapture!(sim::Simulation, periodDurationsIter::Stateful, warmUpDuration::Float)
	stats = sim.stats = deepcopy(sim.backup.stats)
	stats.doCapture = true
	(stats.periodDurationsIter, stats.warmUpDuration) = (periodDurationsIter, warmUpDuration)
	stats.nextCaptureTime = sim.startTime + (stats.warmUpDuration > 0 ? stats.warmUpDuration : first(stats.periodDurationsIter))
	sim.backup.stats = deepcopy(sim.stats)
end

# parameters
configFilename = "sim_config.xml"
outputFolder = "output"
numSearches = 1
nestedCompTableRng = Random.MersenneTwister(0) # useful for reproducing results, if using random restarts
# replication parameters
warmUpDuration = 0.1
periodDuration = 1.0

# parameter checks
@assert(numSearches >= 1)
@assert(warmUpDuration >= 0)
@assert(periodDuration > 0)

# initialise sim
t = time()
@assert(isfile(configFilename))
println("Initialising simulation from config: ", configFilename)
sim = initSim(configFilename, createBackup = true)
@assert(sim.numReps >= 2, "Need at least two simulation replications (sim.reps) for statistics.")

# set sim stats capturing
setSimStatsCapture!(sim, Stateful(periodDuration), warmUpDuration)
periodEndTime = sim.startTime + warmUpDuration + periodDuration
for rep in sim.reps
	# statistics should stop collecting before last call arrives (if not earlier), to avoid cool-down period
	@assert(periodEndTime <= rep.calls[end].arrivalTime)
end

# generate nested compliance tables
stationCapacities = [station.capacity for station in sim.stations]
nestedCompTables = [makeRandNestedCompTable(sim.numAmbs, sim.numStations; stationCapacities = stationCapacities, rng = nestedCompTableRng) for i = 1:numSearches]

# set sim to use nested compliance table for move up
setMoveUpModule!(sim, compTableModule)

makeRepsRunnable!(sim) # copy changes to sim (stats, nested compliance table) to sim.reps

# run
println("Starting local search(es).")
nestedCompTableLocalSearch!(sim, nestedCompTables, outputFolder = outputFolder)
println("Total runtime: ", round(time()-t, digits = 2), " seconds")


text:
"""
Obtain structural properties of a carbon-fiber, tubular spar,
given internal and external diameters, local chord and
displacement from quarter chord
"""
function structural_tube(
    Ri::Real,
    Ro::Real,
    c::Real = 1.0;
    e::Real = 0.0,
    E::Real = 57.9e9, # AeroDesign ITA: carbono em tecido
    G::Real = 2.710e9,
    ρ::Real = 1.8e3
)

    R = (Ri + Ro) / 2
    t = (Ro - Ri)

    A = π * R ^ 2

    GJ = 4 * A ^ 2 * t * G / (2 * π * R)

    A = π * (Ro ^ 2 - Ri ^ 2)

    EA = E * A

    I = π * (Ro ^ 4 - Ri ^ 4) / 4

    EI = E * I

    m = A * ρ

    Ixx = ρ * π * (Ro ^ 4 - Ri ^ 4) / 2 + c ^ 2 * e ^ 2 * m

    mx = e * m

    EA, EI, GJ, m, Ixx, mx

end

"""
Get carbon-fiber fuselage structural data, given
radius and skin thickness
"""
structural_fuselage(
    R::Real,
    t::Real;
    E::Real = 19.3e9, # AeroDesign ITA: carbono em tecido
    G::Real = (1.0 + 0.35) * 19.3e9 / 2,
    ρ::Real = 1.4e3
) = structural_tube(
    R - t, R;
    E = E,
    G = G,
    ρ = ρ
)

"""
Get approximate structural properties of a balsa wood structural skin,
given internal airfoil area (for unit chord) and local chord
"""
function structural_skin(
    Arel::Real,
    t::Real,
    c::Real = 1.0;
    G::Real = 3.4e9 / (2 * (1.0 + 0.009)),
    ρ::Real = 0.16e3
)

    GJ = 4 * Arel ^ 2 * t * G * c / 2

    m = t * 2 * ρ * c
    Ixx = m * c/ 12 + 0.25 ^ 2 * m * c

    EA = EI = 0.0

    mx = 0.25 * m

    EA, EI, GJ, m, Ixx, mx

end

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
"""
Kind of Petri Net.
Petri Net Type Definition (pntd) URI mapped to PnmlType subtype singleton.

# Imports
$(DocStringExtensions.IMPORTS)

# Exports
$(DocStringExtensions.EXPORTS)
"""
module PnmlTypes
using PNML
using DocStringExtensions

export PnmlType,
        PnmlCore, PTNet,
        AbstractHLCore, HLCore, PT_HLPNG, SymmetricNet, HLNet,
                        StochasticNet, TimedNet, OpenNet,
        AbstractContinuousCore, ContinuousNet
export pnmltype, pntd_symbol

"""
Abstract root of a dispatch type based on Petri Net Type Definition (pntd).

Each Petri Net Markup Language (PNML) network element will have a single pntd URI
as a required 'type' XML attribute. That URI should refer to a RelaxNG schema defining
the syntax and semantics of the XML model.

Selected abbreviations, URIs that do not resolve to a valid schema file, are also allowed.

Refer to [`pntd_symbol`](@ref) and [`pnmltype`](@ref)
for how to get from the URI to a singleton.

$(TYPEDEF)
"""
abstract type PnmlType end

"""
The most minimal concrete Petri Net.
Used to implement and test the complete Petri Net Graph infrastructure.
Labels of the graph is where meaning is attached.
Much of the Label infrastructure for High Level Petri Net Graphs is tested at this level.
Subtypes of `PnmlType` should be used to specialize Labels for expressiveness and optimization.

$(TYPEDEF)
"""
struct PnmlCore <: PnmlType end

"""
Place-Transition Petri Nets add small extensions to core.
The grammer file is ptnet.pnml so we name it PTNet.
Note that 'PT' is often the prefix for XML tags specilized for this net type.

$(TYPEDEF)
"""
struct PTNet <: PnmlType end

"""
Base of High Level Petri Net pntds.
See [`SymmetricNet`](@ref), [`PT_HLPNG`](@ref) and others.

$(TYPEDEF)
"""
abstract type AbstractHLCore <: PnmlType end

"""
High-Level Petri Net Graphs (HLPNGs) add large extensions to core.
The `HLCore` can be used for generic high-level nets.
We try to implement and test all function at `PnmlCore` level,
but expect to find use for a concrete type at this level.

$(TYPEDEF)
"""
struct HLCore <: AbstractHLCore end

"""
HLNet is the most intricate High-Level Petri Net schema.
It extends [`SymmetricNet`](@ref), including with
   - declarations for sorts and functions (ArbitraryDeclarations)
   - sorts for Integer, String, and List

$(TYPEDEF)
"""
struct HLNet <: AbstractHLCore end

"""
Place-Transition Net in HLCore notation (HLPNG=High-Level Petri Net Graph).

$(TYPEDEF)
"""
struct PT_HLPNG <: AbstractHLCore end

"""
Symmetric Petri Net is the best-worked use case on the `primer`
and ISO specification part 2.

$(TYPEDEF)
"""
struct SymmetricNet <: AbstractHLCore end

"""
$(TYPEDEF)

Uses floating point numbers for markings, inscriptions, and conditions.
"""
abstract type AbstractContinuousCore <: PnmlType end

"""
TODO: Continuous Petri Net
Concrete type.
$(TYPEDEF)
"""
struct ContinuousNet <: AbstractContinuousCore end

"""
TODO: Stochastic Petri Net

$(TYPEDEF)
"""
struct StochasticNet <: AbstractHLCore end

"""
TODO: Timed Petri Net

$(TYPEDEF)
"""
struct TimedNet <: AbstractHLCore end

"""
TODO: Open Petri Net

$(TYPEDEF)
"""
struct OpenNet <: AbstractHLCore end

"""
$(TYPEDEF)

Map from Petri Net Type Definition (pntd) URI to Symbol.
Allows multiple strings to map to the same pntd.

There is a companion map [`pnmltype_map`](@ref) that takes the symbol to a type object.

The URI is a string and may be the full URL of a pntd schema,
just the schema file name, or a placeholder for a future schema.

# Examples

The 'pntd symbol' should match the name used in the URI with inconvinient characters
removed or replaced. For example, '-' is replaced by '_'.
"""
const default_pntd_map = Dict{String, Symbol}(
    "http://www.pnml.org/version-2009/grammar/ptnet" => :ptnet,
    "http://www.pnml.org/version-2009/grammar/highlevelnet" => :hlnet,
    "http://www.pnml.org/version-2009/grammar/pnmlcoremodel" => :pnmlcore,
    "http://www.pnml.org/version-2009/grammar/pnmlcore" => :pnmlcore,
    "http://www.pnml.org/version-2009/grammar/pt-hlpng" => :pt_hlpng,
    "http://www.pnml.org/version-2009/grammar/symmetricnet" => :symmetric,
    "pnmlcore"   => :pnmlcore,
    "ptnet"      => :ptnet,
    "hlnet"      => :hlnet,
    "hlcore"     => :hlcore,
    "pt-hlpng"   => :pt_hlpng,
    "pt_hlpng"   => :pt_hlpng,
    "symmetric"  => :symmetric,
    "symmetricnet" => :symmetric,
    "stochastic"   => :stochastic,
    "timed"        => :timednet,
    "nonstandard"  => :pnmlcore,
    "open"         => :pnmlcore,
    "continuous"   => :continuous,
    )

"""
$(TYPEDEF)

The key Symbols are the supported kinds of Petri Nets.
"""
const pnmltype_map = Dict{Symbol, PnmlType}(
    :pnmlcore   => PnmlCore(),
    :hlcore     => HLCore(),
    :ptnet      => PTNet(),
    :hlnet      => HLNet(), #HLPNG?
    :pt_hlpng   => PT_HLPNG(),
    :symmetric  => SymmetricNet(),
    :stochastic => StochasticNet(),
    :timednet   => TimedNet(),
    :continuous => ContinuousNet(),
    )

"""
$(TYPEDSIGNATURES)

Add or replace mapping from symbol `s` to nettype dispatch singleton `t`.
"""
add_nettype!(dict::AbstractDict, s::Symbol, pntd::T) where {T<:PnmlType} =
    dict[s] = pntd #TODO test this

"""
$(TYPEDSIGNATURES)

Map string `s` to a pntd symbol using [`default_pntd_map`](@ref).
Any unknown `s` is mapped to `:pnmlcore`.
Returned symbol is suitable for [`pnmltype`](@ref) to use to index into [`pnmltype_map`](@ref).

# Examples

```jldoctest
julia> using PNML: PnmlTypes.pntd_symbol

julia> pntd_symbol("foo")
:pnmlcore
```
"""
pntd_symbol(s::String) = get(default_pntd_map::Dict{String,Symbol}, s, :pnmlcore)

"""
    pnmltype(pntd::T; kw...)
    pnmltype(uri::AbstractString; kw...)
    function pnmltype(s::Symbol; pnmltype_map=pnmltype_map, kw...)

Map either a text string or a symbol to a dispatch type singlton.

While that string may be a URI for a pntd, we treat it as a simple string without parsing.
The [`PnmlTypes.pnmltype_map`](@ref) and [`PnmlTypes.default_pntd_map`](@ref) are both assumed to be correct here.

Unknown or empty `uri` will map to symbol `:pnmlcore`.
Unknown `symbol` throws a [`PNML.MalformedException`](@ref)

# Examples

```jldoctest
julia> using PNML; using PNML.PnmlTypes: pnmltype, pntd_symbol

julia> pnmltype(PnmlCore())
PnmlCore()

julia> pnmltype("nonstandard")
PnmlCore()

julia> pnmltype(:symmetric)
SymmetricNet()
```
"""
function pnmltype end
pnmltype(pntd::T; kw...) where {T<:PnmlType} = pntd
pnmltype(uri::AbstractString; kw...) = pnmltype(pntd_symbol(uri))
function pnmltype(s::Symbol)
    typemap = pnmltype_map::Dict{Symbol, PnmlType}
    !haskey(typemap, s) && throw(DomainError("Unknown PNTD symbol $s"))
    @inbounds typemap[s]
end

#TODO add traits

end # module PnmlTypes


text:
# Computes barycentric weights, can be used for value interpolation
# from finite element mesh

#   defining quadratic function
quadratic(a, sqr_term, b) = (-b + sqr_term) / 2a

#   finding positive and negative components of quadratic function
function quadratic2(a::T, b::T, c::T) where T<:Real
    sqr_term = sqrt(b^2-4a*c)
    r1 = quadratic(a, sqr_term, b)
    r2 = quadratic(a, -sqr_term, b)

    r1, r2  # the two roots
end

function computeBaryWeights!(w::Vector{T},x::Vector{T}, y::Vector{T}, p::Vector{T}) where T<:Real
    #   Computing the Barycentric weights
    #

    #   INPUTS:
    #   x  = (x1, x2, x3, x4) [vector of x values]
    #   y  = (y1, y2, y3, y4) [vector of y values]
    #   z  = (f(x1, y1), f(x2, y2), f(x3, y3), f(x4, y4)) [vector of z values]
    #   xi = x-value of point we want to interpolate
    #   yi = y-value of point we want to interpolate

    #   OUTPUT (IN PLACE):
    #   w  = (w1, w2, w3, w4) [vector of barycentric weights for sanity check]

    # dummy variable
    bad = -1000.0

    a = [x[1] - p[1], y[1]-p[2], 0.0]
    b = [x[2] - x[1], y[2]-y[1], 0.0]
    c = [x[4] - x[1], y[4]-y[1], 0.0]
    d = [x[1] - x[2] - x[4] + x[3], y[1] - y[2] - y[4] + y[3], 0.0]

    # FINDING μ terms
    A = cross(c,d)
    B = cross(c,b)+cross(a,d)
    C = cross(a,b)

    if (abs(A[3]) < 1.0e-10) # if A=0, we have a line
        u1 = -C[3]/B[3]
        u2 = u1
    else
        sqr_term = sqrt(B[3]^2-4*A[3]*C[3])
        if (sqr_term > zero(T)) # only looking at real solutions
            u1, u2 = quadratic2(A[3], B[3], C[3])
        else #complex solutions don't matter
            u1 = bad
            u2 = u1
        end
    end

    μ = bad
    if (u1 >= zero(T) && u1 <= one(T))
        μ = u1
    end
    if (u2 >= zero(T) && u2 <= one(T))
        μ = u2
    end

    # FINDING λ terms
    A = cross(b,d)
    B = cross(b,c)+cross(a,d)
    C = cross(a,c)

    if (abs(A[3]) < 1.0e-10) # if A=0, we have a line
        u3 = -C[3]/B[3]
        u4 = u3
    else
        sqr_term = sqrt(B[3]^2-4*A[3]*C[3])
        if (sqr_term > zero(T)) # only looking at real solutions
            u3, u4 = quadratic2(A[3], B[3], C[3])
        else #complex solutions don't matter
            u3 = bad
            u4 = u3
        end
    end

    λ = bad
    if (u3 >= zero(T) && u3 <= one(T))
        λ = u3
    end
    if (u4 >= zero(T) && u4 <= one(T))
        λ = u4
    end

    # barycentric weights
    w[1] = (1-μ)*(1-λ)
    w[2] = λ*(1-μ)
    w[3] = λ*μ
    w[4] = (1-λ)*μ

    return w #returning vector of weights
end
# w, a, b, c input just to avoid re-allocation
function computeBaryWeights!(w::Vector{T},polygon::Vector{Point2D{T}}, p::Point2D{T},
                                a::Vector{T},b::Vector{T},c::Vector{T},d::Vector{T}) where T<:Real
    #   Computing the Barycentric weights
    #

    #   INPUTS:
    #   x  = (x1, x2, x3, x4) [vector of x values]
    #   y  = (y1, y2, y3, y4) [vector of y values]
    #   z  = (f(x1, y1), f(x2, y2), f(x3, y3), f(x4, y4)) [vector of z values]
    #   xi = x-value of point we want to interpolate
    #   yi = y-value of point we want to interpolate

    #   OUTPUT (IN PLACE):
    #   w  = (w1, w2, w3, w4) [vector of barycentric weights for sanity check]

    # dummy variable
    bad = -1000.0

    a[1] = polygon[1].x - p.x
    a[2] = polygon[1].y-p.y
    a[3] = 0.0
    b[1] = polygon[2].x - polygon[1].x
    b[2] = polygon[2].y-polygon[1].y
    b[3] = 0.0
    c[1] = polygon[4].x - polygon[1].x
    c[2] = polygon[4].y-polygon[1].y
    c[3] = 0.0
    d[1] = polygon[1].x - polygon[2].x - polygon[4].x + polygon[3].x
    d[2] = polygon[1].y - polygon[2].y - polygon[4].y + polygon[3].y
    d[3] = 0.0

    # FINDING μ terms
    A1 = c[2]*d[3]-c[3]*d[2]; B1 = c[2]*b[3]-c[3]*b[2] + a[2]*d[3]-a[3]*d[2]; C1 = a[2]*b[3]-a[3]*b[2]
    A2 = c[3]*d[1]-c[1]*d[3]; B2 = c[3]*b[1]-c[1]*b[3] + a[3]*d[1]-a[1]*d[3]; C2 = a[3]*b[1]-a[1]*b[3]
    A3 = c[1]*d[2]-c[2]*d[1]; B3 = c[1]*b[2]-c[2]*b[1] + a[1]*d[2]-a[2]*d[1]; C3 = a[1]*b[2]-a[2]*b[1]

    #A = cross(c,d)
    #B = cross(c,b) + cross(a,d)
    #C = cross(a,b)

    #if (abs(A[3]) < 1.0e-10) # if A=0, we have a line
    if (abs(A3) < 1.0e-10) # if A=0, we have a line
        #u1 = -C[3]/B[3]
        u1 = -C3/B3
        u2 = u1
    else
        #sqr_term = sqrt(B[3]^2-4*A[3]*C[3])
        sqr_term = sqrt(B3^2-4*A3*C3)
        if (sqr_term > zero(T)) # only looking at real solutions
            #u1, u2 = quadratic2(A[3], B[3], C[3])
            u1, u2 = quadratic2(A3, B3, C3)
        else #complex solutions don't matter
            u1 = bad
            u2 = u1
        end
    end

    μ = bad
    if (u1 >= zero(T) && u1 <= one(T))
        μ = u1
    end
    if (u2 >= zero(T) && u2 <= one(T))
        μ = u2
    end

    # FINDING λ terms
    A1 = b[2]*d[3]-b[3]*d[2]; B1 = b[2]*c[3]-b[3]*c[2] + a[2]*d[3]-a[3]*d[2]; C1 = a[2]*c[3]-a[3]*c[2]
    A2 = b[3]*d[1]-b[1]*d[3]; B2 = b[3]*c[1]-b[1]*c[3] + a[3]*d[1]-a[1]*d[3]; C2 = a[3]*c[1]-a[1]*c[3]
    A3 = b[1]*d[2]-b[2]*d[1]; B3 = b[1]*c[2]-b[2]*c[1] + a[1]*d[2]-a[2]*d[1]; C3 = a[1]*c[2]-a[2]*c[1]

    #A = cross(b,d)
    #B = cross(b,c) + cross(a,d)
    #C = cross(a,c)

    #if (abs(A[3]) < 1.0e-10) # if A=0, we have a line
    if (abs(A3) < 1.0e-10) # if A=0, we have a line
        u3 = -C3/B3
        u4 = u3
    else
        #sqr_term = sqrt(B[3]^2-4*A[3]*C[3])
        sqr_term = sqrt(B3^2-4*A3*C3)
        if (sqr_term > zero(T)) # only looking at real solutions
            #u3, u4 = quadratic2(A[3], B[3], C[3])
            u3, u4 = quadratic2(A3, B3, C3)
        else #complex solutions don't matter
            u3 = bad
            u4 = u3
        end
    end

    λ = bad
    if (u3 >= zero(T) && u3 <= one(T))
        λ = u3
    end
    if (u4 >= zero(T) && u4 <= one(T))
        λ = u4
    end

    # barycentric weights
    w[1] = (1-μ)*(1-λ)
    w[2] = λ*(1-μ)
    w[3] = λ*μ
    w[4] = (1-λ)*μ

    return w #returning vector of weights
end
function computeBaryWeights(x::Vector{T}, y::Vector{T}, p::Vector{T}) where T<:Real
    w = zeros(T,4)
    computeBaryWeights!(w::Vector{T},x::Vector{T}, y::Vector{T}, p::Vector{T})

    return w
end

max_stars_repo_name: pseastham/eFEMpart


text:
# You know it's getting real when this is uncommented
# __precompile__()

module BinaryBuilder
using Compat
using Reexport
@reexport using BinaryProvider

include("Auditor.jl")
include("DockerRunner.jl")
include("Dependency.jl")

end # module


text:
#One action, zero-then-one observation, two timesteps long, +1 reward at the end: If my agent
#can learn the value in (2.) but not this one, it must be that my reward discounting is broken.
import AlphaZero.GI
using Flux: onehot
struct GameSpec <: GI.AbstractGameSpec end

mutable struct GameEnv <: GI.AbstractGameEnv
  terminated :: Bool
  numero :: Int
end

const INITSTATE = (; terminated=false, numero=1)

GI.init(::GameSpec, state=INITSTATE) = GameEnv(state.terminated, state.numero)

GI.spec(::GameEnv) = GameSpec()

GI.two_players(::GameSpec) = false

const ACTIONS = collect(1:1)

GI.actions(::GameSpec) = ACTIONS

GI.actions_mask(g::GameEnv) = [true]

GI.current_state(g::GameEnv) = (; terminated=g.terminated, numero=g.numero)

GI.white_playing(g::GameEnv) = true

GI.game_terminated(g::GameEnv) = g.terminated

GI.white_reward(g::GameEnv) = g.terminated ? 1. : 0.

function GI.play!(g::GameEnv, action)
  if g.numero >= 2
    g.terminated = true
  end
  g.numero += 1
end

GI.vectorize_state(::GameSpec, state) = collect(onehot(state.numero, 1:2))


text:
@testset "laguerre" begin
    @test all(
        laguerrel(n, α, x) ≈ laguerre(n, α)(x)
        for n in 0:100, α in 0:100, x in -20:0.1:20
    )
end


text:
module AffineSpaces

using FixedSizeArrays

export AffineSpace, FixedAffineSpace, AffinePoint

"""
    abstract AffineSpace{T}

An `AffineSpace` object represents points in an affine space, which is similar
to a Cartesian space but where no special meaning can be assigned to the origin.
Subtypes of `AffineSpace` should have a field named `pos` with vector-like
abilities describing the coordinates.

The consequences are that points do not have vector-like abilities such as
multiplication by a scalar (scaling the distance relative to the origin), unary
inversion (reflection about the origin) or addition of the coordinates of two
points (coorinates being relative to the origin).

Outside these restriction, one can take affine combinations of points, finding
for instance the midpoint of two points. Affine combinations are built with
the colon operation `p₁:p₂` and indexed with a scalar `(p₁:p₂)[i]`, where
`i = 0.0` indicates the point `p₁`, and `i = 1.0` indicates the point `p₂`, and
other values of `i` trace out the line connecting the two points. Generally,
the affine combination of `n` points can be generated with `(p₁:...:pₙ)[i₁,...,iₙ]`
so long as the `iⱼ`s sum to `1` (or the last `iₙ` is not given, and is
automatically calculated).

Furthermore, standard vectors can be added and subtracted from affine points as
a displacement of the point. The displacement vector between two points is given
by their subtraction.

The goal of these types are to provide some programmer safety to avoid
conceptual mistakes regarding manipulation of points, by removing the meaning
of the possibly-arbitrary origin.
"""
abstract AffineSpace{T}
abstract FixedAffineSpace{N,T} <: AffineSpace

"""
    AffinePoint{N,T} <: FixedAffineSpace{T}

A concrete type describing a point in an affine space, where the Cartesian
coordinates of the point are stored in a `FixedSizeArrays.Vec{N,T}`. See
`AffineSpace` for a description of the interface for affine points.
"""
immutable AffinePoint{N,T} <: FixedAffineSpace{T}
    pos::Vec{N,T}
end

Base.length(point::AffineSpace) = length(point.pos)
Base.size(point::AffineSpace) = size(point.pos)
Base.size(point::AffineSpace, d) = size(point.pos, d)
Base.endof(point::AffineSpace) = endof(point.pos)
Base.getindex(point::AffineSpace, i) = point.pos[i]
Base.setindex!(point::AffineSpace, v, i) = setindex!(point.pos, v, i)
Base.start(point::AffineSpace) = start(point.pos)
Base.next(point::AffineSpace, i) = next(point.pos, i)
Base.done(point::AffineSpace, i) = done(point.pos, i)

function Base.show(io::IO, point::AffinePoint)
    print(io, "AffinePoint(")
    showcompact(io, point.pos)
    print(io,")")
end
function Base.showcompact(io::IO, point::AffinePoint)
    showcompact(io, point.pos)
end

# Adding and subtracting vectors from points
Base.(:+){P <: AffineSpace}(point::P, vec::Union{AbstractVector,FixedVector}) = P(point.pos + vec)
Base.(:+){P <: AffineSpace}(vec::Union{AbstractVector,FixedVector}, point::P) = P(vec + point.pos)

Base.(:-){P <: AffineSpace}(point::P, vec::Union{AbstractVector,FixedVector}) = P(point.pos + vec)

# Subtract to give vectors - uses underlying vector structure of the .pos field
Base.(:-)(point1::AffineSpace, point2::AffineSpace) = point1.pos - point2.pos

# Nice messages for disallowed operations
Base.(:+)(point1::AffineSpace, point2::AffineSpace) = error("Cannot add affine points. Consider adding a displacement vector instead.")
Base.(:*)(point::AffineSpace, x::Number) = error("Cannot scale affine points, since scaling is relative to the origin.")
Base.(:*)(x::Number, point::AffineSpace) = error("Cannot scale affine points, since scaling is relative to the origin.")
Base.(:-)(point::AffineSpace) = error("The additive inverse of an affine point is not defined, since affine points cannot be added.")
Base.(:-)(vec::Union{AbstractVector,FixedVector}, point::AffineSpace) = error("Cannot subtract an affine point from a vector.")

# Taking affine combinations of points
immutable AffineSubSpace{N, P<:AffineSpace}
    points::NTuple{N,P}
end

Base.start(ss::AffineSubSpace) = 1
Base.next(ss::AffineSubSpace, n::Int) = (ss.points[n], n+1)
Base.done{N}(ss::AffineSubSpace{N}, n::Int) = n > N

function Base.getindex{N, P}(ss::AffineSubSpace{N, P}, x...)
    lx = length(x)
    if lx == N
        @assert sum(x) ≈ 1 # maybe consider this a boundscheck?
        return P(sum(map(*, x, map(p -> p.pos, ss))))
    elseif lx == N - 1
        x = (x..., 1 - sum(x))
        return P(sum(map(*, x,  map(p -> p.pos, ss))))
    else
        error("Must specify either $N or $(N-1) numbers to define affine combination of $N points - got $lx")
    end
end

# Allow syntax like (point1:point2)[0.5] to get midpoint between point1 and point2
Base.colon(points::AffineSpace...) = AffineSubSpace(points)
Base.colon(ss::AffineSubSpace, points::AffineSpace...) = AffineSubSpace((ss.points..., points...))

# Base.mean as special case of affine combination
Base.mean{N,P <: AffineSpace}(points::NTuple{N,P}) = P(sum(map(p -> p.pos, points)) / N)

end # module