text:
using LinearAlgebra, SparseArrays
using Logging
using Plots

# Suitable for any linear transformation on any type of multi-dimensional arrays.
# Implementation details and advantages:
# * has the option to keep all previous residual vectors for orthogonalizing the newest
#   residual vectors obtained.
# * The re-orthogonalization process is computed in parallel because it uses a matrix of zero that
#   resizes to double the number of columns each time.

mutable struct ConjGrad
    A::Function              # Linear opeartor
    b::AbstractArray         # RHS vector
    tensor_size::Tuple       # The size of the tensor the linear operator is acting on.

    x::AbstractVector        # current guess, started with initial guess.
    r::AbstractVector        # previous computed residual
    rnew::AbstractVector     # Current residual
    d::AbstractVector        # Conjugate Direction

    itr::UInt64              # Iteration count.
    Q::Union{AbstractMatrix, Nothing}        # Re-Orthogonalization basis.
    Q_size::UInt64           # Number of orthogonalization vectors.
    over_write::UInt64       # position for overwrite once storage is limit is reached.

    storage_limit::UInt64    # storage limit for the Q vector.
    reorthogonalize::Bool    # Whether to perform reorthogonalization.
    function ConjGrad(
        A::Function,
        b::AbstractArray,
        x0::Union{AbstractArray,Nothing}=nothing
    )
        this = new()
        x0 = x0===nothing ? b .+ 0.1 : x0
        this.A = A
        this.b = b
        this.tensor_size = size(b)
        this.x = x0
        this.r = ComputeResidualVec(this, this.x)
        this.rnew = similar(this.r)
        this.d = this.r
        this.x = x0
        this.itr = 0

        this.Q = zeros(typeof(this.r[1]), length(this.r), 4)
        this.Q[:, 1] = this.r/norm(this.r)
        this.Q_size = 1
        this.over_write = 0

        this.storage_limit = length(this.r) - 1       # Maximal limit, if not, we have a problem.
        this.reorthogonalize = true
        return this
    end

    function ConjGrad(A::AbstractMatrix, b::AbstractVector, x0::Union{AbstractArray,Nothing}=nothing)
        return ConjGrad((x)-> A*x, b, x0)
    end

end

"""
    Turn on the reorthogonalizations using the residual vectors,
    and add current residual to the list of residuals.
"""
function TurnOnReorthgonalize(this::ConjGrad)
    this.reorthogonalize = true
    this.Q = zeros(typeof(this.r[1]), length(this.r), 4)
    this.Q[:, 1] = this.r/norm(this.r)
    this.Q_size = 1
    this.over_write = 0
    return
end

"""
    Turn off the reorthogonalization on the residual vectors. And then
    clear all the stored residual vectors.
"""
function TurnOffReorthgonalize(this::ConjGrad)
    this.reorthogonalize = false
    this.Q = nothing
    return
end

"""
    Change the storage limit for the number of vectors used for
    re-orthogonalization.
"""
function ChangeStorageLimit(this::ConjGrad, storage_limit)
    this.storage_limit = min(length(this.r) - 1, storage_limit)
    return
end

"""
    Apply the linear transformation on the tensor and return the
    residual vector after transformation.
"""
function ComputeResidualVec(this::ConjGrad, x::AbstractArray)
    resVec = this.b - this.A(reshape(x, this.tensor_size))
    return reshape(resVec, length(resVec))
end

"""
    Compute the vector using the linear transformation.
"""
function ComputeVec(this::ConjGrad, x::AbstractArray)
    vec = this.A(reshape(x, this.tensor_size))
    return reshape(vec, length(vec))
end

"""
    Performs One step of conjugate Gradient.
"""
function (this::ConjGrad)()
    r = this.r
    if norm(r) == 0
        return 0 # The problem is solved already.
    end
    d = this.d
    Ad = ComputeVec(this, d)
    a = dot(r, r)/dot(d, Ad)
    if a < 0
        error("CG got a non-definite matrix")
    end

    this.x += a*d
    this.rnew = r - a*Ad                    # update rnew

    if this.reorthogonalize
        this.rnew -= this.Q*this.Q'*this.rnew

        if this.Q_size == this.storage_limit  # starts overwrite
            this.Q[:, this.over_write + 1] = this.rnew/norm(this.rnew)
            this.over_write = (this.over_write + 1)%this.storage_limit

        elseif this.Q_size == size(this.Q, 2) # Resize
            newQ = zeros(
                typeof(this.r[1]),
                size(this.Q, 1),
                max(2*size(this.Q, 2), this.storage_limit)
            )
            newQ[:, 1:this.Q_size] = this.Q
            newQ[:, this.Q_size + 1] = this.rnew/norm(this.rnew)
            this.Q = newQ
            this.Q_size += 1
        else                                  # add one more.
            this.Q[:, this.Q_size + 1] = this.rnew/norm(this.rnew)
            this.Q_size += 1
        end
    end

    b = dot(this.rnew, this.rnew)/dot(r, r)
    # @assert abs(dot(rnew + β*d, Ad)) < 1e-8 "Not conjugate"
    this.d = this.rnew + b*d
    this.r = this.rnew                      # Override
    this.itr += 1
    return convert(Float64, norm(this.rnew))
end

### Experiments.

function GetErrorEnergyNorm(A, b, xs)
    xStar = A\b
    return xs .|> x -> dot(x - xStar, A, x-xStar)/dot(xs[1] - xStar, A, xs[1] - xStar)
end

"""
    Bad eigen value distribution, they are distributed from 0.001 to 1 like a
    geometric series.

    * A lot of small eigenvalues are clustered close to zero, few larger ones are
    far away from the other Eigenvalues.

    Parameters:
        rho: Number:
            An number between 0 and 1 that paramaterize the distribution of the
            eigenvalues for the PSD matrix.
        N=20:
            The size of the PSD matrix.

    returns:
        a diagonal matrix.
"""
function GetNastyPSDMatrix(rho::Number, N=20, inverted=false)
    @assert rho <= 1 && rho >= 0
    EigenValues = zeros(N)
    EigenMin, EigenMax = 0.001, 1    # Min Max Eigenvalues.
    EigenValues[1] = EigenMin
    for IdexI in 2:N
        EigenValues[IdexI] = EigenMin +
            ((IdexI - 1)/(N - 1))*(EigenMax - EigenMin)*rho^(N - IdexI)  # formulas
    end
    if inverted
        return Diagonal(1 + EigenMin - EigenValues)
    end
    return Diagonal(EigenValues)
end

N = 256
A = GetNastyPSDMatrix(0.9, N)
B = sprand(Float64, N, N, 2/N)
A = A + B'*B
b = rand(N)
cg1 = ConjGrad(A, b)
cg2 = ConjGrad(A, b)
cg3 = ConjGrad(A, b)
cg1Soln = Vector{typeof(b)}()
cg2Soln = Vector{typeof(b)}()
cg3Soln = Vector{typeof(b)}()
push!(cg1Soln, cg1.x)
push!(cg2Soln, cg2.x)
push!(cg3Soln, cg3.x)
TurnOffReorthgonalize(cg2)
ChangeStorageLimit(cg3, 188)
for _ in 1:3*N
    cg1(); cg2(); cg3()
    push!(cg1Soln, cg1.x); push!(cg2Soln, cg2.x); push!(cg3Soln, cg3.x)
end

objectiveVals1 = GetErrorEnergyNorm(A, b, cg1Soln)
objectiveVals2 = GetErrorEnergyNorm(A, b, cg2Soln)
objectiveVals3 = GetErrorEnergyNorm(A, b, cg3Soln)
fig = plot(log10.(objectiveVals1), label="With Re-orthonalization")
plot!(fig, log10.(objectiveVals2), label="without Re-orthogonalization")
plot!(fig, log10.(objectiveVals3), label="partial Re-orthogonalization")


text:
function fib(x, y)
  return y, x+y
end
a, b, s = 1, 1, 0
while a+b < 4000_000
  global a, b, s
  a, b = fib(a, b)
  if iseven(b)
    s += b
  end
end
println(s)

ext: jl

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
struct ArgLooper{SA, SB, RI}  <: AbstractArgIter
    dict_list::Vector{SA}
    runs_iter::RI
    stable_arg::SB
    done::Bool
    run_name::String
end

function ArgLooper(dict_list::Vector{<:Dict}, stable_arg::Dict, runs_iter, run_param)
    ArgLooper(dict_list, runs_iter, stable_arg, false, run_param)
end

function ArgLooper(str_list::Vector{Vector{String}}, stable_arg::Vector{String}, runs_iter, run_param)
    @warn "Arg Looper w/ Arg Parse is depricated in favor of passing dictionaries around."
    ArgLooper(str_list, runs_iter, stable_arg, false, run_param)
end

ArgLooper(list, stable_arg, runs_iter, run_name="--run") =
    ArgLooper(list, runs_iter, stable_arg, false, run_name)

set_save_dir!(iter::ArgLooper, path) = iter.stable_arg["save_dir"] = path

function make_arguments(iter::ArgLooper{SA, RI}, state) where {SA<:Vector{String}, RI}
    arg_list = Vector{String}()
    arg_list = [iter.dict_list[state[2][1]]; [iter.run_name, string(iter.runs_iter[state[2][2]])]; iter.stable_arg]
    return arg_list
end

function make_arguments(iter::ArgLooper{SA, RI}, state) where {SA<:Dict, RI}
    arg_list = Vector{String}()
    arg_list = merge(iter.dict_list[state[2][1]], iter.stable_arg)
    arg_list[iter.run_name] = iter.runs_iter[state[2][2]]
    return arg_list
end

function Base.iterate(iter::ArgLooper)
    state = (1, [1, 1])
    arg_list = make_arguments(iter, state)
    return (state[1], arg_list), next_state(iter, state)
end

