text:
log_k_temp(k, T0) = T0 / log(k)

gaussian_sample_neighbor(s, T, D, ranges) = s + T .* randn(D)

function asa_sample_neighbor(s, T, D, ranges)
  u = rand(D)
  y = sign(u - 0.5) .* T .* ((1 + 1 ./ T).^abs(2*u-1) - 1)
  s .+ (y .* ranges)
end

function simulated_annealing(func,
                             s0,
                             numiterations = 10000,
                             sample_neighbor = asa_sample_neighbor,
                             temperature = log_k_temp,
                             ranges = ones(length(s0)))

  sbest = s = s0
  es = func(s)
  ebest = func(sbest)
  T0 = 100.0
  D = length(s)

  for k = 1:numiterations

    t = log_k_temp(k, T0)

    snew = sample_neighbor(s, t, D, ranges)
    enew = func(snew)

    delta_e = enew - es

    if delta_e < 0

      s = snew
      es = enew
      if es < ebest
        ebest = es
        sbest = s
        println("$(k), New best: $(ebest)")
      end

    else

      if rand() <= exp( -delta_e / t )
        s = snew
        es = enew
      end

    end

  end

  return sbest, ebest

end

function rosenbrock2d(x)
  return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

simulated_annealing(rosenbrock2d, randn(2), 100000)


text:
struct SpinState{T<:Number}
    dims :: Vector{Int}
    vals :: SparseVector{T}
    # function SpinState(dims::Vector{Int}, vals::SparseVector{T}) where {T}
    #     prod(dims) == size(vals,1) || error("size for spin state doesn't match", prod(dims), size(vals,1))
    #     new{T}(dims, vals)
    # end
end

# function SpinState(dims::Vector{Int}, vals::SparseVector{T}) where {T}
#     return SpinState{T}(dims, vals)
# end

@inline nspin(s::SpinState) = length(s.dims)
@inline dimension(s::SpinState) = prod(s.dims)

function kron(s1::SpinState{T}, s2::SpinState{T}) where {T}
    return SpinState([s1.dims; s2.dims], sparse(kron(s1.vals, s2.vals)))
end

function *(r::R, s::SpinState{T}) where {R<:Number, T}
    return SpinState(s.dims, r .* s.vals)
end

function +(s1::SpinState{T}, s2::SpinState{T}) where {T}
    @assert(s1.dims == s2.dims)
    return SpinState(s1.dims, s1.vals + s2.vals)
end

function -(s1::SpinState{T}, s2::SpinState{T}) where {T}
    @assert(s1.dims == s2.dims)
    return SpinState(s1.dims, s1.vals - s2.vals)
end

function permutespins(s::SpinState{T}, perm::Vector{Int}) where {T}
    @assert sort(perm) == collect(1:nspin(s))

    nzind = [conf2index(s, index2conf(s, index)[perm]) for index in s.vals.nzind]
    return SpinState(s.dims[perm], sparsevec(nzind, s.vals.nzval, dimension(s)))
end

function index2conf(s::SpinState, index::Int)
    res = zeros(Int, nspin(s))
    idx = index - 1
    for i in reverse(eachindex(s.dims))
        res[i] = idx % s.dims[i]
        idx = div(idx, s.dims[i])
    end
    return res
end

function conf2index(s::SpinState, c::Vector{Int})
    sum(c.* [reverse(cumprod(reverse(s.dims))); 1][2:end]) + 1
end

function confstring(conf::Vector{Int}, dims::Vector{Int};
                    unicode::Bool=true)
    symbs = [
        ["."],
        ["↓", "↑"],
        ["⇓", "⇔", "⇑"],
        ["⤋", "M", "W", "⤊"]
    ]
    res = ""
    for i in eachindex(conf)
        if dims[i] > length(symbs)
            res = res * string(conf[i])
        else
            res = res * symbs[dims[i]][conf[i]+1]
        end
    end
    return res
end

function sshow(s::SpinState; unicode::Bool=true)
    nzind = s.vals.nzind
    for n in eachindex(nzind)
        str = confstring(index2conf(s, nzind[n]), s.dims, unicode=unicode)
        println(str, ":   ", s.vals.nzval[n])
    end
    println()
end

ext: jl

max_stars_repo_name: amiragha/QSpin

max_issues_repo_name: amiragha/QSpin

max_forks_repo_name: amiragha/QSpin


text:
function Base.convert(T::Type{Vector{SVector{N2,T1}}}, M::Matrix) where {N2,T1}
    N1 = size(M)[1]
    @assert size(M)[2] == N2
    Mout = [szeros(T1, N2) for i = 1:N1]
    for i = 1:N1
        Mout[i] = convert(SVector{N2,T1}, M[i,:])
    end
    return Mout
end

Base.:*(u::LinearAlgebra.AdjointAbsVec, v::SVector{1,T}) where T = u * v[1]
Base.:*(u::AbstractVector, v::SVector{1,T}) where T = u * v[1]

Base.hcat(E::UniformScaling, v::SVector{N,T}) where {T, N} = hcat(SMatrix{N,N,T,N*N}(E), v)
Base.hcat(v::SVector{N,T}, E::UniformScaling) where {T, N} = hcat(v, SMatrix{N,N,T,N*N}(E))
Base.hcat(E::UniformScaling, v::Adjoint{T,SVector{N,T}}) where {T, N} = @error "not implemented" # shcat(E.λ, v)
Base.hcat(v::Adjoint{T,SVector{N,T}}, E::UniformScaling) where {T, N} = @error "not implemented" # shcat(v, E.λ)
Base.hcat(E::UniformScaling, A::SMatrix{N1,N2,T,N1N2}) where {T, N1, N2, N1N2} = hcat(SMatrix{N1,N1,T,N1*N1}(E), A)
Base.hcat(A::SMatrix{N1,N2,T,N1N2}, E::UniformScaling) where {T, N1, N2, N1N2} = hcat(A, SMatrix{N1,N1,T,N1*N1}(E))

Base.vcat(E::UniformScaling, v::SVector{N,T}) where {T, N} = svcat(E.λ, v)
Base.vcat(v::SVector{N,T}, E::UniformScaling) where {T, N} = svcat(v, E.λ)
Base.vcat(E::UniformScaling, v::Adjoint{T,SVector{N,T}}) where {T, N} = vcat(SMatrix{N,N,T,N*N}(E), v)
Base.vcat(v::Adjoint{T,SVector{N,T}}, E::UniformScaling) where {T, N} = vcat(v, SMatrix{N,N,T,N*N}(E))
Base.vcat(E::UniformScaling, A::SMatrix{N1,N2,T,N1N2}) where {T, N1, N2, N1N2} = vcat(SMatrix{N2,N2,T,N2*N2}(E), A)
Base.vcat(A::SMatrix{N1,N2,T,N1N2}, E::UniformScaling) where {T, N1, N2, N1N2} = vcat(A, SMatrix{N2,N2,T,N2*N2}(E))

Base.convert(::Type{SMatrix{N,N,T,N2}}, E::UniformScaling) where {T, N, N2} = SMatrix{N,N,T,N2}(E)

@inline svcat(a::T) where T = SA[a]
@inline svcat(a::StaticArray) = a
@inline svcat(a::T, b::T) where T = SA[a; b]
@inline svcat(a::StaticArray, b::StaticArray) = vcat(a, b)
@inline svcat(a::StaticArray{Tuple{N},T,1}, b::T) where {T,N} = vcat(a, SA[b])
@inline svcat(a::T, b::StaticArray{Tuple{N},T,1}) where {T,N} = vcat(SA[a], b)
@inline svcat(a, b, c...) = svcat(svcat(a, b), svcat(c...))

@inline szeros(::Type{T}, N) where T = @SVector zeros(T, N)
@inline szeros(N)= @SVector zeros(N)
@inline szeros(::Type{T}, N1, N2) where T = @SMatrix zeros(T, N1, N2)
@inline szeros(N1, N2)= @SMatrix zeros(N1, N2)

@inline sones(::Type{T}, N) where T = @SVector ones(T, N)
@inline sones(N)= @SVector ones(N)
@inline sones(::Type{T}, N1, N2) where T = @SMatrix ones(T, N1, N2)
@inline sones(N1, N2)= @SMatrix ones(N1, N2)

@inline srand(::Type{T}, N) where T = @SVector rand(T, N)
@inline srand(N)= @SVector rand(N)
@inline srand(::Type{T}, N1, N2) where T = @SMatrix rand(T, N1, N2)
@inline srand(N1, N2)= @SMatrix rand(N1, N2)

sisnan(a::StaticArray) = any(isnan.(a))