function next_state(iter::ArgLooper, _state)
    state = _state[2]
    n_state = _state[1]

    state[2] += 1
    if state[2] > length(iter.runs_iter)
        state[2] = 1
        state[1] += 1
    end
    return (n_state+1, state)

end

function Base.iterate(iter::ArgLooper, state)
    if state[2][1] > length(iter.dict_list)
        return nothing
    end
    arg_list = make_arguments(iter, state)
    return (state[1], arg_list), next_state(iter, state)
end

function Base.length(iter::ArgLooper)
    return length(iter.dict_list)*length(iter.runs_iter)
end


text:
# Autogenerated wrapper script for Libglvnd_jll for aarch64-linux-gnu
export libEGL, libGL, libGLESv1_CM, libGLESv2, libGLX, libGLdispatch, libOpenGL

using Xorg_libX11_jll
using Xorg_libXext_jll
JLLWrappers.@generate_wrapper_header("Libglvnd")
JLLWrappers.@declare_library_product(libEGL, "libEGL.so.1")
JLLWrappers.@declare_library_product(libGL, "libGL.so.1")
JLLWrappers.@declare_library_product(libGLESv1_CM, "libGLESv1_CM.so.1")
JLLWrappers.@declare_library_product(libGLESv2, "libGLESv2.so.2")
JLLWrappers.@declare_library_product(libGLX, "libGLX.so.0")
JLLWrappers.@declare_library_product(libGLdispatch, "libGLdispatch.so.0")
JLLWrappers.@declare_library_product(libOpenGL, "libOpenGL.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_libX11_jll, Xorg_libXext_jll)
    JLLWrappers.@init_library_product(
        libEGL,
        "lib/libEGL.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGL,
        "lib/libGL.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGLESv1_CM,
        "lib/libGLESv1_CM.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGLESv2,
        "lib/libGLESv2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGLX,
        "lib/libGLX.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGLdispatch,
        "lib/libGLdispatch.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libOpenGL,
        "lib/libOpenGL.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module OptimTools

	using KrylovMethods

	include("checkDerivative.jl")
	include("lineSearch.jl")
	include("newton.jl")
	include("dampedNewton.jl")
	include("sd.jl")
	include("bfgs.jl")
	include("nlcg.jl")
	include("modnewton.jl")
	include("newtoncg.jl")
	include("sr1.jl")
	include("tr.jl")
end


text:
#!/usr/bin/env julia
using Base.Filesystem

struct Meta
	artist::Vector{String}
	albumname::Vector{String}
	albumyear::Vector{String}
	trackname::Vector{String}
	tracknum::Vector{String}
end
Meta() = Meta(Vector{String}(), Vector{String}(), Vector{String}(), Vector{String}(), Vector{String}())

"""
    parse_pairfile(path::AbstractString)

Parse a pair file into a Dict.

A pair file consists of lines containing `key<tab>value` pairs.
"""
function parse_pairfile(path::AbstractString)
	dict = Dict{Symbol, String}()
	for line in readlines(path)
		rmatch = match(r"(\w+)\t(.+)", line)
		dict[Symbol(rmatch.captures[1])] = rmatch.captures[2]
	end
	return dict
end

if isempty(ARGS)
	println("Called with no arguments, assuming current working directory.")
	rootdir = pwd()
else
	rootdir = realpath(ARGS[1])
end
meta = Meta()
# Set the default artist, specified by meta file or deduced from rootdir name otherwise.
pushfirst!(meta.artist, get(parse_pairfile(joinpath(rootdir, "meta")), :artist, basename(rootdir)))
globalmeta = parse_pairfile(joinpath(rootdir, "meta"))

for (path, dirs, files) in walkdir(rootdir)
	println("Entering $path")
	dirmatch = match(r"(.+)\s*\((\d{4})\)", basename(path))
	if dirmatch != nothing
		albumname = dirmatch.captures[1]
		albumyear = dirmatch.captures[2]
	else
		albumname = basename(path)
		albumyear = nothing
	end
	@assert albumname != nothing
	@info("albumname: $albumname")
	for file in files
		fullpath = joinpath(path, file)
		println("Processing $fullpath")
		filematch = match(r"(\d+)\s*-\s*([^\.]+)(?:\..+)?", file)
		if filematch != nothing
			tracknum = filematch.captures[1]
			trackname = filematch.captures[2]
		else
			tracknum = nothing
			trackname = match(r"(.+)(?:\..+)?", file).captures[1]
		end
		converter_args = ["--artist", artist, "--song", trackname, "--album", albumname]
		if albumyear != nothing
			push!(converter_args, "--year", albumyear)
		end
		if tracknum != nothing
			push!(converter_args, "--track", tracknum)
		end
		println.(converter_args)
		run(`id3v2 $converter_args "$fullpath"`)
	end
end


text:
@testset "Test Normalisation" begin
    @info "Test: Ellipse Parameter Cooordinate System Conversion"

    @testset "From Unnormalised to Normalised and back" begin
        𝛏₀ = [10, 5, 25, 25, π/4]
        geo_to_alg = GeometricToAlgebraic()
        s = 0.01
        m₁ = 25
        m₂ = 25
        𝐓 = SMatrix{3,3,Float64}(s,0,0,0,s,0,-s*m₁,-s*m₂,1)
        normalise = NormaliseDataContext(tuple(𝐓))
        𝛉₀ = geo_to_alg(𝛏₀)
        𝛉₀ = 𝛉₀ / norm(𝛉₀)
        𝛉′ = normalise(ToNormalisedSpace(), 𝛉₀)
        𝛉₁ = normalise(FromNormalisedSpace(), 𝛉′)
        @test all(𝛉₀ .≈ 𝛉₁)
    end
end


text:
# Generating data from two two-dimensional gaussian processes
using Distributions
using FixedPointAcceleration
using Random
using DataFrames
using Plots
using LinearAlgebra
true_tau = 0.6
nobs_1 = 400
nobs_2 = convert(Int, round(nobs_1 * ((1-true_tau)/true_tau)))
Random.seed!(1234)
mu_1 = [0.0,8.0]
cov_1 = [2.0,0.5,2.0]
covar_1 = Symmetric([cov_1[1] cov_1[2]; cov_1[2] cov_1[3]])
md_1 = MultivariateNormal(mu_1,covar_1)
mu_2 = [-4.0,10.0]
cov_2 = [2.0,-0.75,12.0]
covar_2 = Symmetric([cov_2[1] cov_2[2]; cov_2[2] cov_2[3]])
md_2 = MultivariateNormal(mu_2,covar_2)

rands_from_1 = transpose(rand(md_1, nobs_1))
rands_from_2 = transpose(rand(md_2, nobs_2))
data1 = DataFrame([rands_from_1[:,1], rands_from_1[:,2]], [:x1, :x2])
data2 = DataFrame([rands_from_2[:,1], rands_from_2[:,2]], [:x1, :x2])
dd  = vcat(data1,data2)
# Plotting it:
plot(data1.x1, data1.x2,seriestype=:scatter)
plot!(data2.x1, data2.x2,seriestype=:scatter)

function z_estimate_given_theta(x::Array{Float64,1}, md_1::MultivariateNormal, md_2::MultivariateNormal, tau::Float64)
    pdf_1 = pdf(md_1, x)
    pdf_2 = pdf(md_2, x)
    return tau*pdf_1 / (tau*pdf_1 + (1-tau)*pdf_2)
end

function update_tau(Z::Array{Float64,1})
    return mean(Z)
end

function update_mu(dd::DataFrame, Z::Array{Float64,1})
    X = convert(Array{Float64,2}, dd[[:x1, :x2]])
    sum_Z = sum(Z)
    updated_mu = (transpose(Z) * X) ./sum_Z
    return vec(updated_mu)
end

function update_cov(dd::DataFrame, updated_mu::Array{Float64,1}, Z::Array{Float64,1})
    X_minus_mu = convert(Array{Float64,2}, dd[[:x1, :x2]]) .- transpose(updated_mu)
    sum_Z = sum(Z)
    updated_cov = (transpose(Z .* X_minus_mu) * X_minus_mu) ./sum_Z
    return [updated_cov[1,1], updated_cov[1,2], updated_cov[2,2]]
end

function update_theta(theta::Array{Float64,1}, dd::DataFrame)
    # We will use the convention that theta's 11 entries are (mu_1, cov_1, mu_2, cov_2, tau). First unpacking theta:
    mu_1    = theta[[1,2]]
    cov_1   = theta[[3,4,5]]
    covar_1 = Symmetric([cov_1[1] cov_1[2]; cov_1[2] cov_1[3]])
    md_1 = MultivariateNormal(mu_1,covar_1)
    mu_2    = theta[[6,7]]
    cov_2   = theta[[8,9,10]]
    covar_2 = Symmetric([cov_2[1] cov_2[2]; cov_2[2] cov_2[3]])
    md_2 = MultivariateNormal(mu_2,covar_2)
    tau     = theta[11]
    # Getting Z
    Z = Array{Float64,1}(undef,size(dd)[1])
    for i in 1:size(dd)[1]
        Z[i] = z_estimate_given_theta([dd[i,:x1], dd[i,:x2]], md_1, md_2, tau)
    end

    # Updating Tau
    updated_tau = update_tau(Z)
    # Updating mu1
    updated_mu_1 = update_mu(dd,Z)
    updated_mu_2 = update_mu(dd, 1 .- Z)
    # Updating Cov
    updated_cov_1 = update_cov(dd, updated_mu_1, Z)
    updated_cov_2 = update_cov(dd, updated_mu_2, 1 .- Z)
    # Returning theta
    updated_theta = vcat(updated_mu_1, updated_cov_1, updated_mu_2, updated_cov_2, updated_tau)
    side_product = (z = Z,)
    return (updated_theta, side_product)
end

InitialGuess = [0.5, 7.5, 2.0, 0.0, 2.0, -5.0, 7.5, 2.0, 0.0, 10.0, 0.5]
fp_anderson = fixed_point(x -> update_theta(x,dd), InitialGuess; Algorithm = Anderson, PrintReports = true)
isa(fp_anderson.Other_Output_, NamedTuple)


text:
# Autogenerated wrapper script for rubberband_jll for x86_64-w64-mingw32-cxx11
export librubberband, rubberband

using FFTW_jll
using libsamplerate_jll
using libsndfile_jll
JLLWrappers.@generate_wrapper_header("rubberband")
JLLWrappers.@declare_library_product(librubberband, "librubberband.dll")
JLLWrappers.@declare_executable_product(rubberband)
function __init__()
    JLLWrappers.@generate_init_header(FFTW_jll, libsamplerate_jll, libsndfile_jll)
    JLLWrappers.@init_library_product(
        librubberband,
        "bin\\librubberband.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        rubberband,
        "bin\\rubberband.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Tests formulations in prob/scopf
@testset "test scopf" begin

scopf_cont_cuts_dc_objective = [14642.16, 32952.18]
@testset "scopf contigency cuts dc, from file - $(i)" for (i,network) in enumerate(networks)
    result = run_scopf_contigency_cuts(ini_file, DCPPowerModel, lp_solver, scenario_id=scenarios[i])

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], scopf_cont_cuts_dc_objective[i]; atol = 1e0)
end

@testset "scopf contigency cuts dc - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    result = run_scopf_contigency_cuts(network, DCPPowerModel, lp_solver)

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], scopf_cont_cuts_dc_objective[i]; atol = 1e0)
end

scopf_cont_cuts_ac_objective = [14676.95, 33281.21]
@testset "scopf contigency cuts ac - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    result = run_scopf_contigency_cuts(network, ACPPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], scopf_cont_cuts_ac_objective[i]; atol = 1e0)
end

#scopf_cont_cuts_ac_objective = [14676.95, 33281.21]
@testset "scopf contigency cuts ac - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    result = run_scopf_contigency_cuts(network, ACRPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], scopf_cont_cuts_ac_objective[i]; atol = 1e0)
end

scopf_ptdf_cuts_dc_objective = [14642.16, 37233.43]
@testset "scopf ptdf cuts dc, from file - $(i)" for (i,network) in enumerate(networks)
    result = run_scopf_ptdf_cuts(ini_file, DCPPowerModel, lp_solver, scenario_id=scenarios[i])

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], scopf_ptdf_cuts_dc_objective[i]; atol = 1e0)
end

@testset "scopf ptdf cuts dc - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_scopf_ptdf_cuts!(network, DCPPowerModel, lp_solver)

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], scopf_ptdf_cuts_dc_objective[i]; atol = 1e0)
end

scopf_ptdf_cuts_ac_objective = [14676.95, 37808.75]
@testset "scopf ptdf cuts ac - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    result = run_scopf_ptdf_cuts!(network, ACPPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], scopf_ptdf_cuts_ac_objective[i]; atol = 1e0)
end

#scopf_ptdf_cuts_ac_objective = [14676.95, 37808.75]
@testset "scopf ptdf cuts ac - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    result = run_scopf_ptdf_cuts!(network, ACRPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], scopf_ptdf_cuts_ac_objective[i]; atol = 1e0)
end

scopf_dc_cuts_soft_woc_objective = [14642.16, 26982.17]
@testset "scopf cuts dc soft, without cuts - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_scopf_cuts_soft(network, DCPPowerModel, lp_solver)

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], scopf_dc_cuts_soft_woc_objective[i]; atol = 1e0)
end

@testset "scopf cuts dc soft, infeasible" begin
    network = deepcopy(network_infeasible)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_scopf_cuts_soft(network, DCPPowerModel, lp_solver)
    @test isapprox(result["termination_status"], INFEASIBLE)
end

scopf_dc_cuts_soft_wc_objective = [14642.16, 37403.76]
@testset "scopf cuts dc soft, with cuts - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    result = run_opf_cheap(network, DCPPowerModel, lp_solver)
    @test isapprox(result["termination_status"], OPTIMAL)

    update_active_power_data!(network, result["solution"])

    cuts = check_contingencies_branch_power(network, total_cut_limit=2, gen_flow_cuts=[], branch_flow_cuts=[])

    #println(length(cuts.gen_cuts) + length(cuts.branch_cuts))
    #cuts_found = sum(length(c.gen_cuts)+length(c.branch_cuts) for c in cuts)
    network["gen_flow_cuts"] = cuts.gen_cuts
    network["branch_flow_cuts"] = cuts.branch_cuts

    result = run_scopf_cuts_soft(network, DCPPowerModel, lp_solver)

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], scopf_dc_cuts_soft_wc_objective[i]; atol = 1e0)
end

scopf_ac_cuts_soft_wc_objective = [14676.95, 37904.03]
@testset "scopf cuts ac soft, with cuts - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    result = run_opf_cheap(network, ACPPowerModel, nlp_solver)
    @test isapprox(result["termination_status"], LOCALLY_SOLVED)

    update_active_power_data!(network, result["solution"])

    cuts = check_contingencies_branch_power(network, total_cut_limit=2, gen_flow_cuts=[], branch_flow_cuts=[])

    #println(length(cuts.gen_cuts) + length(cuts.branch_cuts))
    #cuts_found = sum(length(c.gen_cuts)+length(c.branch_cuts) for c in cuts)
    network["gen_flow_cuts"] = cuts.gen_cuts
    network["branch_flow_cuts"] = cuts.branch_cuts

    result = run_scopf_cuts_soft(network, ACPPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], scopf_ac_cuts_soft_wc_objective[i]; atol = 1e0)
end

scopf_dc_cuts_soft_woc_objective = [14642.16, 26982.17]
@testset "scopf cuts dc soft bpv, without cuts - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_scopf_cuts_soft_bpv(network, DCPPowerModel, lp_solver)

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], scopf_dc_cuts_soft_woc_objective[i]; atol = 1e0)
end

@testset "scopf cuts dc soft bpv, infeasible" begin
    network = deepcopy(network_infeasible)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_scopf_cuts_soft_bpv(network, DCPPowerModel, lp_solver)
    @test isapprox(result["termination_status"], INFEASIBLE)
end

scopf_dc_cuts_soft_wc_objective = [14642.16, 27258.16]
@testset "scopf cuts dc soft bpv, with cuts - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    result = run_opf_cheap(network, DCPPowerModel, lp_solver)
    @test isapprox(result["termination_status"], OPTIMAL)

    update_active_power_data!(network, result["solution"])

    cuts = check_contingencies_branch_power_bpv(network, total_cut_limit=2, gen_flow_cuts=[], branch_flow_cuts=[])

    #println(length(cuts.gen_cuts), " ", length(cuts.branch_cuts))
    #cuts_found = sum(length(c.gen_cuts)+length(c.branch_cuts) for c in cuts)
    network["gen_flow_cuts"] = cuts.gen_cuts
    network["branch_flow_cuts"] = cuts.branch_cuts

    result = run_scopf_cuts_soft_bpv(network, DCPPowerModel, lp_solver)

    @test isapprox(result["termination_status"], OPTIMAL)
    @test isapprox(result["objective"], scopf_dc_cuts_soft_wc_objective[i]; atol = 1e0)
end

scopf_ac_cuts_soft_wc_objective = [14676.95, 28318.96]
@testset "scopf cuts ac soft bpv, with cuts - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)

    result = run_opf_cheap(network, ACPPowerModel, nlp_solver)
    @test isapprox(result["termination_status"], LOCALLY_SOLVED)

    update_active_power_data!(network, result["solution"])

    cuts = check_contingencies_branch_power_bpv(network, total_cut_limit=2, gen_flow_cuts=[], branch_flow_cuts=[])

    #println(length(cuts.gen_cuts), " ", length(cuts.branch_cuts))
    #cuts_found = sum(length(c.gen_cuts)+length(c.branch_cuts) for c in cuts)
    network["gen_flow_cuts"] = cuts.gen_cuts
    network["branch_flow_cuts"] = cuts.branch_cuts

    result = run_scopf_cuts_soft_bpv(network, ACPPowerModel, nlp_solver)

    @test isapprox(result["termination_status"], LOCALLY_SOLVED)
    @test isapprox(result["objective"], scopf_ac_cuts_soft_wc_objective[i]; atol = 1e0)
end

end # close test group


text:
"""
    ImplicitDistribution(sample_function, args...)

Define a distribution that can only be sampled from using `rand`, but has no explicit `pdf`.

Each time `rand(rng, d::ImplicitDistribution)` is called,
```julia
sample_function(args..., rng)
```
will be called to generate a new sample.

`ImplicitDistribution` is designed to be used with anonymous functions or the `do` syntax as follows:

# Examples

```julia
ImplicitDistribution(rng->rand(rng)^2)
```

```julia
struct MyMDP <: MDP{Float64, Int} end

function POMDPs.transition(m::MyMDP, s, a)
    ImplicitDistribution(s, a) do s, a, rng
        return s + a + 0.001*randn(rng)
    end
end

td = transition(MyMDP(), 1.0, 1)
rand(td) # will return a number near 2
```
"""
struct ImplicitDistribution{F, A}
    f::F
    args::A

    # internal constructor needed for single argument case
    ImplicitDistribution(f, args...) = new{typeof(f), typeof(args)}(f, args)
end

function Base.rand(rng::AbstractRNG, s::Random.SamplerTrivial{<:ImplicitDistribution})
    d = s[]
    d.f(d.args..., rng)
end


text:
const Vec{N, T} = NTuple{N, T}
const Point{N, T} = NTuple{N, T}