# To fix StaticArray bug
zerodimstaticadjoint(A) = A'
zerodimstaticadjoint(::SMatrix{0,N,T,0}) where {T,N} = SMatrix{N,0,T,0}()
zerodimstaticadjoint(::SMatrix{N,0,T,0}) where {T,N} = SMatrix{0,N,T,0}()
Base.:*(x::LinearAlgebra.AdjointAbsVec, A::SMatrix{0,N,T,0}) where {T,N} = (zerodimstaticadjoint(A)*x')'
Base.:*(A::SMatrix{N,0,T,0},x::SVector{0}) where {T,N} = szeros(T,N)
Base.:*(A::SMatrix{N,0,T,0},B::SMatrix{0,M,T,0}) where {T,N,M} = szeros(T,N,M)


text:
"""
    figsize(width_in_inch, height_width_ratio=0.68)

    Convert figure size in inches (which is often required by publishers) to that in points (which is the unit for CairoMakie). Note that 1 point in CairoMakie is equal to 1/72 inch.
"""
function figsize(width_in_inch, height_width_ratio=HWRATIO)
    width_in_point = floor(Int, 72width_in_inch)
    height_in_point = floor(Int, width_in_point*height_width_ratio)
    return width_in_point, height_in_point
end

savefig(f, fig) = save(f, fig, pt_per_unit=1)


text:
struct LevelSet
    interpolater::Any
    coefficients::Any
    mesh::Any
    dofspernode::Any
end

function LevelSet(distancefunc, mesh, basis::LagrangeTensorProductBasis)
    dofspernode = 1
    interpolater = InterpolatingPolynomial(1, basis)
    coeffs = interpolating_levelset_coefficients(distancefunc, mesh)
    LevelSet(interpolater, coeffs, mesh, dofspernode)
end

function LevelSet(
    distancefunc,
    mesh,
    basis::HermiteTensorProductBasis;
    numqp = 4,
)
    dim = dimension(mesh)
    quad = tensor_product_quadrature(dim, numqp)

    dofspernode = hermite_dofs_per_node(dim)
    interpolater = InterpolatingPolynomial(1, basis)
    coeffs = variational_levelset_coefficients(
        distancefunc,
        basis,
        quad,
        mesh,
        dofspernode,
    )
    LevelSet(interpolater, coeffs, mesh, dofspernode)
end

function Base.show(io::IO, levelset::LevelSet)
    p = order(levelset)
    numcoefficients = length(coefficients(levelset))
    dim = dimension(background_mesh(levelset))
    interpolatertype = typeof(basis(interpolater(levelset)))
    str =
        "LevelSet\n\tDimension : $dim\n\tOrder : $p" *
        "\n\tInterpolater Type : $interpolatertype"
    "\n\tTotal Num. Coefficients = $numcoefficients"
    print(io, str)
end

function order(poly::InterpolatingPolynomial)
    return PolynomialBasis.order(poly.basis)
end

function order(levelset::LevelSet)
    return order(interpolater(levelset))
end

function coefficients(levelset::LevelSet)
    return levelset.coefficients
end

function interpolating_levelset_coefficients(distancefunc, mesh)
    nodalcoordinates = nodal_coordinates(mesh)
    return distancefunc(nodalcoordinates)
end

function variational_levelset_coefficients(
    distancefunc,
    basis,
    quad,
    mesh,
    dofspernode,
)
    sysmatrix = SystemMatrix()
    sysrhs = SystemRHS()

    assemble_hermite_mass_matrix!(sysmatrix, basis, quad, mesh, dofspernode)
    assemble_hermite_linear_form!(
        sysrhs,
        distancefunc,
        basis,
        quad,
        mesh,
        dofspernode,
    )

    K = sparse_operator(sysmatrix, mesh, dofspernode)
    R = rhs_vector(sysrhs, mesh, dofspernode)

    coeffs = K \ R
    return coeffs
end

function coefficients(levelset::LevelSet, cellid)
    mesh = background_mesh(levelset)
    dofspernode = dofs_per_node(levelset)
    nodeids = nodal_connectivity(mesh, cellid)
    edofs = element_dofs(nodeids, dofspernode)
    return coefficients(levelset)[edofs]
end

function dofs_per_node(levelset::LevelSet)
    return levelset.dofspernode
end

function background_mesh(levelset::LevelSet)
    return levelset.mesh
end

function interpolater(levelset::LevelSet)
    return levelset.interpolater
end

function load_coefficients!(levelset::LevelSet, cellid)
    update!(interpolater(levelset), coefficients(levelset, cellid))
end

function update_coefficients!(levelset::LevelSet, cellid, coeffs)
    mesh = background_mesh(levelset)
    nodeids = nodal_connectivity(mesh, cellid)
    dofspernode = dofs_per_node(levelset)
    edofs = element_dofs(nodeids, dofspernode)
    levelset.coefficients[edofs] .= coeffs
end

function update_coefficients!(levelset::LevelSet, coeffs)
    levelset.coefficients[:] .= coeffs
end

function (levelset::LevelSet)(x)
    levelset.interpolater(x)
end


text:
using Plots

# Problem parameters
α = 0.3
δ = 0.1
β = 0.9

# Capital grid
kupper = 2
klower = 0.001
n = 10
kgrid = collect(range(klower, stop = kupper, length = n))

# Iteration parameters
tolerance = 0.001
imax = 1000

# Initialize values
vnew = zeros(length(kgrid)) # initial value function guess
v = vnew .+ 2*tolerance # initialize v in a way that ensures loop will start
global cartesianindex = Array{CartesianIndex{2}, n}
i = 1

while maximum(abs.(v - vnew)) > tolerance && i<=imax
    v = vnew;
    c = zeros(n, n);
    for i in 1:n
        for j in 1:n
            c[i,j] = kgrid[i]^α + (1-δ)*kgrid[i] - kgrid[j];
            if c[i,j] <= 0
                c[i,j] = 0;
            end
        end
        (vnew, cartesianindex) = findmax(log.(c) .+ β*v', dims = 2);
    end
    i += 1
end

scatter(kgrid, vnew, title = "v(k)")

# Policy function
kprimeindex = getindex.(cartesianindex, 2)
kprime = kgrid[kprimeindex]
scatter(kgrid, kprime, title = "k'(k)")


text:
using ChebyshevQuantum
using LinearAlgebra

function go(N)

    D = ChebyshevQuantum.Interp.getD(N);
    pts = ChebyshevQuantum.Interp.getCpts(N);
    #v = eigvals( -(D*D )[2:N,2:N]);
    v, _ = eigs( -(D*D )[2:N,2:N], nev=5, which=:SR, tol=1e-4, maxiter=30000);
    v / (pi^2/4)
end

function go_dense(N)

    D = ChebyshevQuantum.Interp.getD(N);
    pts = ChebyshevQuantum.Interp.getCpts(N);
    v = eigvals( -(D*D )[2:N,2:N]);
    #v, _ = eigs( -(D*D )[2:N,2:N], nev=10, which=:SR);
    v / (pi^2/4)
end


text:
###############################################################################
# Types
###############################################################################
export TraceItem

type TraceItem{T}
    location::String
    exprstr::String
    val::T
    lcount::Int
    ts::Float64 #time stamp
end
TraceItem{T}(location, exprstr, val::T, lcount) = TraceItem{T}(location, exprstr, val, lcount, time())

max_stars_repo_name: JobJob/RickTracy


text:
function generate_defs(model::Model)
  namelist()
  calcmod()
  modpara(num_sites(model))
  locspin(model.S2, num_sites(model))
  interall(model)
end

function calcmod()
  const io = open("calcmod.def", "w")
  println(io,
    """
    CalcType 0
    CalcModel   4
    OutputEigenVec 1 """)
  close(io)
end

function modpara(nsites::Integer)
  const io = open("modpara.def", "w")
  println(io,
  """
  --------------------
  Model_Parameters   0
  --------------------
  HPhi_Cal_Parameters
  --------------------
  CDataFileHead  zvo
  CParaFileHead  zqp
  --------------------
  Nsite          $nsites
  Lanczos_max    3000
  initial_iv     1
  nvec           1
  exct           1
  LanczosEps     14
  LanczosTarget  1""")
  close(io)
end

function locspin(S2::Integer, nsites::Integer)
  const io = open("locspin.def", "w")
  println(io,
  """
  === header
  NlocalSpin $nsites
  === reserved
  === reserved
  === end of header """)
  for i in 0:(nsites-1)
    println(io, i, " ", S2)
  end
  close(io)
end

function namelist()
  const io = open("namelist.def", "w")
  println(io,
  """
  CalcMod calcmod.def
  ModPara modpara.def
  LocSpin locspin.def
  InterAll interall.def """)
  close(io)
end

max_stars_repo_name: yomichi/HPhi-samples

max_issues_repo_name: yomichi/HPhi-samples

max_forks_repo_name: yomichi/HPhi-samples


text:
path = joinpath(@__DIR__, "precompile.csv")

result = String[]
open(path) do io
    for line in eachline(io)
        push!(result, split(line, '\t')[2][2:end-1])
    end
end

expressions = parse.(result)

using MacroTools
findmod2(v::Vector, mods = Set()) = (foreach(x-> findmod2(x, mods), v); mods)
function findmod2(expr, mods = Set())
    if @capture(expr, Type_{args__})
        findmod2(args, mods)
        findmod2(Type, mods)
    elseif @capture(expr, call_(args__) | (args__,) | Tuple{args__})
        findmod2(args, mods)
    elseif @capture(expr, getfield(Mod_, sym_))
        if isa(Mod, Symbol)
            push!(mods, Mod)
        else
            findmod2(Mod, mods)
        end
    elseif @capture(expr, Mod_.field_)
        if isa(Mod, Symbol)
            push!(mods, Mod)
        else
            findmod2(Mod, mods)
        end
    end
    mods
end

modules = findmod2(expressions)

:Pkg in names(Base)

open(joinpath(dirname(path), "precompie.jl"), "w") do io
    inbuilds = (:Logging, :Base, :Core, :Pkg, :Main, :Test)
    add_mods = filter(x-> !(x in inbuilds), modules)
    println(io, "Pkg.pkg\"add $(join(add_mods, " "))\"") # make sure packages are added
    println(io, "using $(join(modules, ", "))")
    for elem in expressions
        if @capture(elem, Tuple{f_, args__})
            println(io, "try")
            if isempty(args)
                println(io, "precompile($f, ())")
            else
                println(io, "precompile($f, ($(join(args, ", ")),))")
            end
            println(io, "catch; end")
        end
    end
end


text:
function postprocessing_integer(df::DataFrame,ndim::Int64,npoints::Int64)
    #res=DataFrame(undef, size(df))
    for i in 1:ndim
        ints = unique(df[:,i])
        for j in 1:length(ints)
            df[!,i][df[:,i] .== ints[j]] .= j
        end
    end
    return df
end

function postprocessing_cutpoints(df::DataFrame,ndim::Int64,npoints::Int64, cutpoints::Array{Float64})
    #res=DataFrame(undef, size(df))
    df=float.(df)
    for i in 1:ndim
        ints = unique(df[:,i])
        for j in 1:length(ints)
            df[!,i][df[:,i] .== ints[j]] .= .5*(sort(cutpoints[i,:])[j]-sort(cutpoints[i,:])[j+1])
        end
    end
    return df
end


text:
module Gruber
export gruber, niggly

using Crystals.Constants: default_tolerance
using Crystals.Utilities: cell_parameters
using Crystals.Structures: Crystal, volume
using Unitful: unit, ustrip, Quantity
using MicroLogging
using ArgCheck

function no_opt_change_test(new, last)
    const m_new = 16e0 * new;
    const difference = new - last;
    const m_new_plus_diff = m_new + difference;
    const m_new_plus_diff_minus_m_new = m_new_plus_diff - m_new;
    m_new_plus_diff_minus_m_new != 0;
end

function def_test(real::Real; tolerance::Real=default_tolerance)
    real ≥ tolerance && return [0, 1]
    real > -tolerance && return [1, 0]
    return [0, 0]
end

function def_test(args; tolerance=default_tolerance)
    result = [0, 0]
    for u in args
        result += def_test(u, tolerance=tolerance)
    end
    result
end

function def_gt_0(args...; tolerance=default_tolerance)
    z₀, positive = def_test(args; tolerance=tolerance)
    positive == 3 || (z₀ == 0 && positive == 1)
end

function n1_action(params::Vector, rinv::Matrix)
    rinv[:, :] = rinv * [0 -1 0; -1 0 0; 0 0 -1]
    params[1], params[2] = params[2], params[1]
    params[4], params[5] = params[5], params[4]
end

function n2_action(params::Vector, rinv::Matrix)
    rinv[:, :] = rinv * [-1 0 0; 0 0 -1; 0 -1 0]
    params[2], params[3] = params[3], params[2]
    params[5], params[6] = params[6], params[5]
end

function n3_action(params::Vector, rinv::Matrix; tolerance=default_tolerance)
    const i = params[4] ≤ -tolerance ? -1 : 1
    const j = params[5] ≤ -tolerance ? -1 : 1
    const k = params[6] ≤ -tolerance ? -1 : 1
    rinv[:, :] = rinv * [i 0 0; 0 j 0; 0 0 k]
    params[4:end] = abs.(params[4:end])
end

function n4_action(params::Vector, rinv::Matrix; tolerance=default_tolerance)
    const i = params[4] ≥ tolerance ? -1 : 1
    const j = params[5] ≥ tolerance ? -1 : 1
    const k = params[6] ≥ tolerance ? -1 : 1
    update = diagm([i, j, k])
    if i * j * k < 0
        if k == 1 && params[6] > -tolerance
            update[3, 3] = -1
        elseif j == 1 && params[5] > -tolerance
            update[3, 3] = -1
        elseif i == 1 && params[4] > -tolerance
            update[1, 1] = -1
        elseif i * j * k == -1
            error("Internal error")
        end
    end
    rinv[:, :] = rinv * update
    params[4:end] = -abs.(params[4:end])
end

function n5_action(params::Vector, rinv::Matrix; tolerance=default_tolerance)
    const pos_or_neg = params[4] > tolerance ? -1 : 1
    rinv[:, :] = rinv * [1 0 0; 0 1 pos_or_neg; 0 0 1]
    params[3] += params[2] + pos_or_neg * params[4]
    params[4] += 2pos_or_neg * params[2]
    params[5] += pos_or_neg * params[6]
end

function n6_action(params::Vector, rinv::Matrix; tolerance=default_tolerance)
    const pos_or_neg = params[5] > tolerance ? -1 : 1
    rinv[:, :] = rinv * [1 0 pos_or_neg; 0 1 0; 0 0 1]
    params[3] += params[1] + pos_or_neg * params[5]
    params[4] += pos_or_neg * params[6]
    params[5] += 2pos_or_neg * params[1]
end

function n7_action(params::Vector, rinv::Matrix; tolerance=default_tolerance)
    const pos_or_neg = params[6] > tolerance ? -1 : 1
    rinv[:, :] = rinv * [1 pos_or_neg 0; 0 1 0; 0 0 1]
    params[2] += params[1] + pos_or_neg * params[6]
    params[4] += pos_or_neg * params[5]
    params[6] += 2pos_or_neg * params[1]
end

function n8_action(params::Vector, rinv::Matrix)
    rinv[:, :] = rinv * [1 0 1; 0 1 1; 0 0 1]
    params[3] += sum(params[1:2]) + sum(params[4:end])
    params[4] += 2params[2] + params[6]
    params[5] += 2params[1] + params[6]
end

"""
    gruber(cell::Matrix;
           tolerance::Real=default_tolerance, itermax::Unsigned=50,
           max_no_change::Unsigned=10)

Determines Gruber cell of an input cell.

The Gruber cell is an optimal parameterization of a lattice, e.g. shortest
cell-vectors and angles closest to 90 degrees. The new cell is in the same basis
as the origin cell: no rotation has been incurred. The cell parameters are
uniquely determined, even though the cell itself is not (certain symmetry
operations leaving the structure unchanged may yield a more recognizable cell).
If you want a unique Cartesian cell (in a different Cartesian basis), use
the `niggly` algorithm.

# Arguments
* `cell::Matrix`: the input lattice cell-vectors. Cannot be singular.
* `itermax::Integer`: maximum number of iterations before bailing out
* `tolerance::Number`: tolerance parameter when comparing real numbers
* `max_no_change::Integer`: Maximum number of times to go through algorithm
  without changes before bailing out
"""
function gruber{T <: Number}(cell::Matrix{T};
                             tolerance::Real=default_tolerance, itermax::Integer=50,
                             max_no_change::Integer=10)
    @argcheck size(cell, 1) == size(cell, 2)
    @argcheck volume(cell) > tolerance
    if itermax ≤ 0
        itermax = typemax(itermax)
    end
    if max_no_change ≤ 0
        max_no_change = typemax(max_no_change)
    end

    const ε = tolerance
    const metric = transpose(cell) * cell
    params = vcat(diag(metric), [2metric[2, 3], 2metric[1, 3], 2metric[1, 2]])
    rinv = eye(size(metric, 1))
    no_change, previous = 0, -params[1:3]
    iteration::Int = 0
    for iteration in 1:itermax
        condition0 =
            (a, b, d, e) -> a ≥ b + ε || (abs(a - b) < ε && abs(d) ≥ abs(e) + ε)

        condition0(params[1], params[2], params[4], params[5]) &&
            n1_action(params, rinv)
        condition0(params[2], params[3], params[5], params[6]) &&
            (n2_action(params, rinv); continue)

        if def_gt_0(params[4:end]...; tolerance=ε)
            n3_action(params, rinv; tolerance=ε)
        else
            n4_action(params, rinv; tolerance=ε)
            if all(abs.(previous .- params[1:3]) .< ε)
                no_change += 1
            else
                no_change = 0
            end
            no_change < max_no_change || break
        end

        condition1 =
            (d, b, e, f) -> abs(d) ≥ b + ε ||
            (abs(d - b) < ε && 2e ≤ f - ε) ||
            (abs(d + b) < ε && f ≤ -ε)
        condition1(params[4], params[2], params[5], params[6]) &&
            (n5_action(params, rinv; tolerance=ε); continue)
        condition1(params[5], params[1], params[4], params[6]) &&
            (n6_action(params, rinv; tolerance=ε); continue)
        condition1(params[6], params[1], params[4], params[5]) &&
            (n7_action(params, rinv; tolerance=ε); continue)

        const sum_no_c = sum(params[1:2]) + sum(params[4:end])
        (
            sum_no_c ≤ -ε ||
            (abs(params[6] - params[1]) < ε && (2params[4] ≤ params[5] - ε)) ||
            (abs(sum_no_c) < ε  && 2params[1] + 2params[5] + params[6] ≥ ε)
        ) || break
        n8_action(params, rinv)
    end
    iteration == itermax && error("Reached end of iteration without converging")
    cell * rinv
end

function gruber{T, D, U}(cell::Matrix{Quantity{T, D, U}};
                         tolerance::Real=default_tolerance, itermax::Integer=50,
                         max_no_change::Integer=10)
    gruber(ustrip.(cell);
           tolerance=tolerance,
           itermax=itermax,
           max_no_change=max_no_change) * unit(Quantity{T, D, U})
end

"""
    niggly(cell::Matrix; kwargs...)

Determines a unique Cartesian cell equivalent to the input, with the shortest
possible vectors and squarest angles. For an explanation of the parameters, see
`gruber`. In practice, this function computes the cell-parameters of a `gruber` cell and
then reconstructs the cell matrix. Hence, the result should be quite unique for any lattice
representation, including any rotation of the underlying Cartesian coordinates.
"""
niggly(cell::Matrix, kwargs...) = cell_parameters(cell_parameters(gruber(cell; kwargs...)))

end


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.
##########################################################################

# returns grid indices for given location
function locationToGridIndex(map::Map, grid::Grid, location::Location)
	# fails if location.x == map.xMin or location.y == map.yMin, can fix this later if needed
	ix = Int(ceil((location.x - map.xMin) / map.xRange * grid.nx))
	iy = Int(ceil((location.y - map.yMin) / map.yRange * grid.ny))
	if location.x == map.xMin ix = 1 end
	if location.y == map.yMin iy = 1 end
	@assert(1 <= ix && ix <= grid.nx && 1 <= iy && iy <= grid.ny)
	return ix, iy
end

# place nodes in rectangles of grid
# can specify whether nodes require offRoadAccess = true, in order to be added
function gridPlaceNodes!(map::Map, grid::Grid, nodes::Vector{Node};
	offRoadAccessRequired::Bool = true)
	for node in nodes
		if node.offRoadAccess || !offRoadAccessRequired
			(ix, iy) = locationToGridIndex(map, grid, node.location)
			push!(grid.rects[ix,iy].nodeIndices, node.index)
		end
	end
end

# searches grid for nearest node to given location,
# starting in grid rectangle of location, then
# extending search outward in a rectangle, until
# borders of search rectangle are further from given
# location than the nearest node found
function findNearestNode(map::Map, grid::Grid, nodes::Vector{Node}, location::Location)

	# indices of grid rectangle for given location
	(ix, iy) = locationToGridIndex(map, grid, location)

	# initialise grid search rectangle
	gsr = grid.searchRect
	gsr.ixSearched[1] = ix
	gsr.ixSearched[2] = ix
	gsr.iySearched[1] = iy
	gsr.iySearched[2] = iy
	gsr.ixSearch[1] = ix
	gsr.ixSearch[2] = ix
	gsr.iySearch[1] = iy
	gsr.iySearch[2] = iy
	gsr.xDist[1] = (location.x - (map.xMin + grid.xRange * (ix-1))) * map.xScale
	gsr.xDist[2] = grid.xRangeDist - gsr.xDist[1]
	gsr.yDist[1] = (location.y - (map.yMin + grid.yRange * (iy-1))) * map.yScale
	gsr.yDist[2] = grid.yRangeDist - gsr.yDist[1]

	nearestNode = nullIndex
	bestDist = Inf
	node = nullIndex # init
	dist = nullDist # init

	found = false
	skipSearch = false
	while !found

		if !skipSearch
			for i = gsr.ixSearch[1] : gsr.ixSearch[2]
				for j = gsr.iySearch[1] : gsr.iySearch[2]
					# find nearest node in grid rect
					nearestNodeInGrid = nullIndex
					bestSqrDist = Inf
					sqrDist = 0.0 # init
					for nodeIndex in grid.rects[i,j].nodeIndices
						sqrDist = squareDist(map, nodes[nodeIndex].location, location)
						if sqrDist < bestSqrDist
							bestSqrDist = sqrDist
							nearestNodeInGrid = nodeIndex
						end
					end
					dist = sqrt(bestSqrDist)

					if dist < bestDist
						bestDist = dist
						nearestNode = nearestNodeInGrid
					end
				end
			end
		end
		nearestBorderDist = min(gsr.xDist[1],gsr.xDist[2],gsr.yDist[1],gsr.yDist[2])
		if bestDist <= nearestBorderDist
			found = true
		else
			# extend rectangle search border
			skipSearch = false
			if gsr.xDist[1] == nearestBorderDist
				if gsr.ixSearched[1] > 1 # extend border
					gsr.ixSearched[1] -= 1
					gsr.xDist[1] += grid.xRangeDist
					gsr.ixSearch[1] = gsr.ixSearched[1]
					gsr.ixSearch[2] = gsr.ixSearched[1]
					gsr.iySearch[1] = gsr.iySearched[1]
					gsr.iySearch[2] = gsr.iySearched[2]
				else # cannot extend border
					gsr.xDist[1] = Inf
					skipSearch = true
				end
			elseif gsr.xDist[2] == nearestBorderDist
				if gsr.ixSearched[2] < grid.nx # extend border
					gsr.ixSearched[2] += 1
					gsr.xDist[2] += grid.xRangeDist
					gsr.ixSearch[1] = gsr.ixSearched[2]
					gsr.ixSearch[2] = gsr.ixSearched[2]
					gsr.iySearch[1] = gsr.iySearched[1]
					gsr.iySearch[2] = gsr.iySearched[2]
				else # cannot extend border
					gsr.xDist[2] = Inf
					skipSearch = true
				end
			elseif gsr.yDist[1] == nearestBorderDist
				if gsr.iySearched[1] > 1 # extend border
					gsr.iySearched[1] -= 1
					gsr.yDist[1] += grid.yRangeDist
					gsr.ixSearch[1] = gsr.ixSearched[1]
					gsr.ixSearch[2] = gsr.ixSearched[2]
					gsr.iySearch[1] = gsr.iySearched[1]
					gsr.iySearch[2] = gsr.iySearched[1]
				else # cannot extend border
					gsr.yDist[1] = Inf
					skipSearch = true
				end
			elseif gsr.yDist[2] == nearestBorderDist
				if gsr.iySearched[2] < grid.ny # extend border
					gsr.iySearched[2] += 1
					gsr.yDist[2] += grid.yRangeDist
					gsr.ixSearch[1] = gsr.ixSearched[1]
					gsr.ixSearch[2] = gsr.ixSearched[2]
					gsr.iySearch[1] = gsr.iySearched[2]
					gsr.iySearch[2] = gsr.iySearched[2]
				else # cannot extend border
					gsr.yDist[2] = Inf
					skipSearch = true
				end
			end
		end

	end

	return nearestNode, bestDist # (gsr.ixSearched[2] - gsr.ixSearched[1] + 1)*(gsr.iySearched[2] - gsr.iySearched[1] + 1)
end


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

# Generate decode table.
const BASE64_CODE_END = 0x40
const BASE64_CODE_PAD = 0x41
const BASE64_CODE_IGN = 0x42
const BASE64_DECODE = fill(BASE64_CODE_IGN, 256)
for (i, c) in enumerate(BASE64_ENCODE)
    BASE64_DECODE[Int(c)+1] = UInt8(i - 1)
end
BASE64_DECODE[Int(encodepadding())+1] = BASE64_CODE_PAD
decode(x::UInt8) = @inbounds return BASE64_DECODE[x + 1]

"""
    Base64DecodePipe(istream)

Return a new read-only I/O stream, which decodes base64-encoded data read from
`istream`.

# Examples
```jldoctest
julia> io = IOBuffer();

julia> iob64_decode = Base64DecodePipe(io);

julia> write(io, "SGVsbG8h")
8

julia> seekstart(io);

julia> String(read(iob64_decode))
"Hello!"
```
"""
struct Base64DecodePipe <: IO
    io::IO
    buffer::Buffer
    rest::Vector{UInt8}

    function Base64DecodePipe(io::IO)
        buffer = Buffer(512)
        return new(io, buffer, UInt8[])
    end
end

Base.isreadable(pipe::Base64DecodePipe) = !isempty(pipe.rest) || isreadable(pipe.io)
Base.iswritable(::Base64DecodePipe) = false

function Base.unsafe_read(pipe::Base64DecodePipe, ptr::Ptr{UInt8}, n::UInt)
    p = read_until_end(pipe, ptr, n)
    if p < ptr + n
        throw(EOFError())
    end
    return nothing
end

# Read and decode as much data as possible.
function read_until_end(pipe::Base64DecodePipe, ptr::Ptr{UInt8}, n::UInt)
    p = ptr
    p_end = ptr + n
    while !isempty(pipe.rest) && p < p_end
        unsafe_store!(p, popfirst!(pipe.rest))
        p += 1
    end

    buffer = pipe.buffer
    i = 0
    b1 = b2 = b3 = b4 = BASE64_CODE_IGN
    while true
        if b1 < 0x40 && b2 < 0x40 && b3 < 0x40 && b4 < 0x40 && p + 2 < p_end
            # fast path to decode
            unsafe_store!(p    , b1 << 2 | b2 >> 4)
            unsafe_store!(p + 1, b2 << 4 | b3 >> 2)
            unsafe_store!(p + 2, b3 << 6 | b4     )
            p += 3
        else
            i, p, ended = decode_slow(b1, b2, b3, b4, buffer, i, pipe.io, p, p_end - p, pipe.rest)
            if ended
                break
            end
        end
        if p < p_end
            if i + 4 ≤ lastindex(buffer)
                b1 = decode(buffer[i+1])
                b2 = decode(buffer[i+2])
                b3 = decode(buffer[i+3])
                b4 = decode(buffer[i+4])
                i += 4
            else
                consumed!(buffer, i)
                read_to_buffer(pipe.io, buffer)
                i = 0
                b1 = b2 = b3 = b4 = BASE64_CODE_IGN
            end
        else
            break
        end
    end
    consumed!(buffer, i)

    return p
end

function Base.read(pipe::Base64DecodePipe, ::Type{UInt8})
    if isempty(pipe.rest)
        unsafe_read(pipe, convert(Ptr{UInt8}, C_NULL), 0)
        if isempty(pipe.rest)
            throw(EOFError())
        end
    end
    return popfirst!(pipe.rest)
end

function Base.readbytes!(pipe::Base64DecodePipe, data::AbstractVector{UInt8}, nb::Integer=length(data))
    require_one_based_indexing(data)
    filled::Int = 0
    while filled < nb && !eof(pipe)
        if length(data) == filled
            resize!(data, min(length(data) * 2, nb))
        end
        p = pointer(data, filled + 1)
        p_end = read_until_end(pipe, p, UInt(min(length(data), nb) - filled))
        filled += p_end - p
    end
    resize!(data, filled)
    return filled
end

Base.eof(pipe::Base64DecodePipe) = isempty(pipe.rest) && eof(pipe.io)::Bool
Base.close(pipe::Base64DecodePipe) = nothing

# Decode data from (b1, b2, b3, b5, buffer, input) into (ptr, rest).
function decode_slow(b1, b2, b3, b4, buffer, i, input, ptr, n, rest)
    # Skip ignore code.
    while true
        if b1 == BASE64_CODE_IGN
            b1, b2, b3 = b2, b3, b4
        elseif b2 == BASE64_CODE_IGN
            b2, b3 = b3, b4
        elseif b3 == BASE64_CODE_IGN
            b3 = b4
        elseif b4 == BASE64_CODE_IGN
            # pass
        else
            break
        end
        if i + 1 ≤ lastindex(buffer)
            b4 = decode(buffer[i+=1])
        elseif !eof(input)
            b4 = decode(read(input, UInt8))
        else
            b4 = BASE64_CODE_END
            break
        end
    end

    # Check the decoded quadruplet.
    k = 0
    if b1 < 0x40 && b2 < 0x40 && b3 < 0x40 && b4 < 0x40
        k = 3
    elseif b1 < 0x40 && b2 < 0x40 && b3 < 0x40 && b4 == BASE64_CODE_PAD
        b4 = 0x00
        k = 2
    elseif b1 < 0x40 && b2 < 0x40 && b3 == b4 == BASE64_CODE_PAD
        b3 = b4 = 0x00
        k = 1
    elseif b1 == b2 == b3 == BASE64_CODE_IGN && b4 == BASE64_CODE_END
        b1 = b2 = b3 = b4 = 0x00
    else
        throw(ArgumentError("malformed base64 sequence"))
    end

    # Write output.
    p::Ptr{UInt8} = ptr
    p_end = ptr + n
    function output(b)
        if p < p_end
            unsafe_store!(p, b)
            p += 1
        else
            push!(rest, b)
        end
    end
    k ≥ 1 && output(b1 << 2 | b2 >> 4)
    k ≥ 2 && output(b2 << 4 | b3 >> 2)
    k ≥ 3 && output(b3 << 6 | b4     )

    return i, p, k == 0
end

"""
    base64decode(string)

Decode the base64-encoded `string` and returns a `Vector{UInt8}` of the decoded
bytes.

See also [`base64encode`](@ref).

# Examples
```jldoctest
julia> b = base64decode("SGVsbG8h")
6-element Vector{UInt8}:
 0x48
 0x65
 0x6c
 0x6c
 0x6f
 0x21

julia> String(b)
"Hello!"
```
"""
function base64decode(s)
    b = IOBuffer(s)
    try
        return read(Base64DecodePipe(b))
    finally
        close(b)
    end
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
#############################################################################
#Copyright (c) 2016-2019 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
#############################################################################

VERSION >= v"1.0.0" && __precompile__()

module Spectra

using StatsBase
using Statistics
using Random
using Statistics
using SparseArrays
using LinearAlgebra
using PyPlot
using LsqFit
using PyCall
using Dierckx
#using NMF
using Polynomials

# For PyCall modules
const preprocessing = PyNULL()
const model_selection = PyNULL()
const decomposition = PyNULL()
const kernel_ridge = PyNULL()
const svm = PyNULL()
const linear_model = PyNULL()
const gaussian_process = PyNULL()
const pygcvspl = PyNULL()
const signal = PyNULL()
const rampy = PyNULL()

function __init__()

	copy!(preprocessing, pyimport_conda("sklearn.preprocessing", "scikit-learn"))
	copy!(model_selection, pyimport_conda("sklearn.model_selection", "scikit-learn"))
	copy!(decomposition, pyimport_conda("sklearn.decomposition", "scikit-learn"))
	copy!(kernel_ridge, pyimport_conda("sklearn.kernel_ridge", "scikit-learn"))
	copy!(svm, pyimport_conda("sklearn.svm","scikit-learn"))
	copy!(gaussian_process, pyimport_conda("sklearn.gaussian_process", "scikit-learn"))
	copy!(linear_model, pyimport_conda("sklearn.linear_model", "scikit-learn"))
	copy!(pygcvspl, pyimport("gcvspline"))
	copy!(rampy, pyimport("rampy"))
	copy!(signal, pyimport_conda("scipy.signal", "scipy"))

end

include("integrale.jl")
include("functions.jl")
include("baseline.jl")
include("bootstrap.jl")
include("tlcorrection.jl")
include("deprecated.jl")
include("ml_regressor.jl")
include("peakmeasurement.jl")

#From integrale.jl
export trapz, bandarea

#From functions.jl
export poly, polyfit, gaussiennes, lorentziennes, pseudovoigts, pearson7, normal_dist, xshift_inversion, xshift_direct,xshift_correction, smooth, flipsp, resample

#From baseline.jl
export baseline

#From bootstrap
export bootsample, bootperf

#From tlcorrection
export tlcorrection

#From ml_regressor
export mlregressor

# From peakmeasurement
export peakmeas

end # module


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 126b1cca-2476-11eb-087b-474174c4c1d0
begin
begin
begin
	import Pkg
	Pkg.add("Plots")
end

using Plots
end
	plotly()
end

# ╔═╡ e0937cd4-239d-11eb-141d-1be13a185dca
md"# _Ideal Gas Laws!_
All gases are affected by several variables: Volume, temperature, amount, and pressure. Together, they all affect how a gas behaves. But, before we look at those variables, lets learn the basic rules of gases."

# ╔═╡ 17827eb0-239f-11eb-3266-2d71e0aedd9c
md"## Kinetic Theory of Gases
1. **Gases are mostly empty space** ~ the total number of molecules in any space is much smaller compared to the volume of that space
2. **Gas molecules are in constant, chaotic motion** ~ the molecules constantly collide with eachother and the container walls
3. **Collisions are elastic** ~ there is no energy loss during collisions nor are there attractive forces that keep the molecules together
4. **Gas pressure is caused by collisions of molecules against the walls of a container** ~ Pressure increases as the energy and frequency of the collisions increase.
"

# ╔═╡ e51f5ba4-239f-11eb-0d73-cbb2f0b727b7
md" Now that we understand how gases behave generally, lets look at specific variables that impact gases"

# ╔═╡ 07010ba2-23a0-11eb-2656-b12946848a01
md"## Relationship: Volume and Temperature
Before we look into the relationship between the two, it is important to understand several distinctions about each variable.

**Volume of a Gas equals Volume of a Container** : Gases will expand uniformly to fill any containers that it was placed in which means the volume of the container equals the volume of the gas."

# ╔═╡ 3b34ea3a-23a1-11eb-31b5-3171d41bc984
md"**Temperature in Kelvin**: Whenever discussing the behavior of gas, it is very important to use kelvin temperature. To relate this back to celsius, lets look at a conversion equation."

# ╔═╡ 1b88d926-23a3-11eb-0d2e-7f8a2fa4eb1b
md"""Move the slider to change the degrees in Celsius and to see how that affects the temperature in kelvin: """

# ╔═╡ 9f4d4afa-23a4-11eb-082b-cf05aa408069
 @bind C html"<input type=range>"

# ╔═╡ a34f586e-23a4-11eb-3ac3-13713ac16268
md"$C degrees Celsius equates to"

# ╔═╡ 65a88f26-23a4-11eb-02af-8998dbbbcf7f
TemperatureK= C + 273.15

# ╔═╡ c3f552c0-23a5-11eb-2174-33ebd9de1cff
begin

	md" The equation to calculate kelvin from celsius is T(k) = C + 273.15"
end

# ╔═╡ cc464fea-23a5-11eb-3e9e-87f62b2a4faf
md" Now that we know the basics of Volume and temperature, lets examine their relationship!"

# ╔═╡ e450fb58-23a5-11eb-2caa-639830646520
md".............................................................................................................................................................................................................."

# ╔═╡ 466d76b8-2475-11eb-32e4-19a694079c2c
md"**VOLUME is DIRECTLY PROPORTIONAL to absolute TEMPERATURE**"

# ╔═╡ a50225d4-2475-11eb-3c89-9b50dc329236
md" Volume does not depend on the type of gas you have because all the gases behave the same. So, for all gases, as the temperature increases, so does the volume."

# ╔═╡ 4a55c8d8-2476-11eb-07d3-9b01d5e70a1e
Temperature = [0,50,100,150,200,250,300]

# ╔═╡ 008ce65e-2490-11eb-3b8a-f91d3cb8e452
Volume = [0,1,2,3,4,5,6]

# ╔═╡ bc1d2f24-2476-11eb-0c81-17c4da48b190

# ╔═╡ 68c8e54c-2490-11eb-3d2b-dd73b9af11cf
plot(Temperature, Volume,legend=false, xlabel="Temperature(K)", ylabel="Volume(L)", title="Correlation between Temperature and Volume")

# ╔═╡ be952f68-2476-11eb-049b-1beb8f592f70
md"Since the relationship is directly proportional, we can use this principle and apply it into the equation `V1/T1 = V2/T2`. This equation is apart of a larger equation which we will cover later on in this section"

# ╔═╡ bf35d8e4-2476-11eb-078b-3b6b1f5fa04e
md"Now, lets move on to a different relationship: the relationship between amount and volume."

# ╔═╡ e5a02fd2-2475-11eb-3f1e-152a7a3a43ce
md"## Relationship: Volume and Amount
We understand a lot about volume, but before we can look into the relationship between amount and volume, we need to understand amount and the general assumptions we make about the relationship between volume and amount."

# ╔═╡ dfb0faac-24e3-11eb-350d-d7b5fed50b49
md"First, amount usually expressed in **moles**, the international unit for measuring the amount of a substance. This can be converted into grams(which depends on the substance you are looking for) or moles(which has the conversion factor of 6.022*10^23). Second, when we examine the relationship between volume and amount, it is assumed that the volume of the container is NOT fixed. This is important because if the volume is fixed, other variables are changed and you are no longer examining the relationship between amount and volume. Now lets talk about their relationship!"

# ╔═╡ dc0e97cc-24e5-11eb-215e-5dd5e41f7946
md"**VOLUME is DIRECTLY PROPORTIONAL to AMOUNT**"

# ╔═╡ f5e00744-24e5-11eb-3760-a75ea6b16b9c
md" The relationship between volume and amount is similar to volume and temperature. Assuming that the volume is not fixed, the more moles of a substance you have, the greater the volume. This is because as you increase the number of particles, the volume has to expand to keep everything constant."

# ╔═╡ 89333174-24e6-11eb-36d3-71f27b566e4f
begin
	plotly()
	NumberofMoles = [0.0,0.05,0.1,0.15,0.2,0.25,0.3]
	Volume2 = [0.0,1.0,2.0,3.0,4.0,5.0,6.0]
end

# ╔═╡ ee6eac6c-24e6-11eb-2df9-1fb31e1f2c96
plot(NumberofMoles, Volume2,legend=false, xlabel="Number of Moles", ylabel="Volume(L)", title="Correlation between Amount and Volume")

# ╔═╡ d504768c-27b9-11eb-2a2c-e5ccce269690
md"Since the relationship is directly proportional, we can use this principle and apply it into the equation `V1/N1 = V2/N2`. This equation is also apart of a larger equation which we will cover later on in this section"

# ╔═╡ e82d8afa-27b9-11eb-0ad1-417d616337dd
md"Now lets discuss the next component: Pressure and Temperature"

# ╔═╡ 075ec916-27ba-11eb-3ea2-a516b703aa90
md"## Relationship: Pressure and Temperature
Before we look at the relationship between temperature and pressure, lets understand some basic facts about pressure"

# ╔═╡ 2f624a46-27ba-11eb-205a-3b78744f49af
md" Pressure is defined as **force per unit area**. In the context of gases, as mentioned in the kinetic theory of gases, pressure in gases is caused by collisions into particles or the walls of a container. There are several ways to measure pressure, and each way depends on the conditions in which the pressure is measured. A **manometer** measures pressure in a closed system. This means that it will measure the pressure of a sealed container. A **barometer** measures pressure in a open system, like the atmosphere."

# ╔═╡ 82beda50-27bb-11eb-11af-c75e38f6fa2a
md"Pressure forms based on the number of collisions in particles. The greater the number of collisions, the higher the pressure. Since more collisions are a direct result of increase in temperature, we can say that the higher the temperature, the higher the pressure.(This is assuming that volume remains constant)"

# ╔═╡ 59e62f4e-2f62-11eb-097e-f32238a77ec4
Temperature2 = [0,50,100,150,200,250,300]

# ╔═╡ 9762122c-2f62-11eb-03f7-ff2dc854fe40
Pressure = [0,1,2,3,4,5,6]

# ╔═╡ a00adb0c-2f62-11eb-2dec-23a060be15cc
plot(Temperature2, Pressure,legend=false, xlabel="Temperature(Kelvin)", ylabel="Pressure(atm)", title="Correlation between Pressure and Temperature")

# ╔═╡ d10de92c-2f62-11eb-15c6-9b30bca82516
md"Pressure is directly proportional to Temperature so when the temperature increases, pressure also increases. The equation for this is `P1/T1 = P2/T2`. Now lets look at the final relationship before putting everything together"

# ╔═╡ bd83ed76-2f66-11eb-08df-338b5985c137
md" ## Relationship: Volume and Pressure
This relationship is different from the rest because the relationship between volume and temperature is inversly proportional to pressure. As the volume increases, the pressure decreases and as the volume decreases, the pressure increases. This happens because as the volume increase, the molecules have a lot more room to move around so less collisions happen which means less pressure. On the flipside, when volume decreases, there is less space for the molecules to move around so they collide more often and those collisions result in an increase in pressure."

# ╔═╡ 8d3afe88-2f67-11eb-308e-ad7440ee187c
Volume3 = [11,6,3,2,1.5,1.25,1]

# ╔═╡ 28b7c798-2f67-11eb-1ee7-ff083e19acaf
plot(Pressure, Volume3 ,legend=false, xlabel="Pressure(atm)", ylabel="Volume(L)", title="Correlation between Volume and Pressure")

# ╔═╡ d3129b2a-2f67-11eb-15b7-c9bac218038c
md"Now we understand each component, Volume, pressure, temperature, and moles, and their relationships with eachother, we can put all the factors together to form one equation."

# ╔═╡ f8bbb5ba-2f67-11eb-32ef-bbc3f87fa010
md" ## Ideal Gas Law Formula
The Ideal Gas formula is
$\ce{PV=nRT}$ , where P stands for pressure, V stands for volume, n stands for moles, and T stands for temperature. The R is a special value called the gas constant, which stays the same no matter the calculation. R = 0.0821. Understanding this, you can now solve for each individual variable. "

# ╔═╡ Cell order:
# ╟─e0937cd4-239d-11eb-141d-1be13a185dca
# ╟─17827eb0-239f-11eb-3266-2d71e0aedd9c
# ╟─e51f5ba4-239f-11eb-0d73-cbb2f0b727b7
# ╟─07010ba2-23a0-11eb-2656-b12946848a01
# ╟─3b34ea3a-23a1-11eb-31b5-3171d41bc984
# ╠═1b88d926-23a3-11eb-0d2e-7f8a2fa4eb1b
# ╠═9f4d4afa-23a4-11eb-082b-cf05aa408069
# ╠═a34f586e-23a4-11eb-3ac3-13713ac16268
# ╠═65a88f26-23a4-11eb-02af-8998dbbbcf7f
# ╟─c3f552c0-23a5-11eb-2174-33ebd9de1cff
# ╟─cc464fea-23a5-11eb-3e9e-87f62b2a4faf
# ╟─e450fb58-23a5-11eb-2caa-639830646520
# ╟─466d76b8-2475-11eb-32e4-19a694079c2c
# ╟─a50225d4-2475-11eb-3c89-9b50dc329236
# ╟─126b1cca-2476-11eb-087b-474174c4c1d0
# ╠═4a55c8d8-2476-11eb-07d3-9b01d5e70a1e
# ╠═008ce65e-2490-11eb-3b8a-f91d3cb8e452
# ╟─bc1d2f24-2476-11eb-0c81-17c4da48b190
# ╟─68c8e54c-2490-11eb-3d2b-dd73b9af11cf
# ╟─be952f68-2476-11eb-049b-1beb8f592f70
# ╟─bf35d8e4-2476-11eb-078b-3b6b1f5fa04e
# ╟─e5a02fd2-2475-11eb-3f1e-152a7a3a43ce
# ╟─dfb0faac-24e3-11eb-350d-d7b5fed50b49
# ╟─dc0e97cc-24e5-11eb-215e-5dd5e41f7946
# ╟─f5e00744-24e5-11eb-3760-a75ea6b16b9c
# ╟─89333174-24e6-11eb-36d3-71f27b566e4f
# ╟─ee6eac6c-24e6-11eb-2df9-1fb31e1f2c96
# ╟─d504768c-27b9-11eb-2a2c-e5ccce269690
# ╟─e82d8afa-27b9-11eb-0ad1-417d616337dd
# ╟─075ec916-27ba-11eb-3ea2-a516b703aa90
# ╟─2f624a46-27ba-11eb-205a-3b78744f49af
# ╟─82beda50-27bb-11eb-11af-c75e38f6fa2a
# ╟─59e62f4e-2f62-11eb-097e-f32238a77ec4
# ╟─9762122c-2f62-11eb-03f7-ff2dc854fe40
# ╟─a00adb0c-2f62-11eb-2dec-23a060be15cc
# ╟─d10de92c-2f62-11eb-15c6-9b30bca82516
# ╟─bd83ed76-2f66-11eb-08df-338b5985c137
# ╟─8d3afe88-2f67-11eb-308e-ad7440ee187c
# ╟─28b7c798-2f67-11eb-1ee7-ff083e19acaf
# ╟─d3129b2a-2f67-11eb-15b7-c9bac218038c
# ╟─f8bbb5ba-2f67-11eb-32ef-bbc3f87fa010


text:
using Documenter, UnitfulLsqFit

makedocs(; sitename="UnitfulLsqFit")
deploydocs(; repo="github.com/gustaphe/UnitfulLsqFit.jl.git", devbranch="main")


text:
# # Auto-tuning Hyperparameters

#md # [![](https://img.shields.io/badge/show-github-579ACA.svg)](@__REPO_ROOT_URL__/docs/src/examples/autotuning-ridge.jl)

# This example shows how to learn a hyperparameter in Ridge Regression using a gradient descent routine.
# Let the regularized regression problem be formulated as:

# ```math
# \begin{equation}
# \min_{w} \quad \frac{1}{2nd} \sum_{i=1}^{n} (w^T x_{i} - y_i)^2 + \frac{\alpha}{2d} \| w \|_2^2
# \end{equation}
# ```

# where
# - `x`, `y` are the data points
# - `w` are the learned weights
# - `α` is the hyperparameter acting on regularization.

# The main optimization model will be formulated with JuMP.
# Using the gradient of the optimal weights with respect to the regularization parameters
# computed with DiffOpt, we can perform a gradient descent on top of the inner model
# to minimize the test loss.

# This tutorial uses the following packages

using JuMP     # The mathematical programming modelling language
import DiffOpt # JuMP extension for differentiable optimization
import Ipopt    # Optimization solver that handles quadratic programs
import Plots   # Graphing tool
import LinearAlgebra: norm, dot
import Random

# ## Generating a noisy regression dataset

Random.seed!(42)

N = 100
D = 20
noise = 5

w_real = 10 * randn(D)
X = 10 * randn(N, D)
y = X * w_real + noise * randn(N)
l = N ÷ 2  # test train split

X_train = X[1:l, :]
X_test  = X[l+1:N, :]
y_train = y[1:l]
y_test  = y[l+1:N];

# ## Defining the regression problem

# We implement the regularized regression problem as a function taking the problem data,
# building a JuMP model and solving it.

function fit_ridge(model, X, y, α)
    JuMP.empty!(model)
    set_silent(model)
    N, D = size(X)
    @variable(model, w[1:D])
    @expression(model, err_term, X * w - y)
    @objective(
        model,
        Min,
        dot(err_term, err_term) / (2 * N * D) + α * dot(w, w) / (2 * D),
    )
    optimize!(model)
    @assert termination_status(model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED, MOI.ALMOST_LOCALLY_SOLVED]
    return w
end

# We can solve the problem for several values of α
# to visualize the effect of regularization on the testing and training loss.

αs = 0.00:0.01:0.50
mse_test = Float64[]
mse_train = Float64[]
model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
(Ntest, D) = size(X_test)
(Ntrain, D) = size(X_train)
for α in αs
    w = fit_ridge(model, X_train, y_train, α)
    ŵ = value.(w)
    ŷ_test = X_test * ŵ
    ŷ_train = X_train * ŵ
    push!(mse_test, norm(ŷ_test - y_test)^2 / (2 * Ntest * D))
    push!(mse_train, norm(ŷ_train - y_train)^2 / (2 * Ntrain * D))
end

# Visualize the Mean Score Error metric

Plots.plot(
    αs, mse_test ./ sum(mse_test),
    label="MSE test", xaxis = "α", yaxis="MSE", legend=(0.8, 0.2),
    width=3,
)
Plots.plot!(
    αs, mse_train ./ sum(mse_train),
    label="MSE train",
    linestyle=:dash,
    width=3,
)
Plots.title!("Normalized MSE on training and testing sets")

# ## Leveraging differentiable optimization: computing the derivative of the solution

# Using DiffOpt, we can compute `∂w_i/∂α`, the derivative of the learned solution `̂w`
# w.r.t. the regularization parameter.

function compute_dw_dα(model, w)
    D = length(w)
    dw_dα = zeros(D)
    MOI.set(
        model,
        DiffOpt.ForwardObjectiveFunction(),
        dot(w, w)  / (2 * D),
    )
    DiffOpt.forward_differentiate!(model)
    for i in 1:D
        dw_dα[i] = MOI.get(
            model,
            DiffOpt.ForwardVariablePrimal(),
            w[i],
        )
    end
    return dw_dα
end

# Using `∂w_i/∂α` computed with `compute_dw_dα`,
# we can compute the derivative of the test loss w.r.t. the parameter α
# by composing derivatives.

function d_testloss_dα(model, X_test, y_test, w, ŵ)
    N, D = size(X_test)
    dw_dα = compute_dw_dα(model, w)
    err_term = X_test * ŵ - y_test
    return sum(eachindex(err_term)) do i
        dot(X_test[i,:], dw_dα) * err_term[i]
    end / (N * D)
end

# We can define a meta-optimizer function performing gradient descent
# on the test loss w.r.t. the regularization parameter.

function descent(α0, max_iters=100; fixed_step = 0.01, grad_tol=1e-3)
    α_s = Float64[]
    ∂α_s = Float64[]
    test_loss = Float64[]
    α = α0
    N, D = size(X_test)
    model = Model(() -> DiffOpt.diff_optimizer(Ipopt.Optimizer))
    for iter in 1:max_iters
        w = fit_ridge(model, X_train, y_train, α)
        ŵ = value.(w)
        err_term = X_test * ŵ - y_test
        ∂α = d_testloss_dα(model, X_test, y_test, w, ŵ)
        push!(α_s, α)
        push!(∂α_s, ∂α)
        push!(test_loss, norm(err_term)^2 / (2 * N * D))
        α -= fixed_step * ∂α
        if abs(∂α) ≤ grad_tol
            break
        end
    end
    return α_s, ∂α_s, test_loss
end

ᾱ, ∂ᾱ, msē = descent(0.10, 500)
iters = 1:length(ᾱ);

# Visualize gradient descent and convergence

Plots.plot(
    αs, mse_test,
    label="MSE test", xaxis = ("α"),
    legend=:topleft,
    width=2,
)
Plots.plot!(
    ᾱ, msē,
    label="learned α", width = 5,
    style=:dot,
)
Plots.title!("Regularizer learning")

# Visualize the convergence of α to its optimal value

Plots.plot(
    iters, ᾱ, label = nothing, color = :blue,
    xaxis = ("Iterations"), legend=:bottom,
    title = "Convergence of α"
)

# Visualize the convergence of the objective function

Plots.plot(
    iters, msē, label = nothing, color = :red,
    xaxis = ("Iterations"), legend=:bottom,
    title = "Convergence of MSE"
)

# Visualize the convergence of the derivative to zero

Plots.plot(
    iters, ∂ᾱ, label = nothing, color = :green,
    xaxis = ("Iterations"), legend=:bottom,
    title = "Convergence of ∂α"
)


text:
using AutoGrad: getval

# default feature set used in parser
FEATS=["s1c","s1v","s1p","s1A","s1a","s1B","s1b",
       "s1rL", # "s1rc","s1rv","s1rp",
       "s0c","s0v","s0p","s0A","s0B","s0a","s0b","s0d",
       "s0rL", # "s0rc","s0rv","s0rp",
       "n0lL", # "n0lc","n0lv","n0lp",
       "n0c","n0v","n0p","n0A","n0a",
       "n1c","n1v","n1p",
       ]
FTYPE=Float32                   # floating point type
GPUFEATURES=false               # Whether to compute gpu features or not
MAXWORD=32                      # truncate long words at this length. length("counterintelligence")=19
MAXSENT=64                      # skip longer sentences during training
MINSENT=2                       # skip shorter sentences during training

function main(o)
    println("opts=",[(k,v) for (k,v) in o]...)

    @msg o[:loadfile]
    d = load(o[:loadfile])

    if o[:seed] > 0
        srand(o[:seed])
    end

    if length(o[:dropout])==1
        o[:dropout]=ntuple(i->o[:dropout][1],2)
    end
    if length(o[:report])==1
        o[:report]=ntuple(i->o[:report][1], length(o[:datafiles]))
    end

    # assign transition based Parser type Abstraction from parser.jl
    o[:arctype] = (o[:arctype] == nothing ? get(d,"arctype",ArcHybridR1) : eval(parse(o[:arctype])))

    if haskey(d,"arctype") && o[:arctype] != d["arctype"]
        error("ArcType mismatch")
    end

    o[:feats] = (o[:feats] == nothing ? get(d,"feats",FEATS) : eval(parse(o[:feats])))
    if haskey(d,"feats") && o[:feats] != d["feats"]
        error("Feats mismatch")
    end
    # we specify three embedding dims for postag,deprel and count/distance for empty models.
    # the word and context embed sizes are given by chmodel dims.
    if isempty(o[:embed]) && !haskey(d,"parserv")
        o[:embed] = (128,32,16) # default embedding sizes
    elseif !isempty(o[:embed]) && haskey(d,"parserv") # check compat
        if !(length(d["postagv"][1])==o[:embed][1] &&
             length(d["deprelv"][1])==o[:embed][2] &&
             length(d["lcountv"][1])==length(d["rcountv"][1])==length(d["distancev"][1])==o[:embed][3])
            error("Embed mismatch")
        end
    end

    vocab = create_vocab(d)
    wmodel = makewmodel(d)
    corpora = []
    lfunc = (VERSION >= v"0.6.0" ? loadcorpus_v6 : loadcorpus)
    for f in o[:datafiles]; @msg f
        c = lfunc(f,vocab)
        push!(corpora,c)
    end

    cc = vcat(corpora...)
    info("Context and word embeddings are being generated")
    ppl = fillvecs!(wmodel, cc, vocab)
    unk = unkrate(cc)
    @msg "perplexity=$ppl unkrate=$(unk[1]/unk[2])"

    @msg :initmodel
    (pmodel,optim) = makepmodel(d,o,corpora[1][1])
    # add-hoc savefile function
    save1(file)=(@msg file; savefile(file, vocab, wmodel, pmodel, optim, o[:arctype], o[:feats]))

    # report function that is going to be used during training
    function report(epoch,beamsize=o[:beamsize])
        las_vals = zeros(length(corpora))
        for i=1:length(corpora)
            if o[:report][i] != 0 || (o[:bestfile] != nothing && i==length(corpora))
                las_vals[i] = beamtest(model=pmodel,
                                  corpus=corpora[i],
                                  vocab=vocab,
                                  arctype=o[:arctype],
                                  feats=o[:feats],
                                  beamsize=beamsize,
                                  batchsize=o[:batchsize])
            end
        end
        println((:epoch,epoch,:beam,beamsize,:las,las_vals...))
        return las_vals[end]
    end
    @msg :parsing
    bestlas = report(0,1)
    bestepoch = 0

    if o[:otrain]>0
        @msg :otrain
    end
    free_KnetArray()
    for epoch=1:o[:otrain]
        oracletrain(model=pmodel,
                    optim=optim,
                    corpus=corpora[1],
                    vocab=vocab,
                    arctype=o[:arctype],
                    feats=o[:feats],
                    batchsize=o[:batchsize],
                    pdrop=o[:dropout])
        currlas = report("oracle$epoch",1)
        if currlas > bestlas
            bestlas = currlas
            bestepoch = epoch
            if o[:bestfile] != nothing
                save1(o[:bestfile])
            end
        end
        #if 5 < bestepoch < epoch - 5
        #    break
        #end
    end

    # savemodel
    if o[:savefile] != nothing
        save1(o[:savefile])
    end

    # output parse for corpora[1]
    if o[:output] != nothing    # TODO: parse all data files?
        @msg o[:output]
        if corpora[1][1].parse == nothing
            @msg :parsing
            beamtest(model=pmodel,corpus=corpora[1],vocab=vocab,arctype=o[:arctype],feats=o[:feats],beamsize=o[:beamsize],batchsize=o[:batchsize])
        end
        writeconllu(corpora[1], o[:datafiles][1], o[:output])
    end
    @msg :done

end

function beamtest(;model=_model, corpus=_corpus, vocab=_vocab, arctype=ArcHybridR1, feats=FEATS, beamsize=4, batchsize=128) # large batchsize does not slow down beamtest
    for s in corpus
        s.parse = nothing
    end

    sentbatches = minibatch(corpus,batchsize)
    for sentences in sentbatches
        beamloss(model, sentences, vocab, arctype, feats, beamsize; earlystop=false)
    end
    las(corpus)
end

function beamloss(pmodel, sentences, vocab, arctype, feats, beamsize; earlystop=true, steps=nothing, pdrop=(0,0))
    parsers = parsers0 = map(arctype, sentences)
    beamends = beamends0 = collect(1:length(parsers)) # marks end column of each beam, initially one parser per beam
    pcosts  = pcosts0  = zeros(Int, length(sentences))
    pscores = pscores0 = zeros(FTYPE, length(sentences))
    totalloss = stepcount = 0
    featmodel,mlpmodel = splitmodel(pmodel)

    while length(beamends) > 0
        # features (vcat) are faster on cpu, mlp is faster on gpu
        fmatrix = features(parsers, feats, featmodel) # nfeat x nparser
        if GPUFEATURES && gpu()>=0 #GPU
            @assert isa(getval(fmatrix),KnetArray{FTYPE,2})
        else #CPU
            @assert isa(getval(fmatrix),Array{FTYPE,2})
            if gpu()>=0; fmatrix = KnetArray(fmatrix); end
        end
        cscores = Array(mlp(mlpmodel, fmatrix; pdrop=pdrop)) # candidate scores: nmove x nparser
        cscores = cscores .+ pscores' # candidate cumulative scores
        parsers,pscores,pcosts,beamends,loss = nextbeam(parsers, cscores, pcosts, beamends, beamsize; earlystop=earlystop)
        totalloss += loss
        stepcount += 1
    end
    # emptyvecs!(sentences)       # if we don't empty, gc cannot clear these vectors; empty if finetuning wvecs
    if steps != nothing
        steps[1] += stepcount
        steps[2] += length(sentences[1])*2-2
    end
    return totalloss / length(sentences)
end

function nextbeam(parsers, mscores, pcosts, beamends, beamsize; earlystop=true)
    n = beamsize * length(beamends) + 1
    newparsers, newscores, newcosts, newbeamends, loss = Array(Any,n),Array(Int,n),Array(Int,n),Int[],0.0
    nmoves,nparsers = size(mscores)                     # mscores[m,p] is the score of move m for parser p
    mcosts = Array(Any, nparsers)                       # mcosts[p][m] will be the cost vector for parser[p],move[m] if needed
    n = p0 = 0
    for p1 in beamends                                  # parsers[p0+1:p1], mscores[:,p0+1:p1] is the current beam belonging to a common sentence
        s0,s1 = 1 + nmoves*p0, nmoves*p1                # mscores[s0:s1] are the linear indices for mscores[:,p0:p1]
        nsave = n                                       # newparsers,newscores,newcosts[nsave+1:n] will the new beam for this sentence
        ngold = 0                                       # ngold, if nonzero, will be the index of the gold path in beam
        sorted = sortperm(getval(mscores)[s0:s1], rev=true)
        for isorted in sorted
            linidx = isorted + s0 - 1
            (move,parent) = ind2sub(size(mscores), linidx) # find cartesian index of the next best score
            parser = parsers[parent]

            # skip illegal move
            if !moveok(parser,move)
                continue
            end
            if earlystop && !isassigned(mcosts, parent) # not every parent may have children, avoid unnecessary movecosts
                mcosts[parent] = movecosts(parser, parser.sentence.head, parser.sentence.deprel)
            end
            n += 1
            newparsers[n] = copy(parser); move!(newparsers[n], move)
            newscores[n] = linidx
            #TEST: no newcosts during test
            if earlystop
                newcosts[n] = pcosts[parent] + mcosts[parent][move]
                if newcosts[n] == 0
                    if ngold==0
                        ngold=n
                    else
                        # @msg("multiple gold moves for $(parser.sentence)")
                    end
                end
            end
            if n-nsave == beamsize; break; end
        end
        if n == nsave
            if parsers[p1].nword == 1                   # single word sentences give us no moves
                s = parsers[p1].sentence
                s.parse = parsers[p1]
            else
                error("No legal moves?")                # otherwise this is suspicious
            end
        #TEST: there will be no ngold during test
        elseif earlystop && ngold == 0                  # gold path fell out of beam, early stop
            gindex = goldindex(parsers,pcosts,mcosts,(p0+1):p1)
            if gindex != 0
                newscores[n+1] = gindex
                loss = loss - mscores[gindex] + logsumexp2(mscores, newscores[nsave+1:n+1])
            end
            n = nsave
        elseif endofparse(newparsers[n])                # all parsers in beam have finished, gold among them if earlystop
            #TEST: there will be no ngold during test, cannot return beamloss, just return the highest scoring parse, no need for normalization
            if earlystop
                gindex = newscores[ngold]
                loss = loss - mscores[gindex] + logsumexp2(mscores, newscores[nsave+1:n])
            end
            s = newparsers[n].sentence
            # @assert s == newparsers[nsave+1].sentence
            # @assert mscores[newscores[nsave+1]] >= mscores[newscores[n]] "s[$(1+nsave)]=$(mscores[newscores[nsave+1]]) s[$n]=$(mscores[newscores[n]])"
            s.parse = newparsers[nsave+1]
            n = nsave                                   # do not add finished parsers to new beam
        else                                            # all good keep going
            push!(newbeamends, n)
        end
        p0 = p1
    end
    return newparsers[1:n], mscores[newscores[1:n]], newcosts[1:n], newbeamends, loss
end

function logsumexp2(a, r)
    amax = getval(a)[r[1]]
    log(sum(exp(a[r] - amax))) + amax
end

function goldindex(parsers, pcosts, mcosts, beamrange)
    parent = findfirst(view(pcosts, beamrange), 0)
    if parent == 0
        error("cannot find gold parent in $beamrange")
    end
    parent += first(beamrange) - 1
    if !isassigned(mcosts, parent)
        p = parsers[parent]
        mcosts[parent] = movecosts(p, p.sentence.head, p.sentence.deprel)
    end
    move = findfirst(mcosts[parent],0)
    if move == 0
        return 0
    else
        msize = (parsers[1].nmove,length(parsers))
        return sub2ind(msize,move,parent)
    end
end

endofparse(p)=(p.sptr == 1 && p.wptr > p.nword)

function minibatch(corpus, batchsize; maxlen=typemax(Int), minlen=1, shuf=false)
    data = Any[]
    sorted = sort(corpus, by=length)
    i1 = findfirst(x->(length(x) >= minlen), sorted)
    if i1==0; error("No sentences >= $minlen"); end
    i2 = findlast(x->(length(x) <= maxlen), sorted)
    if i2==0; error("No sentences <= $maxlen"); end
    for i in i1:batchsize:i2
        j = min(i2, i+batchsize-1)
        push!(data, sorted[i:j])
    end
    if shuf
        data=shuffle(data)
    end
    return data
end

function las(corpus)
    nword = ncorr = 0
    for s in corpus
        p = s.parse
        nword += length(s)
        ncorr += sum((s.head .== p.head) & (s.deprel .== p.deprel))
    end
    ncorr / nword
end

function splitmodel(pmodel)
    # optimization: do all getindex operations outside, otherwise each getindex creates a new node
    # TODO: fix this in general
    mlpmodel = Any[]
    mlptemp = parserv(pmodel)
    for i=1:length(mlptemp); push!(mlpmodel, mlptemp[i]); end
    featmodel = Array(Any,5)
    for k in 1:5 # (:postag,:deprel,:lcount,:rcount,:distance)
        featmodel[k] = Any[]
        pmodel_k = pmodel[k]
        for i in 1:length(pmodel_k)
            push!(featmodel[k], pmodel_k[i])
        end
    end
    return (featmodel,mlpmodel)
end

function oracletrain(;model=_model, optim=_optim, corpus=_corpus, vocab=_vocab, arctype=ArcHybridR1, feats=FEATS, batchsize=16, maxiter=typemax(Int), pdrop=(0,0))
    sentbatches = minibatch(corpus,batchsize; maxlen=MAXSENT, minlen=MINSENT, shuf=true)
    nsent = sum(map(length,sentbatches)); nsent0 = length(corpus)
    nword = sum(map(length,vcat(sentbatches...))); nword0 = sum(map(length,corpus))
    @msg("nsent=$nsent/$nsent0 nword=$nword/$nword0")
    nwords = StopWatch()
    losses = Any[0,0,0]
    niter = 0
    @time for sentences in sentbatches
        grads = oraclegrad(model, sentences, vocab, arctype, feats; losses=losses, pdrop=pdrop)
        update!(model, grads, optim)
        nw = sum(map(length,sentences))
        if (speed = inc(nwords, nw)) != nothing
            date("$(nwords.ncurr) words $(round(Int,speed)) wps $(losses[3]) avgloss")
            free_KnetArray()
        end
        if (niter+=1) >= maxiter; break; end
    end
    println()
end

function oracleloss(pmodel, sentences, vocab, arctype, feats; losses=nothing, pdrop=(0,0))
    parsers = map(arctype, sentences)
    mcosts = Array(Cost, parsers[1].nmove)
    parserdone = falses(length(parsers))
    totalloss = 0
    featmodel,mlpmodel = splitmodel(pmodel)

    nsteps = 0
    while !all(parserdone)
        fmatrix = features(parsers, feats, featmodel)
        #fmatrix = rand(FTYPE, 4664, length(parsers))
        if GPUFEATURES && gpu()>=0 #GPU
            @assert isa(getval(fmatrix),KnetArray{FTYPE,2})
        else #CPU
            @assert isa(getval(fmatrix),Array{FTYPE,2})
            if gpu()>=0
                fmatrix = KnetArray(fmatrix)
            end
        end
        scores = mlp(mlpmodel, fmatrix; pdrop=pdrop)
        logprob = logp(scores, 1)
        for (i,p) in enumerate(parsers)
            if parserdone[i]
                continue
            end
            movecosts(p, p.sentence.head, p.sentence.deprel, mcosts)
            goldmove = indmin(mcosts)
            if mcosts[goldmove] == typemax(Cost)
                parserdone[i] = true
                p.sentence.parse = p
            else
                totalloss -= logprob[goldmove,i]
                ##########
                # This part is done to make dynamic-oracle possible
                mmax = 0; smax = -Inf;
                logprob1 = Array((getval(logprob)));
                for m =1:size(logprob1, 1)
                    if mcosts[m] < typemax(Cost) && logprob1[m, i] > smax
                        smax = logprob1[m, i]; mmax = m;
                    end
                end
                move!(p, deepcopy(mmax))
                #########
                if losses != nothing
                    loss1 = -getval(logprob)[goldmove,i]
                    losses[1] += loss1
                    losses[2] += 1
                    if losses[2] < 1000
                        losses[3] = losses[1]/losses[2]
                    else
                        losses[3] = 0.999 * losses[3] + 0.001 * loss1
                    end
                end
            end
        end
        nsteps +=1;
    end
    return totalloss / length(sentences)
end

oraclegrad = grad(oracleloss)

using AutoGrad
import Base: sortperm, ind2sub
@zerograd sortperm(a;o...)
@zerograd ind2sub(a,i...)

# static oracle approx max ~.81 in Eng-Par-TUT
# Here is the general command to run the parser:
# julia main.jl --load data/english_chmodel.jld --datafiles /mnt/ai/data/nlp/conll17/ud-treebanks-v2.0/UD_English-LinES/en_lines-ud-train.conllu /mnt/ai/data/nlp/conll17/ud-treebanks
# For Hungarian:
# julia main.jl --load /mnt/ai/data/nlp/conll17/competition/chmodel_converted/hungarian_chmodel.jld --datafiles /mnt/ai/data/nlp/conll17/ud-treebanks-v2.0/UD_Hungarian/hu-ud-train.conllu /mnt/ai/data/nlp/conll17/ud-treebanks-v2.0/UD_Hungarian/hu-ud-dev.conllu


text:
# module Bukdu

function plugins
end

function before
end

function after
end


text:
# Run with `julia --project=<path/to/Project.toml> coverage.jl`
# e.g. `julia --project=. coverage.jl`
using Pkg
using Coverage

Pkg.activate(".")

# If Coverage is not installed correctly, it can fail but still exit without a
# proper error code, which causes the calling script to report success. So we
# wrap it in a try.
if !isdir(Pkg.dir("Coverage"))
    throw(LoadError(string("Coverage", 4, "Coverage is not installed correctly")))
end

# defaults to src/; alternatively, supply the folder name as argument
coverage = process_folder()
# Get total coverage for all Julia files
covered_lines, total_lines = get_summary(coverage)
# Or process a single file
# @show get_summary(process_file("src/MyPkg.jl"))
println("Covered lines | ", covered_lines)
println("Total lines   | ", total_lines)
println("Coverage      | ", round(covered_lines / total_lines * 100; digits = 1), "%")


text:
abstract type EventLikelihood{T<:Real} <: Likelihood{T} end

include("poisson.jl")

""" Return the labels in a vector of vectors for multiple outputs"""
function treat_labels!(y::AbstractArray{T,N},likelihood::L) where {T,N,L<:EventLikelihood}
    @assert T<:Integer "For event count target(s) should be integers (Bool,Integer or Float)"
    @assert N <= 2 "Target should be a matrix or a vector"
    if N == 1
        return [y],1,likelihood
    else
        return [y[:,i] for i in 1:size(y,2)],size(y,2),likelihood
    end
end


text:
"""
    smoptimize(f::Function, search_range::Array{Tuple{Float64,Float64},1}; options=Options())

Optimize the function `f` in the range `search_range` using a Radial Basis Function based surrogate model.
"""
function smoptimize(f::Function, search_range::Array{Tuple{Float64,Float64},1}; options::Options=Options())

    #Load some option values
    @unpack num_start_samples, sampling_plan_opt_gens,
            iterations, trace, create_final_surrogate = options

    reset_timer!()

    #Create sampling plan
    @timeit_debug "lhc" lhc_plan = scaled_LHC_sampling_plan(search_range,num_start_samples,sampling_plan_opt_gens;trace=trace)

    #Evaluate sampling plan
    @timeit_debug "f_eval" lhc_samples = f_opt_eval(f,lhc_plan;trace=trace)

    #Initialize variables to be returned
    sm_interpolant = nothing
    infill_type = Array{Symbol,1}(undef,0)
    infill_prediction = Array{Float64,1}(undef,0)
    optres = nothing
    infill_plan = Array{Float64,2}(undef,size(lhc_plan,1),0)
    infill_sample = Array{Float64,2}(undef,1,0)

    #Run the entire optimization iterations number of times
    for i = 1:iterations

        #Create the optimized Radial Basis Function interpolant
        samples_all = [lhc_samples infill_sample]
        plan_all = [lhc_plan infill_plan]
        @timeit_debug "sm_opt" begin
            sm_interpolant, optres = surrogate_model(plan_all, samples_all; options=options)
        end

        #Points to add to the sampling plan to improve the interpolant
        @timeit_debug "model_infill" begin
            infill_plan_new, infill_type_new, infill_prediction_new, options  = model_infill(search_range,plan_all,
                                                                        samples_all,sm_interpolant;options=options)
        end

        #Evaluate the new infill points
        print_iteration(trace,i,iterations)
        @timeit_debug "f_eval" infill_sample_new = f_opt_eval(f,infill_plan_new,samples_all;trace=trace)

        #Add infill points
        infill_plan = [infill_plan infill_plan_new]
        infill_sample = [infill_sample infill_sample_new]
        infill_type = [infill_type; infill_type_new]
        infill_prediction = [infill_prediction; infill_prediction_new]

    end

    if create_final_surrogate
        #Create the optimized Radial Basis Function interpolant
        samples_all = [lhc_samples infill_sample]
        plan_all = [lhc_plan infill_plan]
        sm_interpolant, optres = surrogate_model(plan_all, samples_all; options=options)
    end

    timeit_debug_enabled() && print_timer()

    return SurrogateResult( lhc_samples, lhc_plan, sm_interpolant,
                            optres, infill_sample, infill_type,
                            infill_plan, infill_prediction,options)
end


text:
using DashBootstrapComponents, DashHtmlComponents

row = html_div([
    dbc_row(
        [
            dbc_col(html_div("One of two columns"), width=4),
            dbc_col(html_div("One of two columns"), width=4),
        ],
        justify="start",
    ),
    dbc_row(
        [
            dbc_col(html_div("One of two columns"), width=4),
            dbc_col(html_div("One of two columns"), width=4),
        ],
        justify="center",
    ),
    dbc_row(
        [
            dbc_col(html_div("One of two columns"), width=4),
            dbc_col(html_div("One of two columns"), width=4),
        ],
        justify="end",
    ),
    dbc_row(
        [
            dbc_col(html_div("One of two columns"), width=4),
            dbc_col(html_div("One of two columns"), width=4),
        ],
        justify="between",
    ),
    dbc_row(
        [
            dbc_col(html_div("One of two columns"), width=4),
            dbc_col(html_div("One of two columns"), width=4),
        ],
        justify="around",
    ),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# A quadratic triangle, defined in 3D.

# Summary of intersection methods:
#       Triangulation
#           - Speed varies dramatically with precision
#               - The faster method when precision at or below 2 decimal places is desired
#           - Can provide false positives and false negatives due to approximation of the surface
#           - Constant time for a given number of triangles
#               - The only way to determine if two intersections exist is to test until two unique
#                 points are found. Since this happens infrequently, all points are tested, so time
#                 isn't wasted on control logic.
#       Newton-Raphson (iterative)
#           - Speed varies slightly based upon line segment length and orientation
#               - A shorter line segment will converge faster
#               - Converges based upon Jacobian matrix.
#                   - If derivatives are small, the iteration can become slow
#               - The faster method when precision beyond 2 decimal places is desired
#           - Less accurate
#               - May falsely give one intersection instead of two, especially for longer segments.
#               - This is due to the point of convergence being dependent on the initial guess
#                 point. The two starting points are placed close to the line segment start/stop
#                 to try to mitigate this.
#           - Precision to 6+ decimal places.
#       Overall
#           - Triangulation is predictable in speed, slow for high precision, and largely accurate
#           - Newton-Raphson is unpredictable in speed, fast for high precision, and
#               can be inaccurate for 2 intersections
#           - Consider the difference in timing between intersections for porting to GPU
#               - Newton-Raphson may cause thread divergence

struct Triangle6_3D{T <: AbstractFloat}
    # The points are assumed to be ordered as follows
    # p₁ = vertex A
    # p₂ = vertex B
    # p₃ = vertex C
    # p₄ = point on the quadratic segment from A to B
    # p₅ = point on the quadratic segment from B to C
    # p₆ = point on the quadratic segment from C to A
    points::NTuple{6, Point_3D{T}}
end

# Constructors
# -------------------------------------------------------------------------------------------------
Triangle6_3D(p₁::Point_3D{T},
             p₂::Point_3D{T},
             p₃::Point_3D{T},
             p₄::Point_3D{T},
             p₅::Point_3D{T},
             p₆::Point_3D{T}
            ) where {T <: AbstractFloat} = Triangle6_3D((p₁, p₂, p₃, p₄, p₅, p₆))

# Methods
# -------------------------------------------------------------------------------------------------
# Interpolation
function (tri6::Triangle6_3D{T})(r::R, s::S) where {T <: AbstractFloat,
                                                    R <: Real,
                                                    S <: Real}
    # See The Visualization Toolkit: An Object-Oriented Approach to 3D Graphics, 4th Edition
    # Chapter 8, Advanced Data Representation, in the interpolation functions section
    r_T = T(r)
    s_T = T(s)
    return (1 - r_T - s_T)*(2(1 - r_T - s_T) - 1)*tri6.points[1] +
                                     r_T*(2r_T-1)*tri6.points[2] +
                                     s_T*(2s_T-1)*tri6.points[3] +
                             4r_T*(1 - r_T - s_T)*tri6.points[4] +
                                       (4r_T*s_T)*tri6.points[5] +
                             4s_T*(1 - r_T - s_T)*tri6.points[6]
end

function (tri6::Triangle6_3D{T})(p::Point_2D{T}) where {T <: AbstractFloat,
                                                        R <: Real,
                                                        S <: Real}
    r_T = p[1]
    s_T = p[2]
    return (1 - r_T - s_T)*(2(1 - r_T - s_T) - 1)*tri6.points[1] +
                                     r_T*(2r_T-1)*tri6.points[2] +
                                     s_T*(2s_T-1)*tri6.points[3] +
                             4r_T*(1 - r_T - s_T)*tri6.points[4] +
                                       (4r_T*s_T)*tri6.points[5] +
                             4s_T*(1 - r_T - s_T)*tri6.points[6]
end

function derivatives(tri6::Triangle6_3D{T}, r::R, s::S) where {T <: AbstractFloat,
                                                               R <: Real,
                                                               S <: Real}
    # Let T(r,s) be the interpolation function for tri6
    # Returns ∂T/∂r, ∂T/∂s
    r_T = T(r)
    s_T = T(s)
    ∂T_∂r = (4r_T + 4s_T - 3)*tri6.points[1] +
                   (4r_T - 1)*tri6.points[2] +
            4(1 - 2r_T - s_T)*tri6.points[4] +
                       (4s_T)*tri6.points[5] +
                      (-4s_T)*tri6.points[6]

    ∂T_∂s = (4r_T + 4s_T - 3)*tri6.points[1] +
                   (4s_T - 1)*tri6.points[3] +
                      (-4r_T)*tri6.points[4] +
                       (4r_T)*tri6.points[5] +
            4(1 - r_T - 2s_T)*tri6.points[6]
    return ∂T_∂r, ∂T_∂s
end

function area(tri6::Triangle6_3D{T}; N::Int64=79) where {T <: AbstractFloat}
    # Numerical integration required. Gauss-Legendre quadrature over a triangle is used.
    # Let T(r,s) be the interpolation function for tri6,
    #                             1 1-r                          N
    # A = ∬ ||∂T/∂r × ∂T/∂s||dA = ∫  ∫ ||∂T/∂r × ∂T/∂s|| ds dr = ∑ wᵢ||∂T/∂r(rᵢ,sᵢ) × ∂T/∂s(rᵢ,sᵢ)||
    #      D                      0  0                          i=1
    #
    # N is the number of points used in the quadrature.
    # See tuning/Triangle6_3D_area.jl for more info on how N = 79 was chosen.
    w, r, s = gauss_legendre_quadrature(tri6, N)
    a = T(0)
    for i in 1:N
        ∂T_∂r, ∂T_∂s = derivatives(tri6, r[i], s[i])
        a += w[i] * norm(∂T_∂r × ∂T_∂s)
    end
    return a
end

function triangulate(tri6::Triangle6_3D{T}, N::Int64) where {T <: AbstractFloat}
    # N is the number of divisions of each edge
    triangles = Vector{Triangle_3D{T}}(undef, (N+1)*(N+1))
    if N === 0
        triangles[1] = Triangle_3D(tri6.points[1], tri6.points[2], tri6.points[3])
    else
        i = 1
        for S = 1:N, R = 0:N-S
            triangles[i]   = Triangle_3D(tri6(    R/(N+1),     S/(N+1)),
                                         tri6((R+1)/(N+1),     S/(N+1)),
                                         tri6(    R/(N+1), (S+1)/(N+1)))
            triangles[i+1] = Triangle_3D(tri6(    R/(N+1),     S/(N+1)),
                                         tri6((R+1)/(N+1), (S-1)/(N+1)),
                                         tri6((R+1)/(N+1),     S/(N+1)))
            i += 2
        end
        j = (N+1)*N + 1
        for S = 0:0, R = 0:N-S
            triangles[j] = Triangle_3D(tri6(    R/(N+1),     S/(N+1)),
                                       tri6((R+1)/(N+1),     S/(N+1)),
                                       tri6(    R/(N+1), (S+1)/(N+1)))
            j += 1
        end
    end
    return triangles
end

# Triangulate, then intersect
function intersect(l::LineSegment_3D{T}, tri6::Triangle6_3D{T};
        N::Int64 = 25) where {T <: AbstractFloat}
    triangles = triangulate(tri6, N)
    npoints = 0
    p₁ = Point_3D(T, 0)
    p₂ = Point_3D(T, 0)
    intersections = l .∩ triangles
    bools = map(x->x[1], intersections)
    points = map(x->x[2], intersections)
    npoints = count(bools)
    p₁ = Point_3D(T, 0)
    p₂ = Point_3D(T, 0)
    if npoints === 0
        return false, 0, p₁, p₂
    elseif npoints === 1
        p₁ = points[argmax(bools)]
        return true, 1, p₁, p₂
    elseif npoints === 2
        indices = findall(bools)
        p₁ = points[indices[1]]
        p₂ = points[indices[2]]
        # Check uniqueness
        if p₁ ≈ p₂
            return true, 1, p₁, p₂
        else
            return true, 2, p₁, p₂
        end
    else
        # Account for 3 points and 4 points?
        # If intersection is on edge shared by two triangles on entrance and/or exit 3/4 intersections
        # can be detected
        return true, -1, p₁, p₂
    end
end

function real_to_parametric(p::Point_3D{T}, tri6::Triangle6_3D{T}; N::Int64=30) where {T <: AbstractFloat}
    # Convert from real coordinates to the triangle's local parametric coordinates using the
    # the Newton-Raphson method. N is the max number of iterations
    # If a conversion doesn't exist, the minimizer is returned.
    r = T(1//3) # Initial guess at triangle centroid
    s = T(1//3)
    err₁ = p - tri6(r, s)
    for i = 1:N
        ∂T_∂r, ∂T_∂s = derivatives(tri6, r, s)
        J = hcat(∂T_∂r.x, ∂T_∂s.x)
        Δr, Δs = J \ err₁.x
        r = r + Δr
        s = s + Δs
        err₂ = p - tri6(r, s)
        if norm(err₂ - err₁) < 1e-6
            break
        end
        err₁ = err₂
    end
    return Point_2D(r, s)
end

# A more exact intersection algorithm that triangulation, uses Newton-Raphson.
function intersect_iterative(l::LineSegment_3D{T}, tri6::Triangle6_3D{T};
        N::Int64=30) where {T <: AbstractFloat}
    p₁ = Point_3D(T, 0)
    p₂ = Point_3D(T, 0)
    npoints = 0
    u⃗ = l.points[2] - l.points[1]
    ray_start = real_to_parametric(l(0), tri6; N=10) # closest r,s to the ray start
    ray_stop  = real_to_parametric(l(1), tri6; N=10) # closest r,s to the ray stop
    # The parametric coordinates corresponding to the start of the line segment
    r₁ = ray_start[1]
    s₁ = ray_start[2]
    t₁ = T(0)
    # The parametric coordinates corresponding to the start of the line segment
    r₂ = ray_stop[1]
    s₂ = ray_stop[2]
    t₂ = T(1)
    # Iteration for point 1
    err₁ = tri6(r₁, s₁) - l(t₁)
    for i = 1:N
        ∂r₁, ∂s₁ = derivatives(tri6, r₁, s₁)
        J₁ = hcat(∂r₁.x, ∂s₁.x, -u⃗.x)
        # If matrix is singular, it's probably bad luck. Just perturb it a bit.
        if abs(det(J₁)) < 1e-5
            r₁, s₁, t₁ = [r₁, s₁, t₁] + rand(3)/10
        else
            r₁, s₁, t₁ = [r₁, s₁, t₁] - inv(J₁) * err₁.x
            errₙ = tri6(r₁, s₁) - l(t₁)
            if norm(errₙ - err₁) < 5e-6
                break
            end
            err₁ = errₙ
        end
    end
    # Iteration for point 2
    err₂ = tri6(r₂, s₂) - l(t₂)
    for j = 1:N
        ∂r₂, ∂s₂ = derivatives(tri6, r₂, s₂)
        J₂ = hcat(∂r₂.x, ∂s₂.x, -u⃗.x)
        if abs(det(J₂)) < 1e-5
            r₂, s₂, t₂ = [r₂, s₂, t₂] + rand(3)/10
        else
            r₂, s₂, t₂ = [r₂, s₂, t₂] - inv(J₂) * err₂.x
            errₘ = tri6(r₂, s₂) - l(t₂)
            if norm(errₘ - err₂) < 5e-6
                break
            end
            err₂ = errₘ
        end
    end

    p₁ = l(t₁)
    if (0 ≤ r₁ ≤ 1) && (0 ≤ s₁ ≤ 1) && (0 ≤ t₁ ≤ 1) && (p₁ ≈ tri6(r₁, s₁))
        npoints += 1
    end

    p₂ = l(t₂)
    if (0 ≤ r₂ ≤ 1) && (0 ≤ s₂ ≤ 1) && (0 ≤ t₂ ≤ 1) && (p₂ ≈ tri6(r₂, s₂))
        npoints += 1
        # If only point 2 is valid, return it as p₁
        # If points are duplicate, reduce npoints
        if npoints === 2 && p₁ ≈ p₂
            npoints -= 1
        elseif npoints === 1
            p₁ = p₂
        end
    end
    return npoints > 0, npoints, p₁, p₂
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
using NLopt, Test

# Objective function
function rosenbrock(x::Vector, grad::Vector)
    if length(grad) > 0
        grad[1] = -400 * x[1] * (x[2] - x[1]^2) - 2 * (1 - x[1])
        grad[2] = 200 * (x[2] - x[1]^2)
    end
    (1 - x[1])^2 + 100*(x[2] - x[1]^2)^2
end

function ineq01(x::Vector, grad::Vector)
    if length(grad) > 0
        grad[1] = 1
        grad[2] = 2
    end
    x[1] + 2 * x[2] - 1
end

function ineq02(x::Vector, grad::Vector)
    if length(grad) > 0
        grad[1] = 2*x[1]
        grad[2] = 1
    end
    x[1]^2 + x[2] - 1
end

function ineq03(x::Vector, grad::Vector)
    if length(grad) > 0
        grad[1] = 2*x[1]
        grad[2] = -1
    end
    x[1]^2 - x[2] - 1
end

function eq01(x::Vector, grad::Vector)
    if length(grad) > 0
        grad[1] = 2
        grad[2] = 1
    end
    2 * x[1] + x[2] - 1
end

opt = Opt(:LD_SLSQP, 2)
opt.lower_bounds = [0, -0.5]
opt.upper_bounds = [1, 2]
opt.xtol_rel = 1e-21
opt.min_objective = rosenbrock
opt.inequality_constraint = ineq01
opt.inequality_constraint = ineq02
opt.inequality_constraint = ineq03
opt.equality_constraint = eq01

(minf,minx,ret) = optimize(opt, [0.5, 0])
println("got $minf at $minx with constraints (returned $ret)")
@test minx[1] ≈ 0.4149 rtol=1e-3
@test minx[2] ≈ 0.1701 rtol=1e-3

remove_constraints!(opt)
(minf,minx,ret) = optimize(opt, [0.5, 0])
println("got $minf at $minx after removing constraints (returned $ret)")
@test minx[1] ≈ 1 rtol=1e-5
@test minx[2] ≈ 1 rtol=1e-5


text:
mutable struct Edge
    to::Int
    w::Float64
end

function main()
    N, M = parse.(Int, split(readline()))
    G = [Vector{Edge}() for _ = 1:N]
    for _ = 1:M
        a, b, c = parse.(Int, split(readline()))
        push!(G[a], Edge(b, -c))
    end

    # Bellman-Ford
    has_negative_cycle = false
    dists = fill(Inf, N)
    dists[1] = 0.0
    for i = 1:N
        for v = 1:N
            if !isinf(dists[v])
                for e ∈ G[v]
                    if dists[e.to] > dists[v] + e.w
                        dists[e.to] = dists[v] + e.w
                        if e.to == N && i == N
                            has_negative_cycle = true
                        end
                    end
                end
            end
        end

    end

    if has_negative_cycle
        println("inf")
    else
        println("$(-floor(Int, dists[end]))")
    end
end

main()


text:
using LightGraphs

const input = map(l -> split(l, ""), readlines(length(ARGS) == 0 ? "example1.txt" : ARGS[1]))
const m = length(input) - 4
const n = length(input[1]) - 4

function vertex_id(i, j)
    #return input[i][j] == "." ? (i + j * im) : 0
    return input[i][j] == "." ? (1 + (i - 3)n + j - 3) : 0
end

# Assign numeric names to positions in the maze.
const v = [vertex_id(i + 3,j + 3) for i=0:m-1, j=0:n-1]
if length(v) < 1000
    display(v)
    println()
end

# Construct a simple graph of adjacent positions.
# The graph will be initialized with a much larger capacity for vertices than
# it needs. Hopefully the implementation uses a sparse array.
G = Graph(length(v))

for row in 1:(size(v)[1]-1)
    for col in 1:(size(v)[2]-1)
        current = v[row, col]
        if current > 0
            right = v[row + 1, col]
            down = v[row, col + 1]
            if right > 0
                add_edge!(G, current, right) || error("Could not insert edge")
            end
            if down > 0
                add_edge!(G, current, down) || error("Could not insert edge")
            end
        end
    end
end

# Identify the portals and join them.
const portals = Dict()
for row in 3:length(input)-2
    for col in 3:length(input[1])-2
        if input[row][col] == "."
            label = ""
            # letters above
            if isletter(input[row-1][col][1]) && isletter(input[row-2][col][1])
                label = input[row-2][col][1] * input[row-1][col][1]
            # letters below
            elseif isletter(input[row+1][col][1]) && isletter(input[row+2][col][1])
                label = input[row+1][col][1] * input[row+2][col][1]
            # letters left
            elseif isletter(input[row][col-1][1]) && isletter(input[row][col-2][1])
                label = input[row][col-2][1] * input[row][col-1][1]
            # letters right
            elseif isletter(input[row][col+1][1]) && isletter(input[row][col+2][1])
                label = input[row][col+1][1] * input[row][col+2][1]
            end
            if label != ""
                #print("Found portal $(label) at M[$(row),$(col)]")
                #println(" (id=$(vertex_id(row,col)))")
                if !(haskey(portals, label))
                    portals[label] = Array{Int,1}(undef, 0)
                end
                push!(portals[label], vertex_id(row, col))
            end
        end
    end
end

for kv in portals
    length(last(kv)) == 2 && (add_edge!(G, last(kv)[1], last(kv)[2]) ||
        error("Unable to insert edge"))
end

#println(collect(edges(G)))

spf = dijkstra_shortest_paths(G, first(portals["AA"]))
println("Day 20 Part 1: $(spf.dists[first(portals["ZZ"])])")


text:
function mfb_rgb(r::UInt8, g::UInt8, b::UInt8)::UInt32
    reinterpret(UInt32, [b, g, r, zero(UInt8)])[1]
end

"""
    mfb_rgb(r, g, b)::UInt32

convert 8 bit color channels to a single 32 bit buffer.
Input should be of some integer type in the range 0-255.
`InexactError` will be thrown if the inputs are fractional, or larger than 255.
"""
mfb_rgb(r, g, b) = mfb_rgb(UInt8(r), UInt8(g), UInt8(b))

"""
    mfb_rgb(c::Colorant)::UInt32

convert a color object from Colors.jl into a 32 bit buffer.
"""
function mfb_rgb(c::Colorant)::UInt32
    cRGB24 = convert(RGB24, c)
    return reinterpret(UInt32, cRGB24)
end

foreach(names(@__MODULE__, all=true)) do s
   if startswith(string(s), "mfb_")
       @eval export $s
   end
end

function mfb_set_active_callback(window, active::Function)
    active_c = @cfunction($active, Cvoid, (Ptr{Cvoid}, Bool))
    mfb_set_active_callback(window, active_c)
end

function mfb_set_resize_callback(window, resize::Function)
    resize_c = @cfunction($resize, Cvoid, (Ptr{Cvoid}, Cint, Cint))
    mfb_set_resize_callback(window, resize_c)
end

function mfb_set_keyboard_callback(window, keyboard::Function)
    keyboard_c = @cfunction($keyboard, Cvoid, (Ptr{Cvoid}, mfb_key, mfb_key_mod, Bool))
    mfb_set_keyboard_callback(window, keyboard_c)
end

function mfb_set_char_input_callback(window, char_input::Function)
    char_input_c = @cfunction($char_input, Cvoid, (Ptr{Cvoid}, Cuint, Cuint))
    mfb_set_char_input_callback(window, char_input_c)
end

function mfb_set_mouse_button_callback(window, mouse_btn::Function)
    mouse_btn_c = @cfunction($mouse_btn, Cvoid, (Ptr{Cvoid}, mfb_mouse_button, mfb_key_mod, Bool))
    mfb_set_mouse_button_callback(window, mouse_btn_c)
end

function mfb_set_mouse_move_callback(window, mouse_move::Function)
    mouse_move_c = @cfunction($mouse_move, Cvoid, (Ptr{Cvoid}, Cint, Cint))
    mfb_set_mouse_move_callback(window, mouse_move_c)
end

function mfb_set_mouse_scroll_callback(window, mouse_scroll::Function)
    mouse_scroll_c = @cfunction($mouse_scroll, Cvoid, (Ptr{Cvoid}, mfb_key_mod, Cfloat, Cfloat))
    mfb_set_mouse_scroll_callback(window, mouse_scroll_c)
end


text:
using DataStructures
using Test

# Function to check to make sure we are on the canvas
function inbounds(canvas_size, loc)

    # Make sure we are not beneath or to the left of the canvas
    if minimum(Tuple(loc)) < 1
        return false

    # Make sure we are not to the right of the canvas
    elseif loc[2] > canvas_size[2]
        return false

    # Make sure we are not above the canvas
    elseif loc[1] > canvas_size[1]
        return false
    else
        return true
    end
end

function find_neighbors(canvas, loc::CartesianIndex, old_val, new_val)

    # Finding north, south, east, west neighbors
    possible_neighbors = [loc + CartesianIndex(0, 1),
                          loc + CartesianIndex(1, 0),
                          loc + CartesianIndex(0, -1),
                          loc + CartesianIndex(-1, 0)]

    # Exclusing neighbors that should not be colored
    neighbors =  []
    for possible_neighbor in possible_neighbors
        if inbounds(size(canvas), possible_neighbor) &&
           canvas[possible_neighbor] == old_val
            push!(neighbors, possible_neighbor)
        end
    end

    return neighbors
end

function stack_fill!(canvas, loc::CartesianIndex, old_val, new_val)
    if new_val == old_val
        return
    end

    s = Stack{CartesianIndex}()
    push!(s, loc)

    while length(s) > 0
        current_loc = pop!(s)
        if canvas[current_loc] == old_val
            canvas[current_loc] = new_val
            possible_neighbors = find_neighbors(canvas, current_loc,
                                                old_val, new_val)
            for neighbor in possible_neighbors
                push!(s,neighbor)
            end
        end

    end
end

function queue_fill!(canvas, loc::CartesianIndex, old_val, new_val)
    if new_val == old_val
        return
    end

    q = Queue{CartesianIndex}()
    enqueue!(q, loc)

    # Coloring the initial location
    canvas[loc] = new_val

    while length(q) > 0
        current_loc = dequeue!(q)

        possible_neighbors = find_neighbors(canvas, current_loc,
                                            old_val, new_val)

        # Coloring as we are enqueuing neighbors
        for neighbor in possible_neighbors
            canvas[neighbor] = new_val
            enqueue!(q,neighbor)
        end

    end
end

function recursive_fill!(canvas, loc::CartesianIndex, old_val, new_val)

    if (old_val == new_val)
        return
    end

    canvas[loc] = new_val

    possible_neighbors = find_neighbors(canvas, loc, old_val, new_val)
    for possible_neighbor in possible_neighbors
        recursive_fill!(canvas, possible_neighbor, old_val, new_val)
    end
end

function main()

    # Creation of a 5x5 grid with a single row of 1.0 elements
    grid = zeros(5,5)
    grid[3,:] .= 1

    # Create solution grid
    answer_grid = zeros(5,5)
    answer_grid[1:3, :] .= 1

    # Start filling at 1,1
    start_loc = CartesianIndex(1,1)

    @testset "Fill Methods" begin
        # Use recursive method and reinitialize grid
        recursive_fill!(grid, start_loc, 0.0, 1.0)
        @test grid == answer_grid

        grid[1:2,:] .= 0

        # Use queue method and reinitialize grid
        queue_fill!(grid, start_loc, 0.0, 1.0)
        @test grid == answer_grid

        grid[1:2,:] .= 0

        # Use stack method and reinitialize grid
        stack_fill!(grid, start_loc, 0.0, 1.0)
        @test grid == answer_grid
    end

end

main()