struct ISF{
        IntType, FloatType,
        # can't define the Types as e.g. ArrayType{FloatType, 3}, so we need them as type parameter and calculate the types in the constructor
        CArray, FArray, VIArray, VFArray, # Complex array, Float array, Integer Vec array, Float Vec array
        FFTP, IFFTP # types of the fft plan
    }

    grid_res::Vec{3, IntType}
    physical_size::Vec{3, IntType}
    d::Vec{3, FloatType}
    velocity::VFArray

    hbar::FloatType             # reduced Planck constant
    dt::FloatType               # time step
    mask::CArray # Fourier coefficient for solving Schroedinger eq
    # precalculated values
    # we have a couple of indices precalculated:
    i::VIArray # type is to complex to write down
    i_shifted::VIArray # indices shifted by one
    i_shifted2::VIArray # indices shifted by -2

    # other precalculated values
    positions::VFArray # position on the grid
    fac::CArray
    # temporary allocations for intermediates
    f::FArray
    fc::CArray
    psi1::CArray
    psi2::CArray
    fftplan::FFTP
    ifftplan::IFFTP
end

function Base.zero(::Type{NTuple{N, T}}) where {T, N}
    ntuple(x-> zero(T), Val{N})
end

function (::Type{ISF{ArrayType, IntType, FloatType}})(physical_size, dims, hbar, dt) where {ArrayType, IntType, FloatType}
    VF = Vec{3, FloatType}
    VI = Vec{3, IntType}
    grid_res = VI(dims)
    physical_size = VI(physical_size)
    d = VF(physical_size ./ grid_res)
    fac = zeros(Complex{FloatType}, dims)
    mask = zeros(Complex{FloatType}, dims)
    f = -4 * pi^2 * hbar
    for x = 1:dims[1], y = 1:dims[2], z = 1:dims[3]
        xyz = (x, y, z)
        # fac
        s = sin.(pi .* (xyz .- 1) ./ grid_res) ./ d
        denom = sum(s .^ 2)
        fac[x,y,z] = -0.25f0 ./ denom

        # schroedingers mask
        k = (xyz .- 1 .- grid_res ./ 2) ./ physical_size
        lambda = f .* sum(k .^ 2)
        mask[x,y,z] = exp(1.0im * lambda * dt / 2)
    end
    mask = fftshift(mask)
    fac[1,1,1] = 0

    f = zeros(FloatType, dims)
    velocity = zeros(VF, dims)

    i = collect(((x, y, z) for x=1:dims[1], y=1:dims[2], z=1:dims[3]))
    i_shifted = map(i) do i
        mod.(Vec(i), Vec(dims)) .+ 1
    end
    i_shifted2 = map(i) do i
         mod.(i .- 2, dims) .+ 1
    end
    positions = map(i) do i
        (i .- 1) .* d
    end
    fac_gpu = ArrayType(fac)
    p1, p2 = plan_fft!(fac_gpu), plan_ifft!(fac_gpu)
    ISF{
            IntType, FloatType,
            ArrayType{Complex{FloatType}, 3}, ArrayType{FloatType, 3},
            ArrayType{VI, 3}, ArrayType{VF, 3},
            typeof(p1), typeof(p2)
        }(
        grid_res, physical_size, d,
        velocity,
        hbar, dt,
        mask,
        i, i_shifted, i_shifted2,
        positions,
        fac_gpu, f,
        copy(fac_gpu), # temporaries
        copy(fac_gpu),
        copy(fac_gpu),
        p1, p2
    )
end

@inline function velocity_one_form(i, i2, psi, hbar)
    @inbounds begin
        psi12  = psi[i[1],  i[2],  i[3]]
        psix12 = psi[i2[1], i[2],  i[3]]
        psiy12 = psi[i[1],  i2[2] ,i[3]]
        psiz12 = psi[i[1],  i[2],  i2[3]]
        psi1n = (psix12[1], psiy12[1], psiz12[1])
        psi2n = (psix12[2], psiy12[2], psiz12[2])
        return angle.(
            conj(psi12[1]) .* psi1n .+
            conj(psi12[2]) .* psi2n
        ) .* hbar
    end
end

function velocity_one_form!(isf, psi, hbar = 1.0f0)
    isf.velocity .= velocity_one_form.(
        isf.i,
        isf.i_shifted,
        (psi,),
        hbar
    )
end

@fastmath function gauge_transform(psi, q)
    x = Complex64(0f0, 1f0) * (-q)
    eiq = exp(x)
    (psi[1] * eiq, psi[2] * eiq)
end

function div(i, i2, velocity, res, ds)
    x = i[1]; y = i[2]; z = i[3]
    ix = i2[1]; iy = i2[2]; iz = i2[3]
    v1 = velocity[x, y,  z]
    v2 = (
        velocity[ix, y,  z][1],
        velocity[x,  iy, z][2],
        velocity[x,  y, iz][3]
    )
    t1 = v1 .- v2
    t2 = t1 .* ds
    sum(t2)
end
function div!(isf)
    isf.f .= div.(
        isf.i,
        isf.i_shifted2,
        (isf.velocity,),
        (isf.grid_res,),
        (inv.(isf.d .^ 2f0),)
    )
end

function pressure_project!(isf, psi)
    velocity_one_form!(isf, psi)
    div!(isf)
    q = poisson_solve(isf)
    psi .= gauge_transform.(psi, q)
end

function normalize_psi(psi)
    norm = @fastmath hypot(abs(psi[1]), abs(psi[2]))
    (psi[1] / norm, psi[2] / norm)
end
@inline twotuple(a, b) = (a, b)

function poisson_solve(isf)
    isf.fc .= complex.(isf.f)
    isf.fftplan * isf.fc
    isf.fc .= (*).(isf.fc, isf.fac)
    isf.ifftplan * isf.fc
    isf.fc
end

function schroedinger_flow!(isf, psi)
    # extract single psi values into tmp arrays stored in obj
    psi1 = isf.psi1; psi2 = isf.psi2
    psi1 .= getindex.(psi, 1)
    psi2 .= getindex.(psi, 2)
    isf.fftplan * psi1; isf.fftplan * psi2;
    psi1 .= (*).(psi1, isf.mask)
    psi2 .= (*).(psi2, isf.mask)
    isf.ifftplan * psi1; isf.ifftplan * psi2;
    psi .= twotuple.(psi1, psi2) # convert back
    psi
end

function in_grid(p, physical_size)
    p[1] > 0.0f0 && p[1] < physical_size[1] || return false
    p[2] > 0.0f0 && p[2] < physical_size[2] || return false
    p[3] > 0.0f0 && p[3] < physical_size[3] || return false
    true
end

function staggered_advect(p, velocity, args)
    dt = args[1]; d = args[2]; gridsize = args[3]; res = args[4]
    if in_grid(p, gridsize)

        k1 = staggered_velocity(velocity, p, d, gridsize, res)

        k2 = p .+ k1 .* dt .* 0.5f0
        k2 = staggered_velocity(velocity, k2, d, gridsize, res)

        k3 = p .+ k2 .* dt .* 0.5f0
        k3 = staggered_velocity(velocity, k3, d, gridsize, res)

        k4 = p .+ k3 .* dt
        k4 = staggered_velocity(velocity, k4, d, gridsize, res)

        p .+ dt./6f0 .* (k1 .+ 2f0 .* k2 .+ 2f0 .* k3 .+ k4)
    else
        (-1f0, -1f0, -1f0)
    end
end

@inline function staggered_velocity(velocity, point, d, ph_size, res)
    p   = mod.(point, ph_size)
    i   = UInt32.(floor.(p ./ d)) .+ UInt32(1)
    ip  = mod.(i, res) .+ UInt32(1)

    v0  = velocity[i[1], i[2], i[3]]

    pxp = velocity[ip[1], i[2], i[3]]
    pyp = velocity[i[1], ip[2], i[3]]
    pzp = velocity[i[1], i[2], ip[3]]

    vn = (
        velocity[i[1], ip[2], ip[3]][1],
        velocity[ip[1], i[2], ip[3]][2],
        velocity[ip[1], ip[2], i[3]][3]
    )
    pp  = (pyp[1], pxp[2], pxp[3])
    pp2 = (pzp[1], pzp[2], pyp[3])

    w   = p .- (i .- 1f0) .* d
    w1  = (w[3], w[3], w[2])
    w2  = (w[2], w[1], w[1])

    return (
        (1f0 .- w1) .* ((1f0 .- w2) .* v0 .+ w2 .* pp) .+
        w1 .* ((1f0 .- w2) .* pp2 .+ w2 .* vn)
    )
end

function staggered_advect!(particle, isf)
    velocity, dt, d, res = isf.velocity, isf.dt, isf.d, isf.grid_res
    gridsize = isf.physical_size
    particle .= staggered_advect.(
        particle,
        (velocity,),
        ((
            isf.dt,
            isf.d,
            gridsize,
            isf.grid_res
        ),)
    )
end

export ISF, normalize_psi, pressure_project!
export velocity_one_form!, schroedinger_flow!
export Particles, staggered_advect!


text:
# TODO: look at  the following reference for E-free linesearch:
#  Fletcher R. (2005) On the Barzilai-Borwein Method. Applied Optimization,
#   vol 96. Springer.

# TODO: replace with with LineSearches or alternatively move into LineSearches

"""
`function  parab3p(lambdac, lambdam, ff0, ffc, ffm) -> lambdap`

Apply three-point safeguarded parabolic model for a line search.

### Input Parameters

* `lambdac` : current steplength
* `lambdam` : previous steplength
* `ff0` : value of |F(x_c)|^2
* `ffc` : value of |F(x_c + λc d)|^2
* `ffm` : value of |F(x_c + λm d)|^2

### Output

* lambdap : new value of lambda given parabolic model

### Keyword arguments (internal parameters)

* `sigma0 = .1`, `sigma1 = .5` : safeguarding bounds for the linesearch
"""
function parab3p(lambdac, lambdam, ff0, ffc, ffm; sigma0 = 0.1, sigma1 = 0.5)
   # compute coefficients of interpolation polynomial
   # p(lambda) = ff0 + (c1 lambda + c2 lambda^2)/d1
   # d1 = (lambdac - lambdam)*lambdac*lambdam < 0
   #      so if c2 > 0 we have negative curvature and default to
   #      lambdap = sigam1 * lambda
   c2 = lambdam * (ffc - ff0) - lambdac * (ffm - ff0)
   if c2 >= 0
      return sigma1 * lambdac
   end
   c1 = lambdac^2 * (ffm-ff0) - lambdam^2 * (ffc - ff0)
   lambdap = -c1 * 0.5 / c2
   if lambdap < sigma0 * lambdac
      return sigma0 * lambdac
   elseif lambdap > sigma1 * lambdac
      return sigma1 * lambdac
   end
   return lambdap
end

function parab2p(f0, df0, f1, alpha)
   return - (df0 * alpha^2) / ( 2.0 * (f1 - f0 - df0*alpha) )
end

# """
# at the moment, this is just quadratic back-tracking
#
# TODO:
#  * probably switch to LineSearches.jl
# """
# function lsarmijo(x, p, α_old, E, dE, dE0, P, maxstep, K = nothing;
#                   Ca = 0.2)
#
#    maxα = maxstep / norm(p, Inf)
#    # take one potentially forward step
#    slope = dot(dE0, p)
#    E0 = E(x)
#    #  α1 = 0.5 + 0.33*α_old : works really  well for P = I
#    #  α1 = 0.66*α_old : works really well with a good P
#    #  α1 = 1.0 : decent overall
#    if P == I
#       α1 = 0.5 + 0.33 * α_old
#    else
#       α1 = 0.66 * α_old
#    end
#    α1 = min(α1, maxα)
#    x1 = x + α1 * p
#    E1 = E(x1)
#    numdE = 1
#
#    if E1 <= E0 + Ca * α1 * slope
#       # Armijo is satisfied, we will now try to improve a bit but won't
#       # try too hard and just revert to x1 if needed
#       if E0 + α1 * slope < E1   # above the slope > parab2p is meaningful
#          αt = parab2p(E0, slope, E1, α1)
#       else # below the slope > parab2p would give nonsense
#          αt = 1.5 * max(α1, α_old)
#       end
#       # make sure the step is not too small
#       αt = max(α1/4, αt)
#       # and not too large
#       αt = min(α1*2, maxα, αt)
#       # update configuration
#       xt = x + αt * p
#       Et = E(xt)
#       numdE += 1
#       # check whether to return xt or x1
#       if Et > E1
#          Et, αt, xt = E1, α1, x1
#       end
#    else
#       # CASE 2: the starting guess does not satisfy Armijo
#       # > start a back-tracking loop
#       xt, αt, Et = x1, α1, E1
#       while Et > E0 + Ca * αt * slope
#          αt_old, αt = αt, parab2p(E0, slope, Et, αt)
#          αt = max(αt_old/4, αt)
#          xt = x + αt * p
#          Et = E(xt)
#          numdE += 1
#       end
#    end
#
#    # we also need to return the gradient at the final l.s. point
#    dEt = dE(xt)
#    nft = nkdualnorm(P, dEt)
#    return αt, xt, dE(xt), nft, numdE
# end
#
#
# """
# A very crude step-length selection for the pre-asymptotic regime
# """
# function lsforce(x, p, α_old, E, dE, f0, P, maxstep, K;
#             Ca = 0.2, Cw = 0.9, minα = 1e-6, αinit = :one)
#    # in this case, we do something very crude:
#    #   take same step as before, then do one line-search step
#    #   and pick the better of the two.
#    φ = f_ -> dot(K \ f_, P, p)  # dot(f_, p)   # TODO: not sure which one will be better
#    g0 = φ(f0)
#    tol = Cw * abs(g0)
#    maxα = maxstep / norm(p, Inf)
#    numdE = 0
#    if αinit == :one
#       α1 = 1.0
#    elseif αinit == :heuristic
#       if P == I
#          α1 = 0.5 + 0.33 * α_old
#       else
#          α1 = 0.66 * α_old
#       end
#    else
#       error("unknown `αinit`")
#    end
#    α1 = min(α1, maxα)
#    x1 = x + α1 * p
#    f1 = dE(x1)
#    g1 = φ(f1)
#    numdE += 1
#
#    acceptstep = gg -> 2 * g0 < gg < 0.2 * abs(g0)
#
#    if acceptstep(g1)
#       # put a straight line through (0, g0), (α1, g1) and intersect with 0
#       if abs(g1-g0)/α1 < 1e-3
#          αt = 0.5 * α1
#       else
#          αt = - g0 * α1 / (g1-g0)   # g0 + α * (g1-g0)/α1 = 0
#       end
#       # make sure the step is not too small
#       αt = max(α1/4, αt)
#       # and not too large
#       αt = min(α1*2, maxα, αt)
#
#       # update configuration
#       xt = x + αt * p
#       ft = dE(xt)
#       gt = φ(ft)
#       numdE += 1
#       # check whether to return xt or x1
#       if !acceptstep(gt) || abs(gt) > abs(g1)
#          gt, αt, xt = g1, α1, x1
#       end
#    else
#       αt, gt, ft, xt = α1, g1, f1, x1
#       # reduce step-length until the new gradient is not too bad
#       while (2 * g0 > gt) || (gt > 0.2 * abs(g0))
#          αt *= 0.5
#          xt = x + αt * p
#          ft = dE(xt)
#          numdE += 1
#          gt = φ(ft)
#       end
#    end
#
#    return αt, xt, ft, dualnorm(P, ft), numdE
# end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    AbstractDateToken

A token used in parsing or formatting a date time string. Each subtype must
define the tryparsenext and format methods.

"""
abstract type AbstractDateToken end

"""
    tryparsenext(tok::AbstractDateToken, str::String, i::Int, len::Int, locale::DateLocale)

`tryparsenext` parses for the `tok` token in `str` starting at index `i`.
`len` is the length of the string.  parsing can be optionally based on the
`locale`. If a `tryparsenext` method does not need a locale, it can leave
the argument out in the method definition.

If parsing succeeds, returns a tuple of 2 elements `(res, idx)`, where:

* `res` is the result of the parsing.
* `idx::Int`, is the index _after_ the index at which parsing ended.
"""
function tryparsenext end

"""
    format(io::IO, tok::AbstractDateToken, dt::TimeType, locale)

Format the `tok` token from `dt` and write it to `io`. The formatting can
be based on `locale`.

All subtypes of `AbstractDateToken` must define this method in order
to be able to print a Date / DateTime object according to a `DateFormat`
containing that token.
"""
function format end

# fallback to tryparsenext/format methods that don't care about locale
@inline function tryparsenext(d::AbstractDateToken, str, i, len, locale)
    return tryparsenext(d, str, i, len)
end

function Base.string(t::Time)
    h, mi, s = hour(t), minute(t), second(t)
    hh = lpad(h, 2, "0")
    mii = lpad(mi, 2, "0")
    ss = lpad(s, 2, "0")
    nss = tons(Millisecond(t)) + tons(Microsecond(t)) + tons(Nanosecond(t))
    ns = nss == 0 ? "" : rstrip(@sprintf("%.9f", nss / 1e+9)[2:end], '0')
    return "$hh:$mii:$ss$ns"
end

Base.show(io::IO, x::Time) = print(io, string(x))

@inline function format(io, d::AbstractDateToken, dt, locale)
    format(io, d, dt)
end

# Information for parsing and formatting date time values.
struct DateFormat{S, T<:Tuple}
    tokens::T
    locale::DateLocale
end

### Token types ###

struct DatePart{letter} <: AbstractDateToken
    width::Int
    fixed::Bool
end

@inline min_width(d::DatePart) = d.fixed ? d.width : 1
@inline max_width(d::DatePart) = d.fixed ? d.width : 0

function _show_content(io::IO, d::DatePart{c}) where c
    for i = 1:d.width
        print(io, c)
    end
end

function Base.show(io::IO, d::DatePart{c}) where c
    print(io, "DatePart(")
    _show_content(io, d)
    print(io, ")")
end

### Parse tokens

for c in "yYmdHMS"
    @eval begin
        @inline function tryparsenext(d::DatePart{$c}, str, i, len)
            return tryparsenext_base10(str, i, len, min_width(d), max_width(d))
        end
    end
end

for (tok, fn) in zip("uUeE", [monthabbr_to_value, monthname_to_value, dayabbr_to_value, dayname_to_value])
    @eval @inline function tryparsenext(d::DatePart{$tok}, str, i, len, locale)
        next = tryparsenext_word(str, i, len, locale, max_width(d))
        next === nothing && return nothing
        word, i = next
        val = $fn(word, locale)
        val == 0 && return nothing
        return val, i
    end
end

# 3-digit (base 10) number following a decimal point. For InexactError below.
struct Decimal3 end

@inline function tryparsenext(d::DatePart{'s'}, str, i, len)
    val = tryparsenext_base10(str, i, len, min_width(d), max_width(d))
    val === nothing && return nothing
    ms0, ii = val
    len = ii - i
    if len > 3
        ms, r = divrem(ms0, Int64(10) ^ (len - 3))
        r == 0 || throw(InexactError(:convert, Decimal3, ms0))
    else
        ms = ms0 * Int64(10) ^ (3 - len)
    end
    return ms, ii
end

### Format tokens

for (c, fn) in zip("YmdHMS", [year, month, day, hour, minute, second])
    @eval function format(io, d::DatePart{$c}, dt)
        print(io, string($fn(dt), base = 10, pad = d.width))
    end
end

for (tok, fn) in zip("uU", [monthabbr, monthname])
    @eval function format(io, d::DatePart{$tok}, dt, locale)
        print(io, $fn(month(dt), locale))
    end
end

for (tok, fn) in zip("eE", [dayabbr, dayname])
    @eval function format(io, ::DatePart{$tok}, dt, locale)
        print(io, $fn(dayofweek(dt), locale))
    end
end

@inline function format(io, d::DatePart{'y'}, dt)
    y = year(dt)
    n = d.width

    # the last n digits of y
    # will be 0 padded if y has less than n digits
    str = string(y, base = 10, pad = n)
    l = lastindex(str)
    if l == n
        # fast path
        print(io, str)
    else
        print(io, SubString(str, l - (n - 1), l))
    end
end

function format(io, d::DatePart{'s'}, dt)
    ms = millisecond(dt)
    if ms % 100 == 0
        str = string(div(ms, 100))
    elseif ms % 10 == 0
        str = string(div(ms, 10), pad = 2)
    else
        str = string(ms, pad = 3)
    end

    print(io, rpad(str, d.width, '0'))
end

### Delimiters

struct Delim{T, length} <: AbstractDateToken
    d::T
end

Delim(d::T) where {T<:AbstractChar} = Delim{T, 1}(d)
Delim(d::String) = Delim{String, length(d)}(d)

@inline function tryparsenext(d::Delim{<:AbstractChar, N}, str, i::Int, len) where N
    for j = 1:N
        i > len && return nothing
        next = iterate(str, i)
        @assert next !== nothing
        c, i = next
        c != d.d && return nothing
    end
    return true, i
end

@inline function tryparsenext(d::Delim{String, N}, str, i::Int, len) where N
    i1 = i
    i2 = firstindex(d.d)
    for j = 1:N
        if i1 > len
            return nothing
        end
        next1 = iterate(str, i1)
        @assert next1 !== nothing
        c1, i1 = next1
        next2 = iterate(d.d, i2)
        @assert next2 !== nothing
        c2, i2 = next2
        if c1 != c2
            return nothing
        end
    end
    return true, i1
end

@inline function format(io, d::Delim, dt, locale)
    print(io, d.d)
end

function _show_content(io::IO, d::Delim{<:AbstractChar, N}) where N
    if d.d in keys(CONVERSION_SPECIFIERS)
        for i = 1:N
            print(io, '\\', d.d)
        end
    else
        for i = 1:N
            print(io, d.d)
        end
    end
end

function _show_content(io::IO, d::Delim)
    for c in d.d
        if c in keys(CONVERSION_SPECIFIERS)
            print(io, '\\')
        end
        print(io, c)
    end
end

function Base.show(io::IO, d::Delim)
    print(io, "Delim(")
    _show_content(io, d)
    print(io, ")")
end

### DateFormat construction

abstract type DayOfWeekToken end # special addition to Period types

# Map conversion specifiers or character codes to tokens.
# Note: Allow addition of new character codes added by packages
const CONVERSION_SPECIFIERS = Dict{Char, Type}(
    'y' => Year,
    'Y' => Year,
    'm' => Month,
    'u' => Month,
    'U' => Month,
    'e' => DayOfWeekToken,
    'E' => DayOfWeekToken,
    'd' => Day,
    'H' => Hour,
    'M' => Minute,
    'S' => Second,
    's' => Millisecond,
)

# Default values are needed when a conversion specifier is used in a DateFormat for parsing
# and we have reached the end of the input string.
# Note: Allow `Any` value as a default to support extensibility
const CONVERSION_DEFAULTS = IdDict{Type, Any}(
    Year => Int64(1),
    Month => Int64(1),
    DayOfWeekToken => Int64(0),
    Day => Int64(1),
    Hour => Int64(0),
    Minute => Int64(0),
    Second => Int64(0),
    Millisecond => Int64(0),
    Microsecond => Int64(0),
    Nanosecond => Int64(0),
)

# Specifies the required fields in order to parse a TimeType
# Note: Allows for addition of new TimeTypes
const CONVERSION_TRANSLATIONS = IdDict{Type, Any}(
    Date => (Year, Month, Day),
    DateTime => (Year, Month, Day, Hour, Minute, Second, Millisecond),
    Time => (Hour, Minute, Second, Millisecond, Microsecond, Nanosecond),
)

"""
    DateFormat(format::AbstractString, locale="english") -> DateFormat

Construct a date formatting object that can be used for parsing date strings or
formatting a date object as a string. The following character codes can be used to construct the `format`
string:

| Code       | Matches   | Comment                                                      |
|:-----------|:----------|:-------------------------------------------------------------|
| `y`        | 1996, 96  | Returns year of 1996, 0096                                   |
| `Y`        | 1996, 96  | Returns year of 1996, 0096. Equivalent to `y`                |
| `m`        | 1, 01     | Matches 1 or 2-digit months                                  |
| `u`        | Jan       | Matches abbreviated months according to the `locale` keyword |
| `U`        | January   | Matches full month names according to the `locale` keyword   |
| `d`        | 1, 01     | Matches 1 or 2-digit days                                    |
| `H`        | 00        | Matches hours                                                |
| `M`        | 00        | Matches minutes                                              |
| `S`        | 00        | Matches seconds                                              |
| `s`        | .500      | Matches milliseconds                                         |
| `e`        | Mon, Tues | Matches abbreviated days of the week                         |
| `E`        | Monday    | Matches full name days of the week                           |
| `yyyymmdd` | 19960101  | Matches fixed-width year, month, and day                     |

Characters not listed above are normally treated as delimiters between date and time slots.
For example a `dt` string of "1996-01-15T00:00:00.0" would have a `format` string like
"y-m-dTH:M:S.s". If you need to use a code character as a delimiter you can escape it using
backslash. The date "1995y01m" would have the format "y\\ym\\m".

Creating a DateFormat object is expensive. Whenever possible, create it once and use it many times
or try the `dateformat""` string macro. Using this macro creates the DateFormat object once at
macro expansion time and reuses it later. see [`@dateformat_str`](@ref).

See [`DateTime`](@ref) and [`format`](@ref) for how to use a DateFormat object to parse and write Date strings
respectively.
"""
function DateFormat(f::AbstractString, locale::DateLocale=ENGLISH)
    tokens = AbstractDateToken[]
    prev = ()
    prev_offset = 1

    letters = String(collect(keys(CONVERSION_SPECIFIERS)))
    for m in eachmatch(Regex("(?<!\\\\)([\\Q$letters\\E])\\1*"), f)
        tran = replace(f[prev_offset:prevind(f, m.offset)], r"\\(.)" => s"\1")

        if !isempty(prev)
            letter, width = prev
            typ = CONVERSION_SPECIFIERS[letter]

            push!(tokens, DatePart{letter}(width, isempty(tran)))
        end

        if !isempty(tran)
            push!(tokens, Delim(length(tran) == 1 ? first(tran) : tran))
        end

        letter = f[m.offset]
        width = length(m.match)

        prev = (letter, width)
        prev_offset = m.offset + width
    end

    tran = replace(f[prev_offset:lastindex(f)], r"\\(.)" => s"\1")

    if !isempty(prev)
        letter, width = prev
        typ = CONVERSION_SPECIFIERS[letter]

        push!(tokens, DatePart{letter}(width, false))
    end

    if !isempty(tran)
        push!(tokens, Delim(length(tran) == 1 ? first(tran) : tran))
    end

    tokens_tuple = (tokens...,)
    return DateFormat{Symbol(f),typeof(tokens_tuple)}(tokens_tuple, locale)
end

function DateFormat(f::AbstractString, locale::AbstractString)
    DateFormat(f, LOCALES[locale])
end

function Base.show(io::IO, df::DateFormat)
    print(io, "dateformat\"")
    for t in df.tokens
        _show_content(io, t)
    end
    print(io, '"')
end
Base.Broadcast.broadcastable(x::DateFormat) = Ref(x)

"""
    dateformat"Y-m-d H:M:S"

Create a [`DateFormat`](@ref) object. Similar to `DateFormat("Y-m-d H:M:S")`
but creates the DateFormat object once during macro expansion.

See [`DateFormat`](@ref) for details about format specifiers.
"""
macro dateformat_str(str)
    DateFormat(str)
end

# Standard formats
const ISODateTimeFormat = DateFormat("yyyy-mm-dd\\THH:MM:SS.s")
const ISODateFormat = DateFormat("yyyy-mm-dd")
const ISOTimeFormat = DateFormat("HH:MM:SS.s")
const RFC1123Format = DateFormat("e, dd u yyyy HH:MM:SS")

default_format(::Type{DateTime}) = ISODateTimeFormat
default_format(::Type{Date}) = ISODateFormat
default_format(::Type{Time}) = ISOTimeFormat

### API

const Locale = Union{DateLocale, String}

"""
    DateTime(dt::AbstractString, format::AbstractString; locale="english") -> DateTime

Construct a `DateTime` by parsing the `dt` date time string following the
pattern given in the `format` string.

This method creates a `DateFormat` object each time it is called. If you are
parsing many date time strings of the same format, consider creating a
[`DateFormat`](@ref) object once and using that as the second argument instead.
"""
function DateTime(dt::AbstractString, format::AbstractString; locale::Locale=ENGLISH)
    return parse(DateTime, dt, DateFormat(format, locale))
end

"""
    DateTime(dt::AbstractString, df::DateFormat) -> DateTime

Construct a `DateTime` by parsing the `dt` date time string following the
pattern given in the [`DateFormat`](@ref) object. Similar to
`DateTime(::AbstractString, ::AbstractString)` but more efficient when
repeatedly parsing similarly formatted date time strings with a pre-created
`DateFormat` object.
"""
DateTime(dt::AbstractString, df::DateFormat=ISODateTimeFormat) = parse(DateTime, dt, df)

"""
    Date(d::AbstractString, format::AbstractString; locale="english") -> Date

Construct a `Date` by parsing the `d` date string following the pattern given
in the `format` string.

This method creates a `DateFormat` object each time it is called. If you are
parsing many date strings of the same format, consider creating a
[`DateFormat`](@ref) object once and using that as the second argument instead.
"""
function Date(d::AbstractString, format::AbstractString; locale::Locale=ENGLISH)
    parse(Date, d, DateFormat(format, locale))
end

"""
    Date(d::AbstractString, df::DateFormat) -> Date

Parse a date from a date string `d` using a `DateFormat` object `df`.
"""
Date(d::AbstractString, df::DateFormat=ISODateFormat) = parse(Date, d, df)

"""
    Time(t::AbstractString, format::AbstractString; locale="english") -> Time

Construct a `Time` by parsing the `t` time string following the pattern given
in the `format` string.

This method creates a `DateFormat` object each time it is called. If you are
parsing many time strings of the same format, consider creating a
[`DateFormat`](@ref) object once and using that as the second argument instead.
"""
function Time(t::AbstractString, format::AbstractString; locale::Locale=ENGLISH)
    parse(Time, t, DateFormat(format, locale))
end

"""
    Time(t::AbstractString, df::DateFormat) -> Time

Parse a time from a time string `t` using a `DateFormat` object `df`.
"""
Time(t::AbstractString, df::DateFormat=ISOTimeFormat) = parse(Time, t, df)

@generated function format(io::IO, dt::TimeType, fmt::DateFormat{<:Any,T}) where T
    N = fieldcount(T)
    quote
        ts = fmt.tokens
        loc = fmt.locale
        Base.@nexprs $N i -> format(io, ts[i], dt, loc)
    end
end

function format(dt::TimeType, fmt::DateFormat, bufsize=12)
    # preallocate to reduce resizing
    io = IOBuffer(Vector{UInt8}(undef, bufsize), read=true, write=true)
    format(io, dt, fmt)
    String(io.data[1:io.ptr - 1])
end

"""
    format(dt::TimeType, format::AbstractString; locale="english") -> AbstractString

Construct a string by using a `TimeType` object and applying the provided `format`. The
following character codes can be used to construct the `format` string:

| Code       | Examples  | Comment                                                      |
|:-----------|:----------|:-------------------------------------------------------------|
| `y`        | 6         | Numeric year with a fixed width                              |
| `Y`        | 1996      | Numeric year with a minimum width                            |
| `m`        | 1, 12     | Numeric month with a minimum width                           |
| `u`        | Jan       | Month name shortened to 3-chars according to the `locale`    |
| `U`        | January   | Full month name according to the `locale` keyword            |
| `d`        | 1, 31     | Day of the month with a minimum width                        |
| `H`        | 0, 23     | Hour (24-hour clock) with a minimum width                    |
| `M`        | 0, 59     | Minute with a minimum width                                  |
| `S`        | 0, 59     | Second with a minimum width                                  |
| `s`        | 000, 500  | Millisecond with a minimum width of 3                        |
| `e`        | Mon, Tue  | Abbreviated days of the week                                 |
| `E`        | Monday    | Full day of week name                                        |

The number of sequential code characters indicate the width of the code. A format of
`yyyy-mm` specifies that the code `y` should have a width of four while `m` a width of two.
Codes that yield numeric digits have an associated mode: fixed-width or minimum-width.
The fixed-width mode left-pads the value with zeros when it is shorter than the specified
width and truncates the value when longer. Minimum-width mode works the same as fixed-width
except that it does not truncate values longer than the width.

When creating a `format` you can use any non-code characters as a separator. For example to
generate the string "1996-01-15T00:00:00" you could use `format`: "yyyy-mm-ddTHH:MM:SS".
Note that if you need to use a code character as a literal you can use the escape character
backslash. The string "1996y01m" can be produced with the format "yyyy\\ymm\\m".
"""
function format(dt::TimeType, f::AbstractString; locale::Locale=ENGLISH)
    format(dt, DateFormat(f, locale))
end

# show

function Base.show(io::IO, dt::DateTime)
    if millisecond(dt) == 0
        format(io, dt, dateformat"YYYY-mm-dd\THH:MM:SS")
    else
        format(io, dt, dateformat"YYYY-mm-dd\THH:MM:SS.s")
    end
end

function Base.show(io::IO, dt::Date)
    format(io, dt, dateformat"YYYY-mm-dd")
end

function Base.string(dt::DateTime)
    if millisecond(dt) == 0
        format(dt, dateformat"YYYY-mm-dd\THH:MM:SS", 24)
    else
        format(dt, dateformat"YYYY-mm-dd\THH:MM:SS.s", 26)
    end
end

function Base.string(dt::Date)
    # don't use format - bypassing IOBuffer creation
    # saves a bit of time here.
    y,m,d = yearmonthday(value(dt))
    yy = y < 0 ? @sprintf("%05i", y) : lpad(y, 4, "0")
    mm = lpad(m, 2, "0")
    dd = lpad(d, 2, "0")
    return "$yy-$mm-$dd"
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
import gmsh

if length(ARGS) < 1
    println("Usage: julia explore.jl file.msh")
    exit(0)
end

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.open(ARGS[1])

# get all elementary entities in the model
entities = gmsh.model.getEntities()

for e in entities
    # get the mesh nodes for each elementary entity
    nodeTags, nodeCoords, nodeParams = gmsh.model.mesh.getNodes(e[1], e[2])
    # get the mesh elements for each elementary entity
    elemTypes, elemTags, elemNodeTags = gmsh.model.mesh.getElements(e[1], e[2])
    # report some statistics
    numElem = sum([length(t) for t in elemTags])
    println(length(nodeTags), " mesh nodes and ", numElem, " mesh elements ",
            "on entity ", e, " of type ", gmsh.model.getType(e[1], e[2]))
    partitions = gmsh.model.getPartitions(e[1], e[2])
    if length(partitions) > 0
        println(" - Partition tag(s): ", partitions, " - parent entity",
                gmsh.model.getParent(e[1], e[2]))
    end
    for t in elemTypes
        name, dim, order, numv, parv = gmsh.model.mesh.getElementProperties(t)
        println(" - Element type: ", name, ", order ", order)
        println("   with ", numv, " nodes in param coord: ", parv)
    end
end

gmsh.finalize()

max_stars_repo_name: Poofee/fastFEM

max_issues_repo_name: Christophe-Foyer/pyFEA

max_forks_repo_name: Christophe-Foyer/pyFEA


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using RELOG

@testset "geodb_query (2018-us-county)" begin
    region = RELOG.geodb_query("2018-us-county:17043")
    @test region.centroid.lat == 41.83956
    @test region.centroid.lon == -88.08857
    @test region.population == 922_921
end

# @testset "geodb_query (2018-us-zcta)" begin
#     region = RELOG.geodb_query("2018-us-zcta:60439")
#     @test region.centroid.lat == 41.68241
#     @test region.centroid.lon == -87.98954
# end

@testset "geodb_query (us-state)" begin
    region = RELOG.geodb_query("us-state:IL")
    @test region.centroid.lat == 39.73939
    @test region.centroid.lon == -89.50414
    @test region.population == 12_671_821
end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
using Test

@testset "LossFunctions" begin
    include("clusters_test.jl")
    include("loss_functions_test.jl")
    include("activation_functions_test.jl")
end

max_stars_repo_name: brianxie/automaton-edification

max_issues_repo_name: brianxie/automaton-edification

max_forks_repo_name: brianxie/automaton-edification


text:
# You enter a large cavern full of rare bioluminescent dumbo octopuses! They
# seem to not like the Christmas lights on your submarine, so you turn them off
# for now.
#
# There are 100 octopuses arranged neatly in a 10 by 10 grid. Each octopus
# slowly gains energy over time and flashes brightly for a moment when its
# energy is full. Although your lights are off, maybe you could navigate
# through the cave without disturbing the octopuses if you could predict when
# the flashes of light will happen.
#
# Each octopus has an energy level - your submarine can remotely measure the
# energy level of each octopus (your puzzle input). For example:
#
# 5483143223
# 2745854711
# 5264556173
# 6141336146
# 6357385478
# 4167524645
# 2176841721
# 6882881134
# 4846848554
# 5283751526
#
# The energy level of each octopus is a value between 0 and 9. Here, the
# top-left octopus has an energy level of 5, the bottom-right one has an energy
# level of 6, and so on.
#
# You can model the energy levels and flashes of light in steps. During a
# single step, the following occurs:
#
# - First, the energy level of each octopus increases by 1.
# - Then, any octopus with an energy level greater than 9 flashes. This
#   increases the energy level of all adjacent octopuses by 1, including
#   octopuses that are diagonally adjacent. If this causes an octopus to have
#   an energy level greater than 9, it also flashes. This process continues as
#   long as new octopuses keep having their energy level increased beyond 9.
#   (An octopus can only flash at most once per step.)
# - Finally, any octopus that flashed during this step has its energy level set
#   to 0, as it used all of its energy to flash.
#
# Adjacent flashes can cause an octopus to flash on a step even if it begins
# that step with very little energy. Consider the middle octopus with 1 energy
# in this situation:
#
# Before any steps:
# 11111
# 19991
# 19191
# 19991
# 11111
#
# After step 1:
# 34543
# 40004
# 50005
# 40004
# 34543
#
# After step 2:
# 45654
# 51115
# 61116
# 51115
# 45654
#
# An octopus is highlighted when it flashed during the given step.
#
# Here is how the larger example above progresses:
#
# Before any steps:
# 5483143223
# 2745854711
# 5264556173
# 6141336146
# 6357385478
# 4167524645
# 2176841721
# 6882881134
# 4846848554
# 5283751526
#
# After step 1:
# 6594254334
# 3856965822
# 6375667284
# 7252447257
# 7468496589
# 5278635756
# 3287952832
# 7993992245
# 5957959665
# 6394862637
#
# After step 2:
# 8807476555
# 5089087054
# 8597889608
# 8485769600
# 8700908800
# 6600088989
# 6800005943
# 0000007456
# 9000000876
# 8700006848
#
# After step 3:
# 0050900866
# 8500800575
# 9900000039
# 9700000041
# 9935080063
# 7712300000
# 7911250009
# 2211130000
# 0421125000
# 0021119000
#
# After step 4:
# 2263031977
# 0923031697
# 0032221150
# 0041111163
# 0076191174
# 0053411122
# 0042361120
# 5532241122
# 1532247211
# 1132230211
#
# After step 5:
# 4484144000
# 2044144000
# 2253333493
# 1152333274
# 1187303285
# 1164633233
# 1153472231
# 6643352233
# 2643358322
# 2243341322
#
# After step 6:
# 5595255111
# 3155255222
# 3364444605
# 2263444496
# 2298414396
# 2275744344
# 2264583342
# 7754463344
# 3754469433
# 3354452433
#
# After step 7:
# 6707366222
# 4377366333
# 4475555827
# 3496655709
# 3500625609
# 3509955566
# 3486694453
# 8865585555
# 4865580644
# 4465574644
#
# After step 8:
# 7818477333
# 5488477444
# 5697666949
# 4608766830
# 4734946730
# 4740097688
# 6900007564
# 0000009666
# 8000004755
# 6800007755
#
# After step 9:
# 9060000644
# 7800000976
# 6900000080
# 5840000082
# 5858000093
# 6962400000
# 8021250009
# 2221130009
# 9111128097
# 7911119976
#
# After step 10:
# 0481112976
# 0031112009
# 0041112504
# 0081111406
# 0099111306
# 0093511233
# 0442361130
# 5532252350
# 0532250600
# 0032240000
#
# After step 10, there have been a total of 204 flashes. Fast forwarding, here
# is the same configuration every 10 steps:
#
# After step 20:
# 3936556452
# 5686556806
# 4496555690
# 4448655580
# 4456865570
# 5680086577
# 7000009896
# 0000000344
# 6000000364
# 4600009543
#
# After step 30:
# 0643334118
# 4253334611
# 3374333458
# 2225333337
# 2229333338
# 2276733333
# 2754574565
# 5544458511
# 9444447111
# 7944446119
#
# After step 40:
# 6211111981
# 0421111119
# 0042111115
# 0003111115
# 0003111116
# 0065611111
# 0532351111
# 3322234597
# 2222222976
# 2222222762
#
# After step 50:
# 9655556447
# 4865556805
# 4486555690
# 4458655580
# 4574865570
# 5700086566
# 6000009887
# 8000000533
# 6800000633
# 5680000538
#
# After step 60:
# 2533334200
# 2743334640
# 2264333458
# 2225333337
# 2225333338
# 2287833333
# 3854573455
# 1854458611
# 1175447111
# 1115446111
#
# After step 70:
# 8211111164
# 0421111166
# 0042111114
# 0004211115
# 0000211116
# 0065611111
# 0532351111
# 7322235117
# 5722223475
# 4572222754
#
# After step 80:
# 1755555697
# 5965555609
# 4486555680
# 4458655580
# 4570865570
# 5700086566
# 7000008666
# 0000000990
# 0000000800
# 0000000000
#
# After step 90:
# 7433333522
# 2643333522
# 2264333458
# 2226433337
# 2222433338
# 2287833333
# 2854573333
# 4854458333
# 3387779333
# 3333333333
#
# After step 100:
# 0397666866
# 0749766918
# 0053976933
# 0004297822
# 0004229892
# 0053222877
# 0532222966
# 9322228966
# 7922286866
# 6789998766
#
# After 100 steps, there have been a total of 1656 flashes.
#
# Given the starting energy levels of the dumbo octopuses in your cavern,
# simulate 100 steps. How many total flashes are there after 100 steps?

# load grid
grid = hcat(map(eachline("input.txt")) do l
  parse.(Int8, split(l, ""))
end...)

# simulate 100 times
cnt = 0
for _ in 1:100
  # increment the energy of all octopuses
  grid .+= 1

  # find flashes
  all_flashes = falses(10, 10)
  flashes = grid .> 9
  while any(flashes)
    # increment the energy of the octopuses around flashing octopuses
    @views begin
      grid[1:9, 1:9] .+= flashes[2:10, 2:10]
      grid[1:10, 1:9] .+= flashes[1:10, 2:10]
      grid[2:10, 1:9] .+= flashes[1:9, 2:10]
      grid[2:10, 1:10] .+= flashes[1:9, 1:10]
      grid[2:10, 2:10] .+= flashes[1:9, 1:9]
      grid[1:10, 2:10] .+= flashes[1:10, 1:9]
      grid[1:9, 2:10] .+= flashes[2:10, 1:9]
      grid[1:9, 1:10] .+= flashes[2:10, 1:10]
    end

    # update all_flashes and see if there are any cascading flashes
    all_flashes = all_flashes .|| flashes
    flashes = grid .> 9 .&& .!all_flashes
  end

  # count the number of flashes in this round and set any flashing octopuses
  # back to zero
  global cnt += sum(all_flashes)
  grid[all_flashes] .= 0
end

println("Result: ", cnt)


text:
# workspace() # clear workspace by providing a clean workspace

# Used in many places as the minimum possible difference between two distinct
# numbers. Helps make code stable, tests predictable, etc.
kTolerance = Float64(1e-7)

typealias Feature Int64
typealias Label Int64
typealias NodeId Int64
typealias Value Float64
typealias Weight Float64

# An example consists of a vector of feature values, a label and a weight.
# Note that this is a dense feature representation; the value of every
# feature is contained in the vector, listed in a canonical order.
type Example
  values::Vector # Float64 Vector
  label::Label # Int64
  weight::Weight # Float64 probability (distribution)
  # construct function
  Example() = new()
  Example(values::Vector, # Float64 Vector
          label::Label, # Int64
          weight::Weight) = new(values, label, weight)
end

# a tree node
type Node
  examples::Array{Example, 1} # examples at this node
  split_feature::Int64 # Int64, split feature
  split_value::Float64 # Float64, split value
  left_child_id::Int64 # pointer to left child, if any
  right_child_id::Int64 # pointer to right child id, if any
  positive_weight::Float64 # Total weight of positive examples at this node.
  negative_weight::Float64 # total weight of negative examples at this node
  leaf::Bool # is this node a leaf?
  depth::Int64 # depth of the node in the tree. Root node has depth 0.
  Node() = new()
  Node(examples::Array{Example, 1}, # examples at this node
    split_feature::Int64, # Int64, split feature
    split_value::Float64, # Float64, split value
    left_child_id::Int64, # pointer to left child, if any
    right_child_id::Int64, # pointer to right child id, if any
    positive_weight::Float64, # Total weight of positive examples at this node.
    negative_weight::Float64, # total weight of negative examples at this node
    leaf::Bool, # is this node a leaf?
    depth::Int64) = new(examples, split_feature,split_value,
                        left_child_id, right_child_id,
                        positive_weight, negative_weight,
                        leaf, depth)
end

# a tree is a vector of nodes
typealias Tree Vector{Node}

type classifier
    weight::Weight
    tree::Tree
    classifier() = new()
    classifier(weight, tree) = new(weight, tree)
end

# a model is a vector of (weight, tree) pairs, i.e. a weighted combination of trees.
# by yy: for tuple is immutable
typealias Model Vector{classifier}

# a = Example([1.0,2.0,3.0],1,2)
# b = Node([a,a,a],2, 1.0, 2, 5, 6.0, 6.0, false,2)
# c = Node()
# d = Tree([b,b])
# e = Model([(1.2,d),(1.3,d)])


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m5_opt_rho1_9(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m5_opt_rho1_9!(A_copy)
end

@inline function exp_ps_m5_opt_rho1_9!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=8
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 T2k8
    # Computing Bb2 = x*I+x*A
    coeff1=0.36712607162139604
    coeff2=0.9948915464514331
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.36712607162139604
    coeff2=0.9948915464514331
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.06063859630763159
    coeff2=0.3498122134797794
    coeff3=1.2672705889590858
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.5288285666160359
    coeff2=1.252726870905159
    coeff3=0.08932442256452103
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.01410064578193085
    coeff2=-0.03083709071312985
    coeff3=0.33644565179790364
    coeff4=0.9865345056562885
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=1.0157619505241544
    coeff2=0.3045792025786807
    coeff3=0.04527440501144948
    coeff4=0.0014456607895543397
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-1.7630817992883847e-7
    coeff2=0.0013326770269791607
    coeff3=0.005896551938988777
    coeff4=-0.030425344123507926
    coeff5=1.0394518833422
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=5.021187762869422e-6
    coeff2=-1.3145878019370549e-7
    coeff3=2.157293609615063e-8
    coeff4=4.1404179600695967e-10
    coeff5=1.919036623692458e-11
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.003803985749916236
    coeff2=-0.026537261919315755
    coeff3=-0.1336685712514616
    coeff4=0.9611015405160707
    coeff5=0.7033098980139835
    coeff6=4.952503308762937e-5
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.05924467408855637
    coeff2=0.0038390703455692026
    coeff3=0.0005152684424798582
    coeff4=3.802841862018557e-5
    coeff5=2.051633742112383e-6
    coeff6=1.0000000003684573
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing T2k8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.953241101493861
    coeff2=0.7080362573925665
    coeff3=0.31966241179593696
    coeff4=-0.006733738250646286
    coeff5=-0.020217295003617748
    coeff6=-1.6649475207304e-5
    coeff7=0.9991068489851037
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    return memslots1 # Returning T2k8
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
"""
Samples f : R^d -> R at num points in box specified by lx and ux

INPUTS:
f: function being sampled
lx: length d array of lower bounds
ux: length d array of upper bounds

OUTPUTS:
res: length num array of sampled values
loc: num x d matrix of sample locations
"""
function sample_points(f, lx, ux; num = 10)
    @assert length(lx) == length(ux)
    d = length(lx)
    loc = rand(num, d)
    for j = 1:num
        for i = 1:d
            loc[j, i] = loc[j, i]*(ux[i] - lx[i]) + lx[i]
        end
    end
    res = zeros(num)
    for j  = 1:num
        res[j] = f(loc[j, :])
    end
    return res, loc
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg