text:
using Test
using TightlyBound
using Suppressor

#include("../includes_laguerre.jl")
#include("../Ewald.jl")

TESTDIR=TightlyBound.TESTDIR

function loaddata(dirs; scf=true)
    tbc_list  = []
    dft_list = []

    for t in dirs
        #                println(t*"/qe.save")
        tfull = "$TESTDIR/"*t
        dft = TightlyBound.QE.loadXML(tfull*"/qe.save")
        tbc = []
        tbc_scf = []
        try
            if scf
                tbc_scf = TightlyBound.TB.read_tb_crys("projham_scf.xml.gz", directory=tfull)
            else
                tbc_scf = TightlyBound.TB.read_tb_crys("projham.xml.gz", directory=tfull)
            end
        catch
            tbc = TightlyBound.AtomicProj.projwfx_workf(dft, directory=tfull, writefile="projham.xml", skip_og=true, skip_proj=true, freeze=true, localized_factor = 0.15)
            if scf
                tbc_scf = TightlyBound.SCF.remove_scf_from_tbc(tbc)
                TightlyBound.TB.write_tb_crys(t*"/projham_scf.xml.gz", tbc_scf)
            else
                tbc_scf = tbc
            end

        end

        push!(dft_list, dft)
        push!(tbc_list, tbc_scf)
    end
    return     tbc_list, dft_list
end

function test_force()

    @testset "testing force dimer" begin

        if true
#        @suppress begin
            ft = open("$TESTDIR/data_forces/fil_MgS_dimer", "r");
            dirst = readlines(ft);
            close(ft);

            #        println(dirst)

#            for scf = [false true]
            f_cart = zeros(2,3)
            f_cartFD = 0.0

            for scf = [false, true]
                @suppress begin
                tbc_list, dft_list = loaddata(dirst, scf=scf);
                     database_rec = TightlyBound.FitTB.do_fitting_recursive(tbc_list,dft_list = dft_list,  fit_threebody=false, fit_threebody_onsite=false);

                x = 4;
                smearing = 0.01;

                    #en, tbc_x, flag = scf_energy(tbc_list[x].crys, database = database_rec)
                    #en, f_cart,stress = TightlyBound.Force_Stress.get_energy_force_stress(tbc_x, database_rec,   smearing = smearing);

                    en, f_cart,stress = TightlyBound.Force_Stress.get_energy_force_stress(tbc_list[x].crys, database_rec,   smearing = smearing);

                enFD, f_cartFD = TightlyBound.Force_Stress.finite_diff(tbc_list[x].crys, database_rec,1, 3,   smearing = smearing);
                end

#                println(scf, " xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx  ", f_cartFD, "  ", f_cart[1,3])

                @test abs(f_cartFD - f_cart[1,3]) < 1e-3
#                @test abs(f_cartFD - f_cart_fft[1,3]) < 1e-3
                #            println("SCF $scf TEST1 finite diff: ", f_cartFD , " autodiff:   ", f_cart[1,3], " xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
                #            println("TEST1 dft ref ", dft_list[x].forces[1,3])

            end
            #        x = 3;
            #        smearing = 0.01;
            #        en, f_cart,stress = Force_Stress.get_energy_force_stress(tbc_list[x].crys, database_rec,   smearing = smearing);
            #        enFD, f_cartFD = Force_Stress.finite_diff(tbc_list[x].crys, database_rec,1, 3,   smearing = smearing);

            #        @test abs(f_cartFD - f_cart[1,3]) < 1e-3
            #        println("TEST3 finite diff: ", f_cartFD , " autodiff:   ", f_cart[1,3], " xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
            #        println("TEST3 dft ref ", dft_list[x].forces[1,3])

        end
    end
end

function test_stress()

    @testset "testing force znse" begin
        @suppress begin

            ft = open("$TESTDIR/data_forces/fil_MgS_znse", "r");
            dirst = readlines(ft);
            close(ft);

            #        println(dirst)

            tbc_list, dft_list = loaddata(dirst, scf=false);

            database = TightlyBound.FitTB.do_fitting(tbc_list, fit_threebody=false, fit_threebody_onsite=false, do_plot=false);
            #        database = FitTB.do_fitting_recursive(tbc_list,dft_list = dft_list,  fit_threebody=true, fit_threebody_onsite=false);

            x = 1;
            smearing = 0.01;
            en, f_cart, stress = TightlyBound.Force_Stress.get_energy_force_stress(tbc_list[x].crys, database,   smearing = smearing);

            enFD, f_cartFD = TightlyBound.Force_Stress.finite_diff(tbc_list[x].crys, database,1, 3,   smearing = smearing);

            #        println("TEST force finite diff: ", f_cartFD , " autodiff:   ", f_cart[1,3], " xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
            #        println("TEST dft ref ", dft_list[x].forces[1,3])
            @test abs(f_cartFD - f_cart[1,3]) < 1e-4

            x=1
            enFD, stressFD = TightlyBound.Force_Stress.finite_diff(tbc_list[x].crys, database,1, 1, stress_mode=true,  smearing = smearing);

            #        println("TEST stress11 finite diff: ", stressFD , " autodiff:   ", stress[1,1], " xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
            #        println("TEST dft ref ", dft_list[x].stress[1,1])
            @test abs(stressFD - stress[1,1]) < 1e-5

            x=1
            enFD, stressFD = TightlyBound.Force_Stress.finite_diff(tbc_list[x].crys, database,1, 2, stress_mode=true,  smearing = smearing);

            #        println("TEST stress12 finite diff: ", stressFD , " autodiff:   ", stress[1,2], " xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx")
            #        println("TEST dft ref ", dft_list[x].stress[1,2])
            @test abs(stressFD - stress[1,2]) < 1e-5

            #        println("done")
        end
    end
end

test_force()

#println("sleep ")
#sleep(3)

test_stress()


text:
function minmod(a,b,c)
  if (a > 0 && b > 0 && c > 0)
    min(a,b,c)
  elseif (a < 0 && b < 0 && c < 0)
    max(a,b,c)
  else
    zero(a)
  end
end

function minmod(a,b)
  0.5*(sign(a)+sign(b))*min(abs(a),abs(b))
end

"Apply slope limiter πᴺ to u assuming u an Nth polynomial"
mutable struct DGLimiter{pType}
  problem::AbstractConservationLawProblem
  basis::PolynomialBasis
  limiter::AbstractDGLimiter
  params::pType
end

function DGLimiter(problem, basis, limiter)
    DGLimiter(problem, basis, limiter, params_init(limiter, basis, problem))
end

# Trivial parameter initialization
params_init(limiter::AbstractDGLimiter, basis::PolynomialBasis, problem::AbstractConservationLawProblem) = nothing

function (dgLimiter::DGLimiter)(u, f, t)
    @unpack problem, basis, limiter, params = dgLimiter
    apply_limiter!(u,f,t,problem, basis, params, limiter)
end

# Auxiliary functions
"""
get_cell_averages!(uavg, basis, uh, NN)
Compute cell averages of polynomial system u,
using modal values `uh`, a polynomial `basis` and number of
variables `NC`
"""
function get_cell_averages!(uavg::AbstractArray{T,2}, basis::PolynomialBasis, uh::AbstractArray{T,2}, NC::Int) where {T}
    NN = basis.order + 1
    for k = 1:NC
        uavg[k,:] = [basis.φ[1,1]*uh[(k-1)*NN+1,i] for i in 1:size(uh,2)]
    end
end

#MUSCL Limiter Van Leer
#Reference: Hesthaven, Warburton, Nodal Discontinuous Galerkin Methods Algorithms, Analysis and
#applications, pp 150.
immutable Linear_MUSCL_Limiter <: AbstractDGLimiter end

"""
apply_limiter!(u,f,t,mesh, basis, limiter::Linear_MUSCL_Limiter)
    apply Linear MUSCL limiter on function u
  """
function apply_limiter!(u,f,t,prob, basis, params, imiter::Linear_MUSCL_Limiter)
    # Compute modal coefficients
    mesh = prob.mesh; NC = prob.numvars
    NN = basis.order + 1
    uh = myblock(basis.invφ,NC)*u
    @assert basis.order > 1 "MUSCL limiter requires at least linear polynomials"
    #Extract linear polynomial
    ul = uh[:,:];
    for k = 1:NC;ul[((k-1)*NN+3):k*NN,:] = 0.0;end

    # Compute cell averages
    uavg = zeros(eltype(u), NC, numcells(mesh))
    get_cell_averages!(uavg, basis, uh, NC)

    #compute nodes in real coordinates
    x1 = zeros(basis.order+1,numcells(mesh))
    for k in 1:numcells(mesh)
        x1[:,k] = reference_to_interval(basis.nodes, (cell_faces(mesh)[k],cell_faces(mesh)[k+1]))
    end
    x0 = cell_centers(mesh)
    # Compute derivative
    h = cell_volumes(mesh)
    ux = 2*(1./h)'.*(myblock(basis.dφ,NC)*ul)

    ve = apply_boundary(uavg, mesh)
    uavgp1 = ve[:,3:end];  uavgm1 = ve[:,1:end-2]

    #limit nodal cell values
    for i in 1:size(u,2)
        for k = 1:NC
            for j in 1:NN
                u[(k-1)*NN+j,i] = uavg[k,i]+(x1[j,i]-x0[i])*
                minmod(ux[(k-1)*NN+j,i],(uavgp1[k,i]-uavg[k,i])/h[i],(uavg[k,i]-uavgm1[k,i])/h[i])
            end
      end
    end
end

#WENO limiter by Zhong-Shu (2013)
mutable struct WENO_Limiter <: AbstractDGLimiter end

function params_init(WENO_Limiter::WENO_Limiter, basis::PolynomialBasis, problem::AbstractConservationLawProblem)
    return WENOLimiterWeights(basis.order, basis.invφ)
end

"""
function WENOLimiterWeights(m, iV);
Compute operators to enable evaluation of WENO smoothness
indicator and WENO polynomial of order m.
"""
function WENOLimiterWeights(m,iV);
    Q = zeros(m+1,m+1)
    Pmat = zeros(m+1,m+1)
    Xm = Pmat; Xp = Pmat;

    # Compute quadrature points
    (x,w) = gausslegendre(m+1)
    Λ = diagm(w);

    # Initial matrices of Legendre polynomails
    for n = 0:m
      p = poly_legendre(n)
      # Eval nodes
      Pmat[n+1,:] = polyval(p, x)
      Xm[n+1,:] = polyval(p, x-2)
      Xp[n+1,:] = polyval(p, x+2)
    end

    # Compute matrices corresponding to increasing order of derivative
    for l=1:m
        # Set up operator to recover derivaties
        A = zeros(m+2-l,m+2-l); A[1,1] = 1/sqrt((2*l+1)*(2*l-1))
        A[m+2-l,m+2-l] = 1/(sqrt(2*(m+2)+1)*sqrt(2*(m+2)-1))
        for i=2:m-l+1
            Ah = 1/(sqrt(2*(l-1+i)+1)*sqrt(2*(l-1+i)-1))
            A[i,i] = Ah; A[i+1,i-1] = -Ah
        end

        # Recover derivatives at quadrature points
        Ph1 = A\Pmat[l:m+1,:]
        Pmat[1:l,:]=0; Pmat[l+1:m+1,:] = Ph1[1:m-l+1,:]

        # Compute smoothness operator for order l and update
        Qh = Pmat*Λ*Pmat'
        Q = Q + 2^(2*l-1)*Qh
    end

    # Initialize operator for smoothness indicator in nodal space
    Q = basis.invφ'*Q*basis.invφ;

    # Initialize interpolation matrices
    Xp = basis.invφ'*Xp; Xm = basis.invφ'*Xm
    return Q,Xm,Xp
end

"""
apply_limiter!(u,f,t,mesh, basis, limiter::WENO_Limiter)
    apply Zhong Shu WENO limiter on function u
  """
function apply_limiter!(u,f,t,prob, basis, params, limiter::WENO_Limiter)
    mesh = prob.mesh; NC = prob.numvars; NN = basis.order + 1
    (Q, Xp, Xm) = params
    Xp = myblock(Xp,NC)
    Xm = myblock(Xm,NC)
    iV = myblock(basis.invφ,NC)
    V = myblock(basis.φ,NC)
    eps0=1.0e-6;

    # Set constants for limiting
    eps1 = 1e-10; p=1;
    γm1 =0.001; γ0 = 0.998; γp1 = 0.001;

    # Compute cell averages and cell centers
    uh = iV*u
    uavg = zeros(eltype(u), NC, numcells(mesh))
    get_cell_averages!(uavg, basis, uh, NC)

    # Compute extended polynomials with zero cell averages
    ue = apply_boundary(u, mesh);
    Pm = Xp'*ue; Pp = Xm'*ue;
    Ph = iV*Pm; Ph[1:NN:end,:]=0; Pm = V*Ph;
    Ph = iV*Pp; Ph[1:NN:end,:]=0; Pp = V*Ph;

    # Extend cell averages
    ve = apply_boundary(uavg, mesh)

    # extract end values and cell averages for each element
    uel = u[1:NN:end,:]; uer = u[NN:NN:end,:]
    vj = uavg; vjm = ve[:,1:end-2]; vjp = ve[:,3:end]

    for k = 1:NC
        # Find elements that require limiting
        vel = vj[k,:] - minmod.(vj[k,:]-uel[k,:],vj[k,:]-vjm[k,:],vjp[k,:]-vj[k,:])
        ver = vj[k,:] + minmod.(uer[k,:]-vj[k,:],vj[k,:]-vjm[k,:],vjp[k,:]-vj[k,:])
        ids = union(find(x->abs(x)>eps0, vel-uel[k,:]), find(x->abs(x)>eps0, ver-uer[k,:]))

        # Apply limiting when needed
        if (!isempty(ids))
            # Extract local polynomials
            pm1 = Pm[(k-1)*NN+1:k*NN,ids] + ones(NN,1)*vj[k,ids]'
            p0  =  u[(k-1)*NN+1:k*NN,ids]
            pp1 = Pp[(k-1)*NN+1:k*NN,ids+2] + ones(NN,1)*vj[k,ids]'

            # Compute smoothness indicators and WENO weights
            βm1 = diag(pm1'*Q*pm1); αm1 = γm1./(eps1 + βm1).^(2*p);
            β0 = diag(p0'*Q*p0); α0  =  γ0./(eps1 + β0).^(2*p);
            βp1 = diag(pp1'*Q*pp1); αp1 = γp1./(eps1 + βp1).^(2*p);

            αs = αm1 + α0 + αp1;
            omm1 = αm1./αs; om0  = α0./αs; omp1 = αp1./αs;

            # Compute limited function
            u[(k-1)*NN+1:k*NN,ids] = pm1*diagm(omm1) + p0*diagm(om0) + pp1*diagm(omp1)
        end
    end
end

#TODO: Not Working
# immutable MUSCLV2Limiter <: AbstractDGLimiter end
#
# function apply_limiter(u,f,t,mesh, basis, limiter::MUSCLV2Limiter)
#     # Compute cell averages
#     ̄ū = [0.5*dot(basis.φ*u[:,i], basis.weights) for i in 1:size(u,2)]
#
#     # find end values of each element
#     ue1 = u[1,:];ue2=u[end,:]
#
#     #find cell averages
#     vk = ū; vkm1=[ū[1],ū[1:end-1]...];vkp1=[ū[2:end]...,ū[end]]
#
#     #apply reconstruction to find elements in need of limiting
#     ve1 = vk - minmod.(vk-ue1, vk -vkm1,vkp1-vk)
#     ve2 = vk + minmod.(ue2-vk,vk-vkm1,vkp1-vk)
#     tol = 1e-8
#     idx = (abs.(ve1-ue1).>tol)|(abs.(ve2-ue2).>tol)
#     if (!isempty(idx))
#       h2 = 2./diff(cell_faces(mesh))
#       x1 = zeros(basis.order+1,numcells(mesh))
#       for k in 1:numcells(mesh)
#         x1[:,k] = reference_to_interval(basis.nodes, (cell_faces(mesh)[k],cell_faces(mesh)[k+1]))
#       end
#       x0 = cell_centers(mesh)
#       ux = 0.5*h2.*((basis.dφ*u)'*basis.weights)
#       #println("idx: ", idx)
#       for j in 1:size(u,1)
#         u[j,idx] = vk[idx]+(x1[j,idx]-x0[idx]).*
#         minmod.(ux[idx],h2[idx].*(vkp1[idx]-vk[idx]),h2[idx].*(vk[idx]-vkm1[idx]))
#       end
#     end
# end

ext: jl

max_stars_repo_name: Paulms/DiscontinuousGalerkin

max_issues_repo_name: Paulms/DiscontinuousGalerkin

max_forks_repo_name: Paulms/DiscontinuousGalerkin


text:
#=
A file to hold the helper functions used for working with cubic polynomials
=#

"""A helper function; Returns true if any of the arguments is NaN"""
function any_nan{T}(x::T...)
	for n in x
		if isnan(n)
			return true
		end
	end

	#if the loop finished, then they are all not NAN so return false
	return false
end

function maxcubic{T}(x0::T, x1::T, a::T, b::T, c::T, d::T) #a+bx+cx^2+dx^3
	if any_nan(x0, x1, a, b, c, d)
		throw(ArgumentError("Maximum of a cubic with NaN parameters is not defined"))
	end

  if x0 == x1
    return (x0, @evalpoly(x0, a, b, c, d))
  end
  if d == 0 #then the max is at an endpoint. #FIXME: I don't think so
    endpt1_val = @evalpoly(x0, a, b, c, d)
    endpt2_val = @evalpoly(x1, a, b, c, d)
    return endpt2_val > endpt1_val ? (endpt2_val, x1) : (endpt1_val, x0)
  end

  #Otherwise, we use the quadratic formula on the derivative
	discriminant = (c)^2-(3d*b)
	if discriminant >= 0 #then the max may be  in [x0, x1]
		root1 = (-sqrt(discriminant)-c)/(3*d)
		root2 = (sqrt(discriminant)-c)/(3*d)
	end

	val1 = @evalpoly(x0, a, b, c, d)::T
	if discriminant >= 0
		val2 = @evalpoly(root1, a, b, c, d)::T
		val3 = @evalpoly(root2, a, b, c, d)::T
	else
		val2 = typemin(T)
		val3 = typemin(T)
	end
	val4 = @evalpoly(x1, a, b, c, d)::T

	if discriminant >= 0
		if root1 < x0 || root1 >x1
			val2 = typemin(T)
		end
		if root2 < x0 || root2 > x1
			val3 = typemin(T)
		end
	end

	theMax = max(val1, val2, val3, val4)

	maxPos = typemin(T)
	if val1 == theMax
		maxPos = x0
	elseif val2 == theMax
		maxPos = root1
	elseif val3 == theMax
		maxPos = root2
	elseif val4 == theMax
		maxPos = x1
	else
		println("Failed to find the max in maxcubic. This is probably a bug in Splines.jl")
		return nothing
	end

	return (theMax, maxPos)
end

function maxcubic{T}(x::T, a::T, b::T, c::T, d::T) #search on [0, x]
	return maxcubic(zero(T), x, a, b, c, d)
end

#FIXME: NaN if d = 0
function maxabscubic{T}(x0::T, x1::T, a::T, b::T, c::T, d::T) #a+bx+cx^2+dx^3
  @assert !isnan(x0)
  @assert !isnan(x1)
  @assert !isnan(a)
  @assert !isnan(b)
  @assert !isnan(c)
  @assert !isnan(d)
  #print("x0, x1, a, b, c, d: ", x0, a, b, c, d)
  #print("\n")
  if x0 == x1 #then there is no real max, but we do the sanest possible thing
    return (x0, @evalpoly(x0, a, b, c, d))
  end
  if d == 0 #then the max is at an endpoint. #FIXME: I don't think so
    endpt1_val = @evalpoly(x0, a, b, c, d)
    endpt2_val = @evalpoly(x1, a, b, c, d)
    return endpt2_val > endpt1_val ? (endpt2_val, x1) : (endpt1_val, x0)
  end

	#otherwise, we use the quadratic formula on the derivative
	discriminant = (c)^2-(3d*b)
	if discriminant >= 0 #then the max may be  in [x0, x1]
		root1 = (-sqrt(discriminant)-c)/(3*d)
		root2 = (sqrt(discriminant)-c)/(3*d)
	end

	val1 = abs(@evalpoly(x0, a, b, c, d)::T)
	if discriminant >= 0
		val2 = abs(@evalpoly(root1, a, b, c, d)::T)
		val3 = abs(@evalpoly(root2, a, b, c, d)::T)
	else
		val2 = typemin(T)
		val3 = typemin(T)
	end
	val4 = abs(@evalpoly(x1, a, b, c, d)::T)

	if discriminant >= 0
		if root1 < x0 || root1 >x1
			val2 = typemin(T)
		end
		if root2 < x0 || root2 > x1
			val3 = typemin(T)
		end
	end

	theMax = max(val1, val2, val3, val4)

  #	println(val1, val2, val3, val4, theMax)
	maxPos = typemin(T)
	if val1 == theMax
		maxPos = x0
	elseif val2 == theMax
		maxPos = root1
	elseif val3 == theMax
		maxPos = root2
	elseif val4 == theMax
		maxPos = x1
  #=	else
    println("Failed to find the max in maxabscubic. This is probably a bug in Splines.jl")
		return nothing=#
	end

	return (theMax, convert(T, maxPos))
end

function maxabscubic{T}(x::T, a::T, b::T, c::T, d::T) #search on [0, x]
	return maxcubic(zero(T), x, a, b, c, d)
end

function mincubic{T}(x0::T, x1::T, a::T, b::T, c::T, d::T) #a+bx+cx^2+dx^3
	if any_nan(x0, x1, a, b, c, d)
		throw(ArgumentError("Minimum of a cubic with NaN parameters is not defined"))
	end
  if x0 == x1
    return (x0, @evalpoly(x0, a, b, c, d))
  end
  if d == 0 #then the min is at an endpoint. #FIXME: I don't think so
    endpt1_val = @evalpoly(x0, a, b, c, d)
    endpt2_val = @evalpoly(x1, a, b, c, d)
    return endpt2_val < endpt1_val ? (endpt2_val, x1) : (endpt1_val, x0)
  end

	#otherwise, we use the quadratic formula on the derivative
	discriminant = (c)^2-(3*d*b)
	if discriminant >= 0
		root1 = (-sqrt(discriminant)-c)/(3*d)
		root2 = (sqrt(discriminant)-c)/(3*d)
	end

	val1 = @evalpoly(x0, a, b, c, d)
	if discriminant >= 0
		val2 = @evalpoly(root1, a, b, c, d)
		val3 = @evalpoly(root2, a, b, c, d)
	else
		val2 = typemax(T)
		val3 = typemax(T)
	end
	val4 = @evalpoly(x1, a, b, c, d)

	if discriminant >= 0
		if root1 < x0 || root1 > x1
			val2 = typemax(T)
		end

		if root2 < x0 || root2 > x1
			val3 = typemax(T)
		end
	end

	theMin = min(val1, val2, val3, val4)

	minPos = typemax(T)
	if val1 == theMin
		minPos = x0
	elseif val2 == theMin
	if !isdefined(:root1)
		println("val1: $val1 val2: $val2 val3: $val3 val4: $val4 theMin: $theMin")
	end
		minPos = root1
	elseif val3 == theMin
		minPos = root2
	elseif val4 == theMin
		minPos = x1
	else
		error("Failed to find the minimum. This is probably a bug in Splines.jl")
	end

	return (theMin, minPos)
end

function mincubic{T}(x::T, a::T, b::T, c::T, d::T) #search between 0 and x
	return mincubic(zero(T), x, a, b, c, d)
end

function integrate_cubic{T}(x0::T, x1::T, a::T, b::T, c::T, d::T) #a+bx+cx^2+dx^3
  return @evalpoly(x1, 0, a, b/2, c/3, d/4)-@evalpoly(x0, 0, a, b/2, c/3, d/4)
end


text:
###
# Optimization problems for VoPP Upper and Lower Bounds
###
#Throughout we assume mu = 1, c = 0, and hbar s.t. E[hbar(Vc)] = 0.
#Methods with aleading underscore are not true bounds but approximations for debugging/configuration

module vopp

using JuMP, Gurobi, Ipopt, LinearAlgebra, QuadGK, Roots

include("LambertW.jl")
include("Helpers.jl")
include("closedFormUB.jl")
#include("debugging.jl")
include("upperBounds.jl")
include("UpperBoundsSymmetric.jl")
include("UpperBoundsUnimodal.jl")
include("lowerBoundsSymmetric.jl")
include("lowerBoundsUnimodal.jl")

end #ends module


text:
using CalibrationErrorsDistributions
using Distances
using LinearAlgebra
using PDMats
using Random
using Test

using CalibrationErrorsDistributions:
    Wasserstein, SqWasserstein, MixtureWasserstein, SqMixtureWasserstein

Random.seed!(1234)

@testset "CalibrationErrorsDistributions" begin
    @testset "distances" begin
        @testset "Wasserstein distance" begin
            include("distances/wasserstein.jl")
        end
    end

    @testset "Normal" begin
        include("normal.jl")
    end
    @testset "MvNormal" begin
        include("mvnormal.jl")
    end

    @testset "deprecations" begin
        include("deprecated.jl")
    end
end


text:
#
# Edit `Port for Communication with the Julia Process` to 9999
#

include("port.jl")
using Atom
using Atom; using Juno; Juno.connect(port)


text:
module ManualModelMod

using DataFrames
using StatsBase
using ..DataMod

export ManualModel, predict

struct ManualModel <: RegressionModel
    coefs::Vector{Float64}
    intercept::Bool
    predicted::Symbol
    data::Data
    ManualModel(coefs::Vector{Float64}, intercept::Bool, predicted::Symbol, data::Data) =
        @inbounds new([abs(coefs[i]) > eps(Float64) ? coefs[i] : 0. for i in eachindex(coefs)], intercept, predicted, data)
end

ncoef(m::ManualModel) = sum(m.coefs .!= 0.)

function Base.show(io::IO, m::ManualModel)
    println(io, "Manual linear model predicting :" * string(m.predicted) * " -- " * (m.intercept ? "" : "not ") * "using intercept")
    show(io, CoefTable([m.coefs], ["Coef."], m.intercept ? ["Intercept"; ":" .* string.(m.data.xs)] : ":" .* string.(m.data.xs)))
end

StatsBase.coef(m::ManualModel)::Vector{Float64} = m.coefs
function StatsBase.predict(m::ManualModel, A::Matrix{Float64})
    col = size(A, 2)
    len = length(m.coefs)
    len != (m.intercept ? col + 1 : col) && error("Model size of $len does not match matrix size of $col")
    m.intercept && (A = hcat(ones(Float64, size(A, 1)), A))
    return A * m.coefs
end
StatsBase.predict(m::ManualModel, A::DataFrame) = StatsBase.predict(m, Matrix(A))
StatsBase.predict(m::ManualModel) = predict(m, m.data.xmat)

function StatsBase.loglikelihood(m::ManualModel)
    # This function is derived from the loglikelihood function in GLM.jl
    pred = predict(m)
    act = m.data.ydf[!, m.predicted]::Vector{Float64}
    n = nrow(m.data.df)
    dev = 0.0
    @inbounds @simd for i in 1:n
        dev += abs2(act[i] - pred[i])::Float64
    end
    return -n/2 * (log(2π * dev/n) + 1)
end
StatsBase.dof(m::ManualModel) = ncoef(m) + 1
StatsBase.nobs(m::ManualModel) = nrow(m.data.df)

end

using .ManualModelMod

max_stars_repo_name: KasperNooteboom/thesis-rvfl-fs

max_issues_repo_name: KasperNooteboom/thesis-rvfl-fs

max_forks_repo_name: KasperNooteboom/thesis-rvfl-fs


text:
# Intersection
# ------------

# There are two ways we can tackle intersection: traversal (search the tree for
# Segments that intersect) and iteration (iterate through the tree's Segments
# in sorted order).
#
# Traversal is roughly O(log(n)) while iteration is roughly O(n). Obviously
# if we are doing a single intersection, we would choose traversal, but when
# intersecting two trees it becomes more interesting.
#
# If the second tree is of size m, then intersection by traversal is
# O(min(m * log(n), n * log(m))) while interesction by iteration is O(m + n).
# When m and n are both large, it becomes more efficient to intersect by
# iteration, but deciding when requires some heuristics.
#

# Return true iff two key1 and key2 intersect.
@inline function intersects(key1::AbstractSegment{K}, key2::AbstractSegment{K}) where K
    return first(key1) <= last(key2) && first(key2) <= last(key1)
end

# Return true if the tree has an intersection with the given position
function hasintersection(t::SegmentBTree{K, V, B}, query) where {K, V, B}
    return hasintersection(t.root, convert(K, query))
end

function hasintersection(t::InternalNode{K, V, B}, query::K) where {K, V, B}
    if isempty(t) || t.maxend < query
        return false
    end

    for (i, child) in enumerate(t.children)
        if child.maxend >= query && (i == 1 || t.keys[i-1].first <= query)
            if hasintersection(child, query)
                return true
            end
        elseif minstart(child) > query
            break
        end
    end

    return false
end

function hasintersection(t::LeafNode{K, V, B}, query::K) where {K, V, B}
    if isempty(t) || t.maxend < query
        return false
    end

    for i in 1:length(t)
        if first(t.keys[i]) <= query <= last(t.keys[i])
            return true
        elseif query < first(t.keys[i])
            break
        end
    end

    return false
end

# Represent an intersection in an SegmentTree by pointing to a LeafNode
# and an index within that leaf node. No intersection is represented with
# index == 0.
mutable struct Intersection{K, V, B}
    index::Int
    node::LeafNode{K, V, B}

    Intersection{K,V,B}(index, node) where {K,V,B} = new{K,V,B}(index, node)
    Intersection{K,V,B}() where {K,V,B} = new{K,V,B}(0)
end

# Find the first Segment in the tree that intersects the query and return
# as a (leafnode, index) pair, indicating that leafnode.keys[index] intersects.
# If no intersection is found, index is 0 and leafnode is the last node
# searched.
function firstintersection!(t::SegmentBTree{K, V, B},
                            query::AbstractSegment{K},
                            lower::Union{Nothing, V},
                            out::Intersection{K, V, B},
                            filter::F) where {F, K, V, B}
    return firstintersection!(t.root, query, lower, out, filter)
end

function firstintersection!(t::InternalNode{K, V, B},
                            query::AbstractSegment{K},
                            lower::Union{Nothing, V},
                            out::Intersection{K, V, B},
                            filter::F) where {F, K, V, B}
    if isempty(t) || t.maxend < first(query)
        out.index = 0
        return
    end

    (query_first, query_last) = (first(query), last(query))

    i = lower === nothing ? 1 : searchsortedfirst(t.keys, notnothing(lower))

    @inbounds while i <= length(t.children)
        if i > 1 && t.keys[i-1].first > query_last
            break
        end
        if t.maxends[i] >= query_first
            firstintersection!(t.children[i], query, lower, out, filter)
            if out.index > 0
                return
            end
        end
        i += 1
    end

    out.index = 0
    return
end

function firstintersection!(t::LeafNode{K, V, B},
                            query::AbstractSegment{K},
                            lower::Union{Nothing, V},
                            out::Intersection{K, V, B},
                            filter::F) where {F, K, V, B}
    if isempty(t) || t.maxend < first(query)
        out.index = 0
        return
    end

    # TODO: search keys
    i = lower === nothing ? 1 :
        searchsortedfirst(t.entries, notnothing(lower), 1, Int(t.count),
                          Base.ord(basic_isless, identity, false))

    while i <= length(t)
        if intersects(t.entries[i], query) && filter(t.entries[i], query)
            out.index = i
            out.node = t
            return
        elseif last(query) < first(t.keys[i])
            break
        end
        i += 1
    end

    out.index = 0
    return
end

function firstintersection(t::InternalNode{K, V, B},
                           query::Segment{K},
                           lower::Union{Nothing, Segment{K}}) where {K, V, B}
    if isempty(t) || t.maxend < first(query)
        return (nothing, 1)
    end

    i = lower === nothing ? 1 : searchsortedfirst(t.keys, notnothing(lower))

    @inbounds while i <= length(t.children)
        if i > 1 && t.keys[i-1].first > last(query)
            break
        end
        if t.maxends[i] >= first(query)
            w, k = firstintersection(t.children[i], query, lower)
            w !== nothing && return notnothing(w), k
        end
        i += 1
    end

    return (nothing, 1)
end

function firstintersection(t::LeafNode{K, V, B}, query::Segment{K},
                           lower::Union{Nothing, Segment{K}}) where {K, V, B}

    if isempty(t) || t.maxend < first(query)
        return (nothing, 1)
    end

    i = lower === nothing ? 1 : searchsortedfirst(t.keys, notnothing(lower))

    while i <= length(t)
        if intersects(t.keys[i], query) &&
           (lower === nothing || !basic_isless(t.keys[i], notnothing(lower)))
            return (t, i)
        elseif last(query) < first(t.keys[i])
            break
        end
        i += 1
    end

    return (nothing, 1)
end

# Find the first key with an end point >= query
function firstfrom(t::SegmentBTree{K, V, B}, query::K) where {K, V, B}
    return firstfrom(t.root, query)
end

function firstfrom(t::InternalNode{K, V, B}, query::K) where {K, V, B}
    if isempty(t)
        return (t, 0)
    end

    for (i, child) in enumerate(t.children)
        if child.maxend >= query
            return firstfrom(child, query)
        end
    end

    return (t, 0)
end

function firstfrom(t::LeafNode{K, V, B}, query::K) where {K, V, B}
    if isempty(t)
        return (t, 0)
    end

    for i in 1:length(t)
        if last(t.keys[i]) >= query
            return (t, i)
        end
    end

    return (t, 0)
end

# If query intersects node.keys[i], return the next intersecting key as
# a (leafnode, index) pair.
function nextintersection!(t::LeafNode{K, V, B}, i::Integer,
                           query::AbstractSegment{K},
                           out::Intersection{K, V, B},
                           filter::F) where {F, K, V, B}
    j = i + 1
    while true
        while j <= length(t)
            if intersects(t.keys[j], query) && filter(t.entries[j], query)
                out.index = j
                out.node = t
                return
            end
            j += 1
        end
        j = 1

        t.right === nothing && break

        t = notnothing(t.right)

        if minstart(t) > last(query)
            break
        end

        # When we hit a leaf node with no possible intersections, we start over
        # from the top of the leaf. This is a heuristic to avoid linear search
        # behavior in the presence of extremely long Segments.
        if t.maxend < first(query)
            return firstintersection!(root(t), query,
                                      t.entries[t.count], out, filter)
        end
    end

    out.index = 0
    return
end

mutable struct SegmentIntersectionIterator{F, K, V, B}
    filter::F
    intersection::Intersection{K, V, B}
    t::SegmentBTree{K, V, B}
    query::AbstractSegment{K}

    function SegmentIntersectionIterator{F,K,V,B}(filter, intersection, t, query) where {F,K,V,B}
        return new{F,K,V,B}(filter, intersection, t, query)
    end

    function SegmentIntersectionIterator{F,K,V,B}() where {F,K,V,B}
        return new{F,K,V,B}(Intersection{F, K, V, B}())
    end
end

Base.eltype(::Type{SegmentIntersectionIterator{F,K,V,B}}) where {F,K,V,B} = V
Base.IteratorSize(::Type{SegmentIntersectionIterator{F,K,V,B}}) where {F,K,V,B} = Base.SizeUnknown()

# Intersect an Segment tree t with a single Segment, returning an iterator
# over the intersecting (key, value) pairs in t.
function Base.intersect(t::SegmentBTree{K, V, B}, query0::Tuple{Any, Any},
                        filter::F=true_cmp) where {F, K, V, B}
    query = Segment{K}(query0[1], query0[2])
    return intersect(t, query, filter)
end

function Base.intersect(t::SegmentBTree{K, V, B}, first::K, last::K,
                        filter::F=true_cmp) where {F, K, V, B}
    query = Segment{K}(first, last)
    return intersect(t, query, filter)
end

function Base.intersect(t::SegmentBTree{K, V, B}, query::AbstractSegment{K},
                        filter::F=true_cmp) where {F, K, V, B}
    return SegmentIntersectionIterator{F, K, V, B}(filter, Intersection{K, V, B}(), t, query)
end

function Base.iterate(it::SegmentIntersectionIterator{F, K, V, B}) where {F, K, V, B}
    firstintersection!(it.t, it.query, nothing, it.intersection, it.filter)
    return iterate(it, nothing)
end

function Base.iterate(it::SegmentIntersectionIterator{F, K, V, B}, _) where {F, K, V, B}
    intersection = it.intersection
    if intersection.index == 0
        return nothing
    end
    entry = intersection.node.entries[intersection.index]
    nextintersection!(intersection.node, intersection.index,
                      it.query, intersection, it.filter)
    return entry, nothing
end

mutable struct IntersectionIterator{F, K, V1, B1, V2, B2}
    filter::F

    t1::SegmentBTree{K, V1, B1}
    t2::SegmentBTree{K, V2, B2}

    # if true, use successive intersection, if false iterative.
    successive::Bool

    # true if done
    isdone::Bool

    # intersection state
    u::Union{Nothing, LeafNode{K, V1, B1}}
    v::Union{Nothing, LeafNode{K, V2, B2}}
    w::Union{Nothing, LeafNode{K, V2, B2}}
    i::Int
    j::Int
    k::Int

    function IntersectionIterator{F,K,V1,B1,V2,B2}(filter::F, t1, t2, successive::Bool) where {F,K,V1,B1,V2,B2}
        return new{F,K,V1,B1,V2,B2}(filter, t1, t2, successive)
    end
end

Base.eltype(::Type{IntersectionIterator{F,K,V1,B1,V2,B2}}) where {F,K,V1,B1,V2,B2} = Tuple{V1,V2}
Base.IteratorSize(::Type{IntersectionIterator{F,K,V1,B1,V2,B2}}) where {F,K,V1,B1,V2,B2} = Base.SizeUnknown()

function Base.iterate(it::IntersectionIterator, _=iterinitstate(it))
    if it.isdone
        return nothing
    end
    if it.successive
        u = it.u
        w = it.w
        value = (u.entries[it.i], w.entries[it.k])
        @nextleafkey(w, it.w, it.k)
        successive_nextintersection!(it)
    else
        u = it.u
        w = it.w
        value = (u.entries[it.i], w.entries[it.k])
        @nextleafkey(w, it.w, it.k)
        iterative_nextintersection!(it)
    end
    return value, nothing
end

function iterinitstate(it::IntersectionIterator)
    it.isdone = true
    if it.successive
        it.u = isempty(it.t1) ? nothing : firstleaf(it.t1)
        it.w = isempty(it.t2) ? nothing : firstleaf(it.t2)
        it.i, it.k = 1, 1
        successive_nextintersection!(it)
    else
        # Iterative Intersection: Intersect by iterating through the two
        # collections in unison.
        #
        # Notation notes. We proceed by finding all the the Segments
        # in t2 that intersect u.keys[i], before proceeding to the next
        # key in t1.
        #
        # We use v to keep track of the first leaf node in t2 that might contain an
        # intersecting Segment, # and with (w, j) the current intersecting entry in
        # t2. We need to hang onto v since we will need to jump back when the
        # entry in t1 gets incremented.
        #
        # The thing to keep in mind is that this is just like the merge operation in
        # mergesort, except that some backtracking is needed since intersection
        # isn't as simple as ordering.

        it.u = isempty(it.t1) ? nothing : firstleaf(it.t1)
        it.v = isempty(it.t2) ? nothing : firstleaf(it.t2)
        it.w = it.v
        it.i, it.j, it.k = 1, 1, 1
        iterative_nextintersection!(it)
    end
    return nothing  # iterator is stateful
end

function iterative_nextintersection!(
    it::IntersectionIterator{F, K, V1, B1, V2, B2}) where {F, K, V1, B1, V2, B2}

    u, v, w, i, j, k = it.u, it.v, it.w, it.i, it.j, it.k

    it.isdone = true

    while true
        u === nothing && return
        unode = notnothing(u)
        ukey = unode.keys[i]

        v === nothing && return
        vnode = notnothing(v)
        vkey = vnode.keys[j]

        # find next intersection w
        while w !== nothing
            wnode = notnothing(w)
            wkey = wnode.keys[k]

            if first(wkey) <= last(ukey)
                if first(ukey) <= last(wkey) && it.filter(unode.entries[i], wnode.entries[k])
                    it.isdone = false
                    break
                end
                @nextleafkey(wnode, w, k)
            else
                break
            end
        end

        # intersection was found
        if !it.isdone
            break
        end

        # if no intersection found, advance u and start again
        @nextleafkey(unode, u, i)

        u === nothing && break

        # find new v corresponding to new u
        unode = notnothing(u)
        ukey = unode.keys[i]
        while v !== nothing
            vnode = notnothing(v)
            vkey = vnode.keys[j]

            if last(vkey) < first(ukey)
                @nextleafkey(vnode, v, j)
            else
                break
            end
        end

        w, k = v, j
    end

    it.u = u
    it.v = v
    it.w = w
    it.i = i
    it.j = j
    it.k = k
    return
end

function successive_nextintersection!(
    it::IntersectionIterator{F, K, V1, B1, V2, B2}) where {F, K, V1, B1, V2, B2}

    u, w, i, k = it.u, it.w, it.i, it.k

    it.isdone = true

    while true
        u === nothing && return
        unode = notnothing(u)
        ukey = unode.keys[i]

        # find next intersection w
        while w !== nothing
            wnode = notnothing(w)
            wkey = wnode.keys[k]

            if first(wkey) <= last(ukey)
                if first(ukey) <= last(wkey) && it.filter(unode.entries[i], wnode.entries[k])
                    it.isdone = false
                    break
                end

                # next w key
                if k < length(wnode)
                    k += 1
                else
                    k = 1
                    w = wnode.right

                    # When we hit a leaf node with no possible intersections, we
                    # start over from the top of the leaf. This is a heuristic
                    # to avoid linear search behavior in the presence of
                    # extremely long Segments.
                    if w !== nothing && w.maxend < first(ukey)
                        wnode = notnothing(w)
                        w, k = firstintersection(it.t2.root, ukey,
                                                 wnode.keys[wnode.count])
                    end
                end
            else
                break
            end
        end

        # intersection was found
        if !it.isdone
            break
        end

        # if no intersection found, advance u and start again
        @nextleafkey(unode, u, i)

        u === nothing && break

        # find u's first intersection
        unode = notnothing(u)
        ukey = unode.keys[i]

        w, k = firstintersection(it.t2.root, ukey, nothing)
    end

    it.u = u
    it.w = w
    it.i = i
    it.k = k
    return
end

# Intersect two Segment trees, returning an iterator yielding values of the
# form:
#   ((key1, value1), (key2, value2))
#
# Where key1 is from the first tree and key2 from the second, and they
# intersect.
function Base.intersect(t1::SegmentBTree{K, V1, B1},
                        t2::SegmentBTree{K, V2, B2},
                        filter::F=true_cmp; method=:auto) where {F, K, V1, B1, V2, B2}
    # We decide heuristically which intersection algorithm to use.
    m = length(t1)
    n = length(t2)

    if method == :auto
        iterative_cost  = n + m
        successive_cost = 0.25 * m * log(1 + n) + 1e5
        if iterative_cost < successive_cost
            return IntersectionIterator{F, K, V1, B1, V2, B2}(filter, t1, t2, false)
        else
            return IntersectionIterator{F, K, V1, B1, V2, B2}(filter, t1, t2, true)
        end
    elseif method == :successive
        return IntersectionIterator{F, K, V1, B1, V2, B2}(filter, t1, t2, true)
    elseif method == :iterative
        return IntersectionIterator{F, K, V1, B1, V2, B2}(filter, t1, t2, false)
    else
        error("No such intersection method: $method")
    end
end


text:
export numberofstates, constants

abstract type AbstractBalanceLaw{FT, D, S, C} end

Base.eltype(::AbstractBalanceLaw{FT}) where {FT} = FT
Base.ndims(::AbstractBalanceLaw{FT, D}) where {FT, D} = D
numberofstates(::AbstractBalanceLaw{FT, D, S}) where {FT, D, S} = S
constants(::AbstractBalanceLaw{FT, D, S, C}) where {FT, D, S, C} = C

auxiliary(law::AbstractBalanceLaw, x⃗) = SVector(nothing)

function flux end
function wavespeed end
boundarystate(::AbstractBalanceLaw, n⃗, q⁻, aux⁻, tag) = q⁻, aux⁻
source!(::AbstractBalanceLaw, dq, q, aux) = nothing
nonconservative_term!(::AbstractBalanceLaw, dq, q, aux) = nothing

function entropy end
function entropyvariables end


text:
"""
    Lon <: XDim <: Dimension
    Lon(val=:)

Longitude [`Dimension`]($DDdimdocs).

## Example:
```julia
longdim = Lon(10:10:100)
# Or
val = A[Lon(1)]
# Or
mean(A; dims=Lon)
```
"""
@dim Lon XDim "Longitude"

"""
    Lat <: YDim <: Dimension
    Lat(val=:)

Latitude [`Dimension`]($DDdimdocs).

## Example:
```julia
vertdim = Lat(10:10:100)
# Or
val = A[Lat(1)]
# Or
mean(A; dims=Lat)
```
"""
@dim Lat YDim "Latitude"

"""
    Vert <: ZDim <: Dimension
    Vert(val=:)

Vertical [`Dimension`]($DDdimdocs).

## Example:
```julia
vertdim = Vert(10:10:100)
# Or
val = A[Vert(1)]
# Or
mean(A; dims=Vert)
```
"""
@dim Vert ZDim "Vertical"

"""
    Band <: Dimension
    Band(val=:)

Band [`Dimension`]($DDdimdocs) for multi-band rasters.

## Example:
```julia
banddim = Band(10:10:100)
# Or
val = A[Band(1)]
# Or
mean(A; dims=Band)
```
"""
@dim Band

"""
    mappedbounds(x)

Get the bounds converted to the `usercrs` value.

Whithout ArchGDAL loaded, this is just the regular bounds.
"""
function mappedbounds end

mappedbounds(dims::Tuple) = map(mappedbounds, dims)
mappedbounds(dim::Dimension) = bounds(dim)
mappedbounds(dim::Union{Lat,Lon}) = mappedbounds(mode(dim), dim)
mappedbounds(::Mapped, dim) = bounds(dim)
@noinline mappedbounds(mode::IndexMode, dim) =
    if mode isa Projected
        error("Load ArchGDAL to convert Projected mode bounds to mapped")
    else
        error("cannot get mapped bounds of a $(nameof(typeof(mode))) mode dim")
    end

projectedbounds(dims::Tuple) = map(projectedbounds, dims)
projectedbounds(dim::Dimension) = bounds(dim)
projectedbounds(dim::Union{Lat,Lon}) = projectedbounds(mode(dim), dim)
projectedbounds(::Projected, dim) = bounds(dim)
@noinline projectedbounds(mode::IndexMode, dim) =
    if mode isa Mapped
        error("Load ArchGDAL to convert Mapped mode dim to projected")
    else
        error("cannot get projected bounds of a $(nameof(typeof(mode))) mode dim")
    end

"""
    mappedindex(x)

Get the index value of a dimension converted to the `usercrs` value.

Whithout ArchGDAL loaded, this is just the regular dim value.
"""
function mappedindex end

mappedindex(dims::Tuple) = map(mappedindex, dims)
mappedindex(dim::Dimension) = index(dim)
mappedindex(dim::Union{Lat,Lon}) = mappedindex(mode(dim), dim)
mappedindex(::Mapped, dim) = index(dim)
@noinline mappedindex(mode::IndexMode, dim) =
    if mode isa Projected
        error("Load ArchGDAL to convert Projected mode index to mapped")
    else
        error("cannot get mapped index of a $(nameof(typeof(mode))) mode dim")
    end

projectedindex(dims::Tuple) = map(projectedindex, dims)
projectedindex(dim::Dimension) = index(dim)
projectedindex(dim::Union{Lat,Lon}) = projectedindex(mode(dim), dim)
projectedindex(::Projected, dim) = index(dim)
@noinline projectedindex(mode::IndexMode, dim) =
    if mode isa Mapped
        error("Load ArchGDAL to convert Mapped mode index to projected")
    else
        error("cannot get projected index of a $(nameof(typeof(mode))) mode dim")
    end


text:
const TINY_HANDMADE_INSTANCES = Tuple{String, Int64, String}[
	("Plate 1x1 and 1 piece 1x1 with profit 1.", 1, """
		1 1
		1
		1 1 1 1
	""")
	, ("Pieces do not fit plate.", 0, """
		10 10
		3
		15 15 1 1
		15 10 1 1
		10 15 1 1
	""")
	, ("One piece fit the plate others do not.", 10, """
		10 10
		4
		10 10 10 1
		15 15 1 1
		15 10 1 1
		10 15 1 1
	""")
	, ("Two halfs instead slightly larger very efficient piece.", 20, """
		10 20
		2
		10 10 10 2
		10 11 19 1
	""")
	, ("The whole plate could be used if the plates were rotated.", 20, """
		16 20
		2
		10 8 10 2
		10 16 19 1
	""")
	, ("Plate 100x100 and 100 pieces 10x10 with profit 7.", 700, """
		100 100
		1
		10 10 7 100
	""")
	, ("Plate 100x100 and 50 pieces 10x10 with profit 7.", 350, """
		100 100
		1
		10 10 7 50
	""")
	, ("Plate 50x200 and 200 pieces 10x10 with profit 7.", 700, """
		100 100
		1
		10 10 7 200
	""")
]

const EASY_LITERATURE_INSTANCES = Tuple{String, Int64, String}[
  ("gcut1", 48368, """
		250 250
		10
		167 184 30728 1
		114 118 13452 1
		167 152 25384 1
		83 140 11620 1
		70 86 6020 1
		143 166 23738 1
		120 160 19200 1
		66 148 9768 1
		87 141 12267 1
		69 165 11385 1
	""")
	, ("cgcut1", 244, """
		15 10
		7
		8 4 66 2
		3 7 35 1
		8 2 24 3
		3 4 17 5
		3 3 11 2
		3 2  8 2
		2 1  2 1
	""")
	, ("CHL5", 390, """
		20 20
		10
		14  2 28 1
		1   5  5 1
		20  4 80 2
		12  3 36 3
		11  8 88 2
		11  6 66 2
		7   9 63 1
		17  5 85 1
		7  14 98 2
		1   7  7 3
	""")
]

import JuMP
import MathOptInterface
const MOI = MathOptInterface
import GLPK # to load glue code in GuillotineModels
import Cbc
using GuillotineModels: build_model, get_cut_pattern, to_pretty_str, simplify!
using GuillotineModels.Data: read_from_string
using GuillotineModels.Utilities.Args: accepted_arg_list,
	create_normalized_arg_subset
using GuillotineModels.CommandLine.SolversArgs: empty_configured_model

function test_obj_val(
	model :: Symbol, solver :: Symbol, instance :: String, obj_val :: Int64
)
	instance = read_from_string(Val(:Classic_G2KP), instance)
	N, L, W, l, w, p, d = length(instance.l), instance.L, instance.W,
		instance.l, instance.w, instance.p, instance.d

	solver_type = Val(solver)
	solver_pp = create_normalized_arg_subset(
		Dict{String, Any}("no-output" => true),
		accepted_arg_list(solver_type)
	)
	m = empty_configured_model(solver_type, solver_pp)
	#println(instance)
	bmr = build_model(Val(model), m, d, p, l, w, L, W, Dict{String, Any}())
	#JuMP.write_to_file(m, "Cbc_false_infeasible.lp")
	#JuMP.write_to_file(m, "Cbc_false_infeasible.mof.json")
	#JuMP.write_to_file(m, "Cbc_false_infeasible.mps")
	# TODO: uncomment the real tests below after finishing the pricing.
	JuMP.optimize!(m)
	@test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@test JuMP.objective_value(m) ≈ obj_val rtol=1e-6 atol=1e-6
	# #=
	if model === :PPG2KP
		solution = get_cut_pattern(Val(model), m, eltype(d), eltype(l), bmr)
		#@show solution
		#println("to_pretty_str(solution) = $(to_pretty_str(solution))")
		#println("to_pretty_str(simplify!(solution)) = $(to_pretty_str(simplify!(solution)))")
	end
	# =#
end

function test_obj_val_of_all_combinations(
	models :: Vector{Symbol},
	solvers :: Vector{Symbol},
	instances :: Vector{Tuple{String, Int64, String}}
)
	@testset "Model: $model" for model in models
		@testset "Solver: $solver" for solver in solvers
			for (description, obj_val, instance) in instances
				@testset "Instance: $description" begin
					#println(description)
					test_obj_val(model, solver, instance, obj_val)
				end
			end
		end
	end
end

test_obj_val_of_all_combinations(
	[:PPG2KP, #=:Flow=#],
	[:GLPK, #=:Cbc=#], # Cbc is kinda bugged
	[TINY_HANDMADE_INSTANCES; EASY_LITERATURE_INSTANCES]
)


text:
# Faster to not load precompiled
#VERSION >= v"0.4.0-dev+6521" && __precompile__()

module ZChop
using Compat
export zchop, zchop!

const zeps = 1e-14

zchop{T<:Real}(x::T, eps=zeps) = abs(x) > convert(T,eps) ? x : zero(T)
zchop{T<:Integer}(x::T, eps=zeps) = abs(x) > eps ? x : zero(T)
zchop{T<:Complex}(x::T, eps=zeps) = complex(zchop(real(x),eps),zchop(imag(x),eps))
zchop{T<:AbstractArray}(a::T, eps=zeps) = (b = similar(a); for i in 1:length(a) b[i] = zchop(a[i],eps) end ; b)
zchop!{T<:AbstractArray}(a::T, eps=zeps) = (for i in 1:length(a) a[i] = zchop(a[i],eps) end ; a)
zchop{T<:Union{AbstractString,Char}}(x::T,eps=zeps) = x
zchop{T<:Irrational}(x::T,eps=zeps) = zchop(float(x),eps)
zchop(x::Expr,eps=zeps) = Expr(x.head,zchop(x.args)...)
zchop(x,eps) =  applicable(start,x) ? map((x)->zchop(x,eps),x) : x
zchop(x) = applicable(start,x) ? map(zchop,x) : x

end # module ZChop


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

using Markdown
using InteractiveUtils

# ╔═╡ 2a641766-809f-11eb-274a-8fdb234c0d08
begin
  using Pkg
  Pkg.activate("MLJ_env", shared=true)

  using CSV
  using DataFrames
  using PlutoUI
  using Test
  using Printf
  using Plots
  using LinearAlgebra
  using SparseArrays
  using TextAnalysis
  using Random
  using Distances
  using Statistics
  using Dates
end

# ╔═╡ f701bf4a-809e-11eb-2f1f-bb7dad75bfcc
md"""
## k-means with text data

In this assignment we will
  - Cluster Wikipedia documents using k-means
  - Explore the role of random initialization on the quality of the clustering
  - Explore how results differ after changing the number of clusters
  - Evaluate clustering, both quantitatively and qualitatively

When properly executed, clustering uncovers valuable insights from a set of unlabeled documents.
"""

# ╔═╡ 692a39c8-80a4-11eb-2077-e9f391787416
# Pkg.add("TextAnalysis")
# Pkg.add("Distances")

# ╔═╡ 5b3b148e-809f-11eb-3631-cf4682cbd959
begin
  """compute norm of a sparse vector"""
  norm_sv(v::AbstractVector) = sqrt.(dot(v, v'))

  function norm_sv(m::AbstractMatrix)
    n = size(m)[1]
    reshape([norm_sv(m[ix, :]) for ix in 1:n], n, 1)
  end
end

# ╔═╡ 97d6cb18-809f-11eb-075a-7b3c5354ce4c
## Norm of a sparse vector
begin
  v = sparsevec([1 2 0 0 0 3 0 4 0])
  nv = norm_sv(v)

  @test issparse(v)
  @test nv ≈ √(1^2 + 2^2 + 3^2 + 4^2)
end

# ╔═╡ 97d2bd70-8239-11eb-2a8c-efea2584e6d6
## Norm of a dense Matrix (row only)
begin
  m = [1 2 3; 4 1 2; 1 3 3; 1 0 2]  # 4×3
  dim = size(m)[1]

  norm_sv(m)

  ## then normalize matrix
  norm_m = m ./ norm_sv(m)

  ## test all rows have norm 1
  @test all(v -> norm_sv(v) ≈ 1., eachrow(norm_m))
end

# ╔═╡ 54ce09fa-8246-11eb-35fc-a9ba218d8c76
## Norm of a sparse Matrix (row only)
begin
  mₛ = sparse([1 2 0 3; 4 0 0 0; 0 0 3 1; 1 0 2 0])  # 4×4
  dimₛ = size(mₛ)[1]

  norm_sv(mₛ)

  ## then normalize matrix
  norm_mₛ = mₛ ./ norm_sv(mₛ)

  ## test all rows have norm 1
  @test all(v -> norm_sv(v) ≈ 1., eachrow(norm_mₛ))
end

# ╔═╡ 6ada72c6-80a0-11eb-2346-c3554265a4de
md"""
#### Load in the Wikipedia dataset

**Just a small subset of the data to validate our implementation.**

**Later we will scale and see...**
"""

# ╔═╡ 9ee529a8-80a0-11eb-23cc-b78ead6503e4
begin
  # just a subset
  wiki = train = CSV.File("../../ML_UW_Spec/data/people_wiki_500samples.csv";
                          header=true) |> DataFrame;
  size(wiki)
end

# ╔═╡ 16431f28-80a4-11eb-12a2-f3104b93f673
describe(wiki, :eltype, :nmissing, :first => first)

# ╔═╡ e98e9a5e-80a3-11eb-1d38-db4ff0f9575b
md"""
Let us add a unique id to each row of the wiki DataFrame
"""

# ╔═╡ fbebf142-80a3-11eb-0913-4b43d543a681
insertcols!(wiki, 1, :id => collect(1:size(wiki)[1]));

# ╔═╡ 637c61ac-80a4-11eb-21be-f1661d34b568
names(wiki)

# ╔═╡ a30ee624-81e8-11eb-0c75-a5715165ff36
md"""
#### Extract features
"""

# ╔═╡ 4543b590-80a5-11eb-3215-b9bb37e1e766
begin
  crps = Corpus(StringDocument.(wiki.text))
  update_lexicon!(crps)
end

# ╔═╡ fd36c16e-80a6-11eb-3a80-5d7b8cb48b38
dtm = DocumentTermMatrix(crps);

# ╔═╡ 4e94e7a2-80a7-11eb-0fb3-1955f645d0ca
ya_tf_idf = tf_idf(dtm);

# ╔═╡ 107cca2e-80ad-11eb-19b8-bf015a43b184
lexicon_size(crps)

# ╔═╡ 69c34f00-80ac-11eb-0b03-0362d6f5e84e
crps

# ╔═╡ b0d25bf0-80ae-11eb-00de-d57e0d22a000
dtm[1, :], ya_tf_idf[1, :]

# ╔═╡ 07757e58-8247-11eb-299a-97cbb8ce1460
md"""
As discussed in the previous assignment, Euclidean distance can be a poor metric of similarity between documents, as it unfairly penalizes long articles. For a reasonable assessment of similarity, we should disregard the length information and use length-agnostic metrics, such as cosine distance.

The k-means algorithm does not directly work with cosine distance, so we take an alternative route to remove length information: we normalize all vectors to be unit length. It turns out that Euclidean distance closely mimics cosine distance when all vectors are unit length. In particular, the squared Euclidean distance between any two vectors of length one is directly proportional to their cosine distance.

We can prove this as follows. Let $\mathbf{x}$ and $\mathbf{y}$ be normalized vectors, i.e. unit vectors, so that $\|\mathbf{x}\|=\|\mathbf{y}\|=1$. Write the squared Euclidean distance as the dot product of $(\mathbf{x} - \mathbf{y})$ to itself:

```math
\begin{align*}
\|\mathbf{x} - \mathbf{y}\|^2 &= (\mathbf{x} - \mathbf{y})^T(\mathbf{x} - \mathbf{y})\\
                              &= (\mathbf{x}^T \mathbf{x}) - 2(\mathbf{x}^T \mathbf{y}) + (\mathbf{y}^T \mathbf{y})\\
                              &= \|\mathbf{x}\|^2 - 2(\mathbf{x}^T \mathbf{y}) + \|\mathbf{y}\|^2\\
                              &= 2 - 2(\mathbf{x}^T \mathbf{y})\\
                              &= 2(1 - (\mathbf{x}^T \mathbf{y}))\\
                              &= 2\left(1 - \frac{\mathbf{x}^T \mathbf{y}}{\|\mathbf{x}\|\|\mathbf{y}\|}\right)\\
                              &= 2\left[\text{cosine distance}\right]
\end{align*}
```

This tells us that two **unit vectors** that are close in Euclidean distance are also close in cosine distance. Thus, the k-means algorithm (which naturally uses Euclidean distances) on normalized vectors will produce the same results as clustering using cosine distance as a distance metric.

"""

# ╔═╡ 99a6e9c4-8238-11eb-3acc-63a515a49b64
norm_ya_tf_idf = ya_tf_idf ./ norm_sv(ya_tf_idf);

# ╔═╡ 3f2f94e4-80af-11eb-32d2-a5e101797f53
size(ya_tf_idf), lexicon(crps)

# ╔═╡ fcfefb98-81e7-11eb-014a-b72b0477acab
# FIXME: disabled for now
# begin
#   update_inverse_index!(crps);
#   map_index_to_word = inverse_index(crps)
# end

# ╔═╡ 6c7e5bb4-823b-11eb-3833-254b4f5dafa7
## all rows are unit vector (vector with unit norm)
@test all(v -> norm_sv(v) ≈ 1., eachrow(norm_ya_tf_idf))

# @test all(ix -> abs(norm_sv(norm_ya_tf_idf[ix, :]) - 1.) ≤ 1e-8, collect(size(ya_tf_idf)[1]))

# ╔═╡ aaea9810-81e9-11eb-0452-83d979dca4c3
md"""
### Implement k-means

Let us implement the k-means algorithm. First, we choose an initial set of centroids. A common practice is to choose randomly from the data points.

Note: In practice, we highly recommend to use different seeds every time (for instance, by using the current timestamp).
"""

# ╔═╡ aad0780c-81e9-11eb-2344-63996006b0cb
function get_init_centroids(df, k; seed=42)
  """
  Randomly choose k data points as initial centroids
  """
    rng = MersenneTwister(seed);
    n = size(df)[1]
    rand_ix = randperm(rng, n)[1:k]  ## Pick k indices from range [1, n].

    ## Keep centroids as dense format, as many entries will be nonzero
    ## due to averaging.
    ## As long as at least one document in a cluster contains a word,
    ## it will carry a nonzero weight in the TF-IDF vector of the centroid.
    Matrix(df[rand_ix, :])
end

# ╔═╡ 446d44fa-81eb-11eb-0daa-83db97d8eec9
md"""
After initialization, the k-means algorithm iterates between the following two steps:

  1. Assign each data point to the closest centroid.
$$z_i \gets \mathrm{argmin}_j \|\mu_j - \mathbf{x}_i\|^2$$

  2. Revise centroids as the mean of the assigned data points.
$$\mu_j \gets \frac{1}{n_j}\sum_{i:z_i=j} \mathbf{x}_i$$

In pseudocode, we iteratively do the following:
```
cluster_assignment = assign_clusters(data, centroids)
centroids = revise_centroids(data, k, cluster_assignment)
```
"""

# ╔═╡ f1494764-81eb-11eb-02e8-e9a44cbbd905
md"""
#### Assigning clusters

How do we implement Step 1 of the main k-means loop above?

For the sake of demonstration, let's look at documents 100 through 102 as query documents and compute the distances between each of these documents and every other document in the corpus.

In the k-means algorithm, we will have to compute pairwise distances between the set of centroids and the set of documents.

"""

# ╔═╡ 444fc56a-81eb-11eb-27d2-1b249cae06d9
begin
    queries = norm_ya_tf_idf[100:102, :] ## Get the TF-IDF vectors for documents 100 through 102.
    ## Compute pairwise distances from every data point to each query vector.
    dist = pairwise(Euclidean(), norm_ya_tf_idf, queries, dims=1)
    @assert size(dist) == (size(norm_ya_tf_idf)[1], size(queries)[1])
end

# ╔═╡ 9716ee8a-81fb-11eb-3279-c19f1a5736a2
closest_cluster = argmin(dist; dims=2)

# ╔═╡ 95d6fcdc-81ff-11eb-0d5b-5b161dba559f
@test length(closest_cluster) == size(dist)[1]

# ╔═╡ 68962f52-8215-11eb-1cf2-1763b3ea2fcf
const ϵ = 1e-9 # tolerance for distances

# ╔═╡ aab3edea-81e9-11eb-2104-f18b12260e09
function assign_clusters(data, centroids)
    ## Compute distances between each data point and the set of centroids:
    dists_from_centroids = pairwise(Euclidean(ϵ), data, centroids, dims=1)

    ## Compute cluster assignments for each data point:
    argmin(dists_from_centroids; dims=2)
end

# ╔═╡ aa994aa0-81e9-11eb-1fca-f1b45c0220f0
md"""
**Checkpoint.** Let us check if Step 1 was implemented correctly.

With rows 1, 3, 5, and 7 of corpus as an initial set of centroids, we assign cluster labels to rows 30, 37, ..., and 99 of corpus.
"""

# ╔═╡ 8f71f718-8200-11eb-2dc1-4bc972d149f2
begin
  sample_data = norm_ya_tf_idf[50:4:99, :]
  sample_centroids = norm_ya_tf_idf[31:3:39, :]

  @assert size(sample_centroids)[1] == 3

  sample_cluster_assign = assign_clusters(sample_data, sample_centroids)
end

# ╔═╡ 20c45328-8237-11eb-1cf3-a9e5a4ee496f
size(sample_centroids), size(sample_data)

# ╔═╡ f812155c-8232-11eb-0113-5f5c9c8044e2
dists_from_centroids = pairwise(Euclidean(ϵ), sample_data, sample_centroids, dims=1)

# ╔═╡ 19ac0ab6-8212-11eb-27ed-132248a13585
size(sample_centroids), size(sample_data)

# ╔═╡ 756d72ec-8208-11eb-378d-2f9b2b288791
md"""
#### Revising clusters

Let's turn to Step 2, where we compute the new centroids given the cluster assignments.

To develop intuition about filtering, let's look at a toy example consisting of 3 data points and 2 clusters.

"""

# ╔═╡ 9a0f5d4a-8208-11eb-11d0-9742f32e4705
begin
  toy_data = [1. 2. 0.; 0. 0. 0.; 2. 2. 0.]
  toy_centroids = [0.5 0.5 0.; 0. -0.5 0.]

  size(toy_data), size(toy_centroids)
end

# ╔═╡ 23ea71ba-822e-11eb-0086-493cd5a4ae95
pairwise(Euclidean(ϵ), toy_data, toy_centroids, dims=1)

# ╔═╡ 4f06ce52-822e-11eb-1c3f-5777aee2de37
pairwise(Euclidean(ϵ), sparse(toy_data), sparse(toy_centroids), dims=1)

# ╔═╡ e58d0b00-8208-11eb-0414-019ec7c74eec
toy_cluster_assign = assign_clusters(toy_data, toy_centroids)

# => 1st and 3rd datapoint assigned to cluster defined by 1st centroid
# => 2nd datapoint assigned to cluster defined by 2nd centroid

# ╔═╡ f00143fa-820d-11eb-2983-43333f967308
function find_cluster(data, cluster_assign, cluster_num)
  """
  Return datapoints assigned to cluster cluster_num
  """
  ary = map(c -> data[c[1], :],
    filter(c -> c[2] == cluster_num,  cluster_assign)
  )
  Matrix(hcat(ary...)')
end

# ╔═╡ 6ceb9cf8-820c-11eb-1652-99c659469ff9
## find datapoint in 1st cluster
find_cluster(toy_data, toy_cluster_assign, 1)

# ╔═╡ 13440200-8250-11eb-0026-c91c8fdff383
## how many datapoints?
size(find_cluster(toy_data, toy_cluster_assign, 1))[1]

# ╔═╡ 86d3d792-8250-11eb-36c0-33630e0c38ba
function count_cluster_assign(cluster_assign)
	"""
	How many data points per cluster
	"""
	ids = map(c -> c[2], cluster_assign)
	[length(ids[ids .== id]) for id ∈ sort(unique(ids))]
end

# ╔═╡ 509e66ca-8252-11eb-28da-dd665ec5575a
count_cluster_assign(sample_cluster_assign)

# ╔═╡ 58c012f2-8252-11eb-3b9a-6b736099ad9d
count_cluster_assign(toy_cluster_assign)

# ╔═╡ c5137efa-820c-11eb-3773-057e25e6c242
## find datapoint in 2nd cluster
find_cluster(toy_data, toy_cluster_assign, 2)

# ╔═╡ 91866190-8213-11eb-00f2-036a1e31e428
## using sparse matrix -> cluster 3
begin
  _s = find_cluster(sample_data, sample_cluster_assign, 3)
  @assert size(_s)[1] ≥ 0

  size(_s), _s
end

# ╔═╡ 33b957ac-820e-11eb-17d6-a35b4a60fb2d
md"""
Given all the data points in a cluster, it only remains to compute the mean.

Use `mean()` from `Statistics` package.
"""

# ╔═╡ 338394f8-820e-11eb-0db4-1f86c79eacbc
mean(find_cluster(toy_data, toy_cluster_assign, 1), dims=1)

# ╔═╡ ab856686-8213-11eb-17e5-730b0a96ad12
## change cluster_num
mean(find_cluster(sample_data, sample_cluster_assign, 3), dims=1)

# ╔═╡ 14fc4e02-8211-11eb-3f9d-d37995fb2522
function revise_centroids(data, k, cluster_assign)
  	n_centroids = []
  	n = size(data)[2]

  	for ix ∈ 1:k
		## 1. Select all data points that belong to cluster ix
        datapts_clu_ix = find_cluster(data, cluster_assign, ix)

    	# size(datapts_clu_ix)[1] == 0 && continue

    	## 2. Compute the mean of the data points.
        centroid = mean(datapts_clu_ix, dims=1)

    	push!(n_centroids, centroid)
  	end
  	reshape(vcat(n_centroids...), k, :)
end

# ╔═╡ 7d692dae-821c-11eb-2aef-97b108c8d8a5
begin
  xv = get_init_centroids(ya_tf_idf, 3; seed=42)
  size(xv), typeof(xv), xv
end

# ╔═╡ c6bb718c-821f-11eb-072c-bf89e4b16e6a
size(toy_data), size(sample_data)

# ╔═╡ e51692fe-8208-11eb-1fc9-1f8f3a9d5db6
begin
  vs1 = revise_centroids(toy_data, 2, toy_cluster_assign)
  size(vs1), typeof(vs1), vs1
end

# ╔═╡ 9330b2f2-8211-11eb-205e-7db725dfdf85
begin
  vs2 = revise_centroids(sample_data, 3, sample_cluster_assign)
  size(vs2), typeof(vs2), vs2
end

# ╔═╡ 92f9c26c-8211-11eb-1f71-818d8fd916fb
md"""
#### Assessing the convergence

How can we tell if the k-means algorithm is converging?

We can look at the cluster assignments and see if they stabilize over time. In fact, we'll be running the algorithm until the cluster assignments stop changing at all.

To be extra safe, and to assess the clustering performance, we'll be looking at an additional criteria: *the sum of all squared distances between data points and centroids*.

This is defined as:

$$J(\mathcal{Z},\mu) = \sum_{j=1}^k \sum_{i:z_i = j} \|\mathbf{x}_i - \mu_j\|^2.$$

The smaller the distances, the more homogeneous the clusters are. In other words, we'd like to have "tight" clusters.
"""

# ╔═╡ 92e15d44-8211-11eb-2977-8d3b271553e1
function compute_heterogeneity(data, k, centroids, cluster_assign)
    heterogeneity = zero(eltype(data))

    for ix ∈ 1:k
    	## Select all data points that belong to cluster ix
        datapts_clu_ix = find_cluster(data, cluster_assign, ix)

    	## check if i-th cluster is non-empty
        if size(datapts_clu_ix)[1] > 0
      		## Compute distances from centroid to data points
      		m_centroid_ix = reshape(centroids[ix, :], 1, :)
        	dists = pairwise(Euclidean(ϵ), datapts_clu_ix, m_centroid_ix, dims=1)
        	heterogeneity += sum(dists .^ 2)
    	end
  	end
    heterogeneity
end

# ╔═╡ 92c7d752-8211-11eb-1546-cd510e3d8659
compute_heterogeneity(sample_data, 2, sample_centroids, sample_cluster_assign)

# ╔═╡ 92accbd8-8211-11eb-0703-0bb21268d6bc
md"""
#### Combining into a single function

Once the two k-means steps have been implemented, as well as our heterogeneity metric we wish to monitor, it is only a matter of putting these functions together to write a k-means algorithm that

  -  Repeatedly performs Steps 1 and 2
  - Tracks convergence metrics
  -  Stops if either no assignment changed or we reach a certain number of iterations.
"""

# ╔═╡ c43df0a0-8216-11eb-3be2-0f69c1faaba2
function kmeans(data, k, init_centroids;
    maxiter=100, record_heterogeneity=nothing, verbose=false)
    """
  This function runs k-means on given data and initial set of centroids.

  - maxiter: maximum number of iterations to run.
    - record_heterogeneity: (optional) a list, to store the history of heterogeneity as function of iterations if nothing, do not store the history.
    - verbose: if true, print how many data points changed their cluster labels in each iteration
  """
    centroids = copy(init_centroids)
    prev_cluster_assign, cluster_assign = nothing, nothing

    for ix ∈ 1:maxiter
        verbose && print("\niter: $(ix) ")

        ## 1. Make cluster assignments using nearest centroids
        cluster_assign = assign_clusters(data, centroids)

    	## 2. Compute a new centroid for each of the k clusters,
    	##    averaging all data points assigned to that cluster.
        centroids = revise_centroids(data, k, cluster_assign)

    	## 3. Check for convergence: if none of the assignments changed, stop
        if !isnothing(prev_cluster_assign) &&
          all(t -> t[1] == t[2], zip(prev_cluster_assign, cluster_assign))
      		break
    	end

    	## 4. Print number of new assignments
        if !isnothing(prev_cluster_assign)
            num_changed = sum(prev_cluster_assign .≠ cluster_assign)
            verbose &&
			  @printf("\t%5d elements changed their cluster assignment.\n", num_changed)
    	end

        if !isnothing(record_heterogeneity)
      		## 5. Record heterogeneity convergence metric
            score = compute_heterogeneity(data, k, centroids, cluster_assign)
            push!(record_heterogeneity, score)
    	end

        prev_cluster_assign = copy(cluster_assign) ## Copy (and not a ref.)
  	end
    (centroids, cluster_assign)
end

# ╔═╡ c422dff2-8216-11eb-2f63-235d58c77e35
md"""
##### Plotting the convergence
"""

# ╔═╡ c406c684-8216-11eb-3aba-65c9a55b7441
function plot_heterogeneity(heterogeneity, k)
    # plt.figure(figsize=(7,4))
 	plot(heterogeneity,
    	linewidth=4,
    	xlabel="# Iterations",
    	ylabel="Heterogeneity",
    	title="Heterogeneity of clustering over time, K=$(k)"
 	)
end

# ╔═╡ c3ebe820-8216-11eb-06e1-538f0cda84b7
begin
  k, heterogeneity = 3, []
  init_centroids = get_init_centroids(norm_ya_tf_idf, k; seed=42)
  centroids, cluster_assign = nothing, nothing

  with_terminal() do
    global centroids, cluster_assign = kmeans(norm_ya_tf_idf, k, init_centroids;
        maxiter=400, record_heterogeneity=heterogeneity, verbose=true);
  end
end

# ╔═╡ 31a82b96-8234-11eb-2795-bbadf977d8a3
centroids

# ╔═╡ bf1fd5b6-8234-11eb-0fbc-91da7f6dbe28
heterogeneity

# ╔═╡ 3633fa2a-8234-11eb-315d-d93094e52cab
size(find_cluster(norm_ya_tf_idf, cluster_assign, 1))

# ╔═╡ 67ecd7b2-8234-11eb-03f5-41fd8adda1fd
size(find_cluster(norm_ya_tf_idf, cluster_assign, 2))

# ╔═╡ 81b59776-8234-11eb-03d8-9736ea84e866
size(find_cluster(norm_ya_tf_idf, cluster_assign, 3))

# ╔═╡ a6d93074-8234-11eb-3cba-f1d9d8dbdfe3
plot_heterogeneity(heterogeneity, k)

# ╔═╡ 30eb64d4-824a-11eb-06c6-3393ca6a5780
md"""
#### Beware of local maxima

One weakness of k-means is that it tends to get stuck in a local minimum. To see this, let us run k-means multiple times, with different initial centroids created using different random seeds.

Note: Again, in practice, you should set different seeds for every run. We give you a list of seeds for this assignment so that everyone gets the same answer.

This may take several minutes to run.
"""

# ╔═╡ 2c29b028-8253-11eb-2781-f3de2eb49ba3
now()

# ╔═╡ 4c3ed812-824a-11eb-356a-6594a7b2a65e
function run_kmeans_with_multi_seeds(tfidf; k=10, init=:random)
 	heterogeneity = Dict()
  	cluster_assign_h = Dict()
  	# start = now()
  	m_ix_largest_clu, m_larg_clu_sz = -1, -1

  	for seed ∈ [0, 2000, 4000, 6000, 8000, 10000, 12000]
    	init_centroids = init == :random ? get_init_centroids(tfidf, k; seed) :
			smart_initialize(tfidf, k, seed)

    	centroids, cluster_assign = kmeans(tfidf, k, init_centroids;
			maxiter=400, record_heterogeneity=nothing, verbose=false)

		## To save time, compute heterogeneity only once in the end
    	heterogeneity[seed] = compute_heterogeneity(tfidf, k, centroids,
			cluster_assign)
    	cluster_assign_h[seed] = count_cluster_assign(cluster_assign)

    	ix_larg_clu = argmax(cluster_assign_h[seed])
    	larg_clu_sz = cluster_assign_h[seed][ix_larg_clu]

    	@printf("%s seed=%6d, heterogeneity=%.5f ", string(now()),
			seed, heterogeneity[seed])
        println("cluster_distribution=$(cluster_assign_h[seed]) / largest cluster: $(ix_larg_clu) size: $(larg_clu_sz)")

    	if larg_clu_sz > m_larg_clu_sz
       		m_ix_largest_clu, m_larg_clu_sz = ix_larg_clu, larg_clu_sz
		end
	end
  	@printf("took: %s\n", string(now()))
  	(cluster_assign_h, heterogeneity)
end

# ╔═╡ 175aa3f2-824c-11eb-01aa-f789353345f4
## with random init
begin
	cluster_assign_wri, heterogeneity_wri = nothing, nothing

	with_terminal() do
		global cluster_assign_wri, heterogeneity_wri =
			run_kmeans_with_multi_seeds(norm_ya_tf_idf)
	end
end

# ╔═╡ Cell order:
# ╟─f701bf4a-809e-11eb-2f1f-bb7dad75bfcc
# ╠═2a641766-809f-11eb-274a-8fdb234c0d08
# ╠═692a39c8-80a4-11eb-2077-e9f391787416
# ╠═5b3b148e-809f-11eb-3631-cf4682cbd959
# ╠═97d6cb18-809f-11eb-075a-7b3c5354ce4c
# ╠═97d2bd70-8239-11eb-2a8c-efea2584e6d6
# ╠═54ce09fa-8246-11eb-35fc-a9ba218d8c76
# ╟─6ada72c6-80a0-11eb-2346-c3554265a4de
# ╠═9ee529a8-80a0-11eb-23cc-b78ead6503e4
# ╠═16431f28-80a4-11eb-12a2-f3104b93f673
# ╟─e98e9a5e-80a3-11eb-1d38-db4ff0f9575b
# ╠═fbebf142-80a3-11eb-0913-4b43d543a681
# ╠═637c61ac-80a4-11eb-21be-f1661d34b568
# ╟─a30ee624-81e8-11eb-0c75-a5715165ff36
# ╠═4543b590-80a5-11eb-3215-b9bb37e1e766
# ╠═fd36c16e-80a6-11eb-3a80-5d7b8cb48b38
# ╠═4e94e7a2-80a7-11eb-0fb3-1955f645d0ca
# ╠═107cca2e-80ad-11eb-19b8-bf015a43b184
# ╠═69c34f00-80ac-11eb-0b03-0362d6f5e84e
# ╠═b0d25bf0-80ae-11eb-00de-d57e0d22a000
# ╟─07757e58-8247-11eb-299a-97cbb8ce1460
# ╠═99a6e9c4-8238-11eb-3acc-63a515a49b64
# ╠═3f2f94e4-80af-11eb-32d2-a5e101797f53
# ╠═fcfefb98-81e7-11eb-014a-b72b0477acab
# ╠═6c7e5bb4-823b-11eb-3833-254b4f5dafa7
# ╟─aaea9810-81e9-11eb-0452-83d979dca4c3
# ╠═aad0780c-81e9-11eb-2344-63996006b0cb
# ╟─446d44fa-81eb-11eb-0daa-83db97d8eec9
# ╟─f1494764-81eb-11eb-02e8-e9a44cbbd905
# ╠═444fc56a-81eb-11eb-27d2-1b249cae06d9
# ╠═9716ee8a-81fb-11eb-3279-c19f1a5736a2
# ╠═95d6fcdc-81ff-11eb-0d5b-5b161dba559f
# ╠═68962f52-8215-11eb-1cf2-1763b3ea2fcf
# ╠═aab3edea-81e9-11eb-2104-f18b12260e09
# ╟─aa994aa0-81e9-11eb-1fca-f1b45c0220f0
# ╠═8f71f718-8200-11eb-2dc1-4bc972d149f2
# ╠═20c45328-8237-11eb-1cf3-a9e5a4ee496f
# ╠═f812155c-8232-11eb-0113-5f5c9c8044e2
# ╠═19ac0ab6-8212-11eb-27ed-132248a13585
# ╟─756d72ec-8208-11eb-378d-2f9b2b288791
# ╠═9a0f5d4a-8208-11eb-11d0-9742f32e4705
# ╠═23ea71ba-822e-11eb-0086-493cd5a4ae95
# ╠═4f06ce52-822e-11eb-1c3f-5777aee2de37
# ╠═e58d0b00-8208-11eb-0414-019ec7c74eec
# ╠═f00143fa-820d-11eb-2983-43333f967308
# ╠═6ceb9cf8-820c-11eb-1652-99c659469ff9
# ╠═13440200-8250-11eb-0026-c91c8fdff383
# ╠═86d3d792-8250-11eb-36c0-33630e0c38ba
# ╠═509e66ca-8252-11eb-28da-dd665ec5575a
# ╠═58c012f2-8252-11eb-3b9a-6b736099ad9d
# ╠═c5137efa-820c-11eb-3773-057e25e6c242
# ╠═91866190-8213-11eb-00f2-036a1e31e428
# ╟─33b957ac-820e-11eb-17d6-a35b4a60fb2d
# ╠═338394f8-820e-11eb-0db4-1f86c79eacbc
# ╠═ab856686-8213-11eb-17e5-730b0a96ad12
# ╠═14fc4e02-8211-11eb-3f9d-d37995fb2522
# ╠═7d692dae-821c-11eb-2aef-97b108c8d8a5
# ╠═c6bb718c-821f-11eb-072c-bf89e4b16e6a
# ╠═e51692fe-8208-11eb-1fc9-1f8f3a9d5db6
# ╠═9330b2f2-8211-11eb-205e-7db725dfdf85
# ╟─92f9c26c-8211-11eb-1f71-818d8fd916fb
# ╠═92e15d44-8211-11eb-2977-8d3b271553e1
# ╠═92c7d752-8211-11eb-1546-cd510e3d8659
# ╟─92accbd8-8211-11eb-0703-0bb21268d6bc
# ╠═c43df0a0-8216-11eb-3be2-0f69c1faaba2
# ╟─c422dff2-8216-11eb-2f63-235d58c77e35
# ╠═c406c684-8216-11eb-3aba-65c9a55b7441
# ╠═c3ebe820-8216-11eb-06e1-538f0cda84b7
# ╠═31a82b96-8234-11eb-2795-bbadf977d8a3
# ╠═bf1fd5b6-8234-11eb-0fbc-91da7f6dbe28
# ╠═3633fa2a-8234-11eb-315d-d93094e52cab
# ╠═67ecd7b2-8234-11eb-03f5-41fd8adda1fd
# ╠═81b59776-8234-11eb-03d8-9736ea84e866
# ╠═a6d93074-8234-11eb-3cba-f1d9d8dbdfe3
# ╟─30eb64d4-824a-11eb-06c6-3393ca6a5780
# ╠═2c29b028-8253-11eb-2781-f3de2eb49ba3
# ╠═4c3ed812-824a-11eb-356a-6594a7b2a65e
# ╠═175aa3f2-824c-11eb-01aa-f789353345f4


text:
function perm_samp(data::PopData)
    popcounts = countmap(data.meta.population)
    pops, counts = keys(popcounts), values(popcounts)
    gdf = groupby(data.loci, :name)
    perm_idx = shuffle(Xoroshiro128Star(), keys(gdf))

    #return gdf

end


text:
using OrdinaryDiffEq, Plots
gr()

#Half-life of Carbon-14 is 5,730 years.
C₁ = 5.730

#Setup
u₀ = 1.0
tspan = (0.0, 1.0)

#Define the problem
radioactivedecay(u,p,t) = -C₁*u

#Pass to solver
prob = ODEProblem(radioactivedecay,u₀,tspan)
sol = solve(prob,Tsit5())

#Plot
plot(sol,linewidth=2,title ="Carbon-14 half-life", xaxis = "Time in thousands of years", yaxis = "Percentage left", label = "Numerical Solution")
plot!(sol.t, t->exp(-C₁*t),lw=3,ls=:dash,label="Analytical Solution")

max_stars_repo_name: sergev/vak-opensource

max_issues_repo_name: sergev/vak-opensource

max_forks_repo_name: sergev/vak-opensource


text:
function get_trends_place(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/trends/place.json", options))
end

function get_trends_available(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/trends/available.json", options))
end

function get_trends_closest(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/trends/closest.json", options))
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
using BenchmarkTools
using FLOWFarm; const ff=FLOWFarm
using DelimitedFiles
using ForwardDiff
using Profile
# using ProfileView

function benchmark_point_in_polygon()

    # vertices = [0.0 0.0; 0.0 10.0; 10.0 10.0; 10.0 0.0] # square
    # vertices = [0.0 0.0; 0.0 10.0; 10.0 10.0; 10.0 6.0; 5.0 6.0; 5.0 4.0; 10.0 4.0; 10.0 0.0] # square with cut out
    # vertices = [0.0 0.0; 0.0 10.0; 4.0 4.0; 10.0 2.0] # almost a triangle
    vertices = ff.star_boundary(5, 5.0, 10.0) .+ 5.0
    # println(vertices)
    normals = ff.boundary_normals_calculator(vertices)
    step = 5E-0
    pointx = -5:step:15.0
    pointy = -5:step:15.0

    # set up function for getting AD derivatives
    pointinpolygon_diff(x) = ff.pointinpolygon(x, vertices, normals, return_distance=true)
    # @benchmark($pointinpolygon_diff([$pointx[1], $pointy[1]]))
    points = ([pointx[i], pointy[j]] for i = 1:length(pointx) for j = 1:length(pointx))
    # @benchmark [$pointinpolygon_diff(point) for point in $points]
    @profile [pointinpolygon_diff(point) for point in points]
    # @benchmark [ForwardDiff.gradient($pointinpolygon_diff, point) for point in $points]
    # distance = zeros((length(pointx), length(pointy)))
    # xderiv = zeros((length(pointx), length(pointy)))
    # yderiv = zeros((length(pointx), length(pointy)))

    # for i = 1:length(pointx)
    #     for j = 1:length(pointy)
    #         point = [pointx[i], pointy[j]]
    #         distance[j,i] = ff.pointinpolygon(point, vertices, normals, return_distance=true)
    #         deriv = ForwardDiff.gradient(pointinpolygon_diff, point)
    #         xderiv[j,i] = deriv[1]
    #         yderiv[j,i] = deriv[2]
    #     end
    # end

    # fig, ax = plt.subplots(2,2)
    # ff.plotboundary!(ax[1,1], vertices)
    # maxmin = maximum(abs.(distance))
    # cplot1 = ax[1,1].contourf(pointx, pointy, distance, cmap="PuOr", vmin=-maxmin, vmax=maxmin, levels=200)
    # cbar1 = plt.colorbar(cplot1)

    # ff.plotboundary!(ax[1,2], vertices)
    # maxmin = maximum(abs.(xderiv))
    # cplot2 = ax[1,2].contourf(pointx, pointy, xderiv, cmap="PuOr", vmin=-maxmin, vmax=maxmin, levels=200)
    # cbar2 = plt.colorbar(cplot2)

    # ff.plotboundary!(ax[2,1], vertices)
    # maxmin = maximum(abs.(yderiv))
    # cplot3 = ax[2,1].contourf(pointx, pointy, yderiv, cmap="PuOr", vmin=-maxmin, vmax=maxmin, levels=200)
    # cbar3 = plt.colorbar(cplot3)

end


text:
include("testsystem1.jl")

U = @inferred trace(LinearEvolution(Φ, Gaussian(b, Q)), 0 => Gaussian(x0, P0), endtime(10000))

X = @inferred trace(DynamicIterators.Sampled(LinearEvolution(Φ, Gaussian(b, Q))), 0 => x0, endtime(10000))

@show norm(mean(values(X)) - b) < eps()
@show cov(values(X))
@show cov(values(U)[end])
@test norm(cov(values(X)) - cov(values(U)[end])) < 0.2

using DynamicIterators: dyniterate

Random.seed!(11)
X = @inferred trace(DynamicIterators.Sampled(LinearEvolution(Φ, Gaussian(b, Q))), 0 => x0, endtime(10))
Y = Any[0 => x0]
Random.seed!(11)
Y = collecttrace(LinearEvolution(Φ, Gaussian(b, Q)), Trace(Sample(0 => x0), Y, endtime(10)))
@test X == trajectory(Y)


text:
using JuLIP, StaticArrays, NearestNeighbors, LinearAlgebra
using JuLIP: AbstractAtoms, AbstractCalculator

import Base: *

export strains

const Vec2{T} = SArray{Tuple{2},T,1,2}
const Mat2{T} = SArray{Tuple{2,2},T,2,4}
const Ten32{T} = SArray{Tuple{2,2,2},T,3,8}
const Ten42{T} = SArray{NTuple{4,2},T,4,16}
const MVec2{T} = MVector{2,T}
const MMat2{T} = MMatrix{2,2,T}
const MTen32{T} = MArray{Tuple{2,2,2},T,3,8}
const MTen42{T} = MArray{NTuple{4,2},T,4,16}

const Vec3{T} = SArray{Tuple{3},T,1,3}
const Mat3{T} = SArray{Tuple{3,3},T,2,9}
const Ten33{T} = SArray{Tuple{3,3,3},T,3,27}
const Ten43{T} = SArray{NTuple{4,3},T,4,81}
const MVec3{T} = MVector{3, T}
const MMat3{T} = MMatrix{3,3,T}
const MTen33{T} = MArray{Tuple{3,3,3},T,3,27}
const MTen43{T} = MArray{NTuple{4,3},T,4,81}

const _EE = Mat3{Float64}([1 0 0; 0 1 0; 0 0 1])
ee(i::Integer) = _EE[:,i]

"""
`strains(U, at; rcut = cutoff(calculator(at)))`

returns maximum strains :  `maximum( du/dr over all neighbours )` at each atom
"""
strains(U, at; rcut = cutoff(calculator(at))) =
   [ maximum(norm(u - U[i]) / s for (u, s) in zip(U[j], r))
                           for (i, j, r, _1) in sites(at, rcut) ]

struct ForceConstantMatrix1{T <: AbstractFloat}
   R::Vector{Vec3{T}}
   H::Vector{Mat3{T}}
end

function *(fcm::ForceConstantMatrix1{T}, t::Tuple{<:Atoms, Vector{<: Vec3}}) where {T}
   at, U = t
   V = zeros(Vec3{T}, length(at))
   rcut = maximum(norm.(fcm.R)) * 1.01
   # WARNING: this is an extremely naive implementation which needs to be fixed
   #          as soon as possible
   for (i, j, R) in sites(at, rcut)
      for n = 1:length(R)
         for m = 1:length(fcm.R)
            if norm(R[n] - fcm.R[m]) < 1e-5
               V[i] += fcm.H[m] * (U[j[n]] - U[i])
               break
            end
         end
      end
   end
   return V
end

ForceConstantMatrix1(calc::AbstractCalculator, at::AbstractAtoms; kwargs...) =
      ForceConstantMatrix1(force_constants(calc, at; kwargs...)...)

function force_constants(calc::AbstractCalculator, at::Atoms{T};
                         h = 1e-5, rcut = 2.01*cutoff(calc) + 1
                         ) where {T}
   @assert length(at) == 1 # assume no basis
   cl = cluster(at, rcut)
   set_calculator!(cl, calc)
   x = cl[1]::Vec3{T}  # x ≡ X[1]
   ∂f_∂xi = []
   for i = 1:3
      cl[1] = x + h * ee(i)
      fip = forces(cl)
      cl[1] = x - h * ee(i)
      fim = forces(cl)
      cl[1] = x
      push!(∂f_∂xi, (fip - fim) / (2*h))
   end
   # convert to dynamical matrix entries
   H = [ - [∂f_∂xi[1][n] ∂f_∂xi[2][n] ∂f_∂xi[3][n]]'  for n = 1:length(cl) ]
   # extract the non-zero entries
   Inz = setdiff(findall( norm.(H) .> 1e-8 ), [1])
   # . . . and return
   R = Vec3{T}[y - x for y in positions(cl)[Inz]]
   return R, H[Inz]
end

"""
`findin(Xsm, Xlge)`

Assuming that `Xsm ⊂ Xlge`, this function
returns `Ism, Ilge`, both `Vector{Int}` such that
* Xsm[i] == Xlge[Ism[i]]
* Xlge[i] == Xsm[Ilge[i]]  if Xlge[i] ∈ Xsm; otherwise Ilge[i] == 0
"""
function findin2(Xsm::Vector{SVec{T}}, Xlge::Vector{SVec{T}}) where T <: AbstractFloat
   # find the nearest neighbours of Xsm points in Xlge
   tree = KDTree(Xlge)
   # construct the Xsm -> Xlge mapping
   Ism = zeros(Int, length(Xsm))
   Ilge = zeros(Int, length(Xlge))
   for (n, x) in enumerate(Xsm)
      i = inrange(tree, Xsm, sqrt(eps(T)))
      if isempty(i)
         Ism[n] = 0         # - Ism[i] == 0   if  Xsm[i] ∉ Xlge
      elseif length(i) > 1
         error("`inrange` found two neighbours")
      else
         Ism[n] = i[1]      # - Ism[i] == j   if  Xsm[i] == Xlge[j]
         Ilge[i[1]] = n     # - Ilge[j] == i  if  Xsm[i] == Xlge[j]
      end
   end
   # - if Xlge[j] ∉ Xsm then Ilge[j] == 0
   return Ism, Ilge
end

# """
# `relative_displacement(Xsm, Xrefsm, Xlge, Xreflge)`
#
# Assumes that Xrefsm ⊂ Xreflge;
# compute a relative displacement `U` such that
# * length(U) == length(Xlge)
# * Xsm[i] = Xlge[j] + U[j]  when  Xrefsm[i] == Xreflge[i]
# * U[j] = 0 otherwise
# """
# function relative_displacement(Xsm, Xrefsm, Xlge, Xreflge)
#    @assert length(Xsm) == length(Xrefsm)
#    @assert length(Xlge) == length(Xreflge)
#    Ism, _ = findin(Xrefsm, Xreflge)
#    @assert maximum(norm.(Xrefsm - Xreflge[Ism])) < 1e-7
#    U = Xlge - Xreflge
#    U[Ism] = Xsm - Xlge[Ism]  # = (Xsm - Xrefsm) - (Xlge - Xreflge)
#    return U
# end      # TODO: U should not be zero outside Xsm, but U = Xlge - Xlgeref


text:
################################################################################
# Spring Impulses
################################################################################

function spring_force(relative::Symbol, joint::Rotational{T},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion;
        rotate::Bool=true,
        unitary::Bool=false) where T

    spring = unitary ? 1.0 : joint.spring
    distance = joint.spring_offset .- minimal_coordinates(joint, xa, qa, xb, qb)
    force = -spring * zerodimstaticadjoint(nullspace_mask(joint)) * distance # force in offset frame

    if relative == :parent
        rotate ? (output = vector_rotate(force, joint.axis_offset)) : (output = force) # rotate into a frame
    elseif relative == :child
        rotate ? (output = vector_rotate(-force, inv(qb) * qa  * joint.axis_offset)) : (output = -force) # rotate back to b frame
    end

    return [szeros(T, 3); output]
end

function spring_impulses(relative::Symbol, joint::Rotational, pbody::Node, cbody::Node, timestep;
    unitary::Bool=false)
    spring_impulses(relative, joint,
        current_configuration(pbody.state)...,
        current_configuration(cbody.state)...,
        timestep, unitary=unitary)
end

function spring_impulses(relative::Symbol, joint::Rotational,
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion,
    timestep;
    unitary::Bool=false)
    timestep * spring_force(relative, joint, xa, qa, xb, qb; unitary=unitary)
end

spring_impulses(relative::Symbol, joint::Rotational{T,3}, pbody::Node, cbody::Node, timestep; unitary::Bool=false) where T = szeros(T, 6)

################################################################################
# Spring Jacobians
################################################################################

function spring_jacobian_configuration(relative::Symbol, jacobian::Symbol, joint::Rotational{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion,
    timestep::T;
    rotate::Bool=true,
    unitary::Bool=false,
    attjac=true) where T

    Z = attjac ? szeros(T, 3, 6) : szeros(T, 3, 7)
    force = spring_force(relative, joint, xa, qa, xb, qb, rotate=false)
    spring = unitary ? 1.0 : joint.spring

    if relative == :parent
        J = timestep * rotation_matrix(joint.axis_offset) * spring * zerodimstaticadjoint(nullspace_mask(joint)) * minimal_coordinates_jacobian_configuration(jacobian, joint, xa, qa, xb, qb, attjac=attjac)
    elseif relative == :child
        X = szeros(T, 3, 3)

        if rotate
            J1 = timestep * rotation_matrix(inv(qb) * qa * joint.axis_offset) * -spring * zerodimstaticadjoint(nullspace_mask(joint)) * minimal_coordinates_jacobian_configuration(jacobian, joint, xa, qa, xb, qb, attjac=attjac)

            if jacobian == :parent
                Q2 = timestep * ∂vector_rotate∂q(force[SVector{3,Int}(4,5,6)], inv(qb) * qa * joint.axis_offset) * Rmat(joint.axis_offset) * Lmat(inv(qb))
                attjac && (Q2 *= LVᵀmat(qa))
                J2 = [X Q2]
            elseif jacobian == :child
                Q2 = timestep * ∂vector_rotate∂q(force[SVector{3,Int}(4,5,6)], inv(qb) * qa * joint.axis_offset) * Rmat(qa * joint.axis_offset) * Tmat()
                attjac && (Q2 *= LVᵀmat(qb))
                J2 = [X Q2]
            end
            J = J1 + J2
        else
            J = timestep * -spring * zerodimstaticadjoint(nullspace_mask(joint)) * minimal_coordinates_jacobian_configuration(jacobian, joint, xa, qa, xb, qb, attjac=attjac)
        end
    end

    return [Z; J]
end

function spring_jacobian_configuration(relative::Symbol, jacobian::Symbol,
    joint::Rotational, pbody::Node, cbody::Node,
    timestep::T;
    attjac::Bool = true) where T

    spring_jacobian_configuration(relative, jacobian, joint,
        current_configuration(pbody.state)...,
        current_configuration(cbody.state)...,
        timestep; attjac=attjac)
end

spring_jacobian_configuration(relative::Symbol, jacobian::Symbol, joint::Rotational{T,3}, pbody::Node, cbody::Node, timestep::T; attjac::Bool = true) where T = attjac ? szeros(T, 6, 6) : szeros(T, 6, 7)
spring_jacobian_velocity(relative::Symbol, jacobian::Symbol, joint::Rotational, pbody::Node, cbody::Node, timestep::T) where T = szeros(T, 6, 6)


text:
"""
  function fasta2matrix(filename::AbstractString; max_gap_fraction = 1)

  Reads an aligned fasta file and returns a matrix with amino acids converted
  to numbers following the convention A --> 1, C --> 2, ... '-' --> 21.

    "filename": full path of the fasta file.
    "max_gap_fraction": to be removed.

"""

function fasta2matrix(filename::AbstractString; max_gap_fraction = 1)
    f = FastaReader(filename)

    max_gap_fraction = Float64(max_gap_fraction)

    # pass 1

    seqs = Int[]
    inds = Int[]
    fseqlen = 0

    for (name, seq) in f
        ngaps = 0
        if f.num_parsed == 1
            ls = length(seq)
            resize!(inds, ls)
            for i = 1:ls
                c = seq[i]
                if c != '.' && c == uppercase(c)
                    fseqlen += 1
                    inds[fseqlen] = i
                    c == '-' && (ngaps += 1)
                end
            end
        else
            ls = length(seq)
            ls == length(inds) || error("inputs are not aligned")
            tstfseqlen = 0
            for i = 1:ls
                c = seq[i]
                if c != '.' && c == uppercase(c)
                    tstfseqlen += 1
                    inds[tstfseqlen] == i || error("inconsistent inputs")
                    c == '-' && (ngaps += 1)
                end
            end
            tstfseqlen == fseqlen || error("inconsistent inputs")
        end
        ngaps / fseqlen <= max_gap_fraction && push!(seqs, f.num_parsed)
    end

    length(seqs) > 0 || error("Out of $(f.num_parsed) sequences, none passed the filter (max_gap_fraction=$max_gap_fraction)")

    # pass 2
    Z = Array{Int16}(undef, fseqlen, length(seqs))
    seqid = 1
    for (name, seq) in f
        seqs[end] < f.num_parsed && break
        seqs[seqid] == f.num_parsed || continue
        for i = 1:fseqlen
            c = seq[inds[i]]
            Z[i, seqid] = letter2num(c)
        end
        seqid += 1
    end
    @assert seqid == length(seqs) + 1

    close(f)

    return Z'
end

let alphabet = [ 1, 21, 2, 3, 4, 5, 6, 7, 8, 21, 9, 10, 11, 12, 21, 13, 14, 15, 16, 17, 21, 18, 19, 21, 20]
               # A, B,  C, D, E, F, G, H, I, J,  K, L,  M,  N,  O,  P,  Q,  R,  S,  T,  U,  V,  W,  X,  Y
    global letter2num
    function letter2num(c::Union{Char,UInt8})
        i = UInt8(c) - 0x40
        1 <= i <= 25 && return alphabet[i]
        return 21
     end
end

"""
	read_par_BM(path::AbstractString; q::Integer = 21)

	Reads the parameters of a Potts model in format

	J i j a b
	...
	h i a

	and returns them in tensor format.
	J[a, b, i, j] is such that  1 == "A", ... "21" == "-",
	same for h[a, i].

    "path": full path of the parameter file.
    "q": number of symbols allowed (21 for protein sequences, gaps included).

"""

function read_par_BM(path::AbstractString, q::Integer = 21)
    params = readdlm(path,' ', use_mmap = true)[:, 2:6]
    l_file = size(params, 1)
    N = Integer(((q - 2) + sqrt( (q-2)^2 + 8*l_file))/(2*q))
	J = Array{Float64}(undef, q, q, N, N)
	h = Array{Float64}(undef, q, N)
	n_J = Int(q*q*N*(N-1)/2)
	n_h = q*N

	for k in 1:n_J
		i, j, a, b, par_j = params[k, :]
		i += 1
		j += 1
		a == 0 && (a = 21)
		b == 0 && (b = 21)
		J[a, b, i, j] = par_j
	end

	for l in (n_J + 1): n_h + n_J
		i, a, par_h = params[l, :]
		i += 1
		a == 0 && (a = 21)
		h[a, i] = par_h
	end

	return h, J
end

function read_par_BM(path::GZipStream, q::Integer = 21)
    params = readdlm(path,' ', use_mmap = true)[:, 2:6]
    l_file = size(params, 1)
    N = Integer(((q - 2) + sqrt( (q-2)^2 + 8*l_file))/(2*q))
	J = Array{Float64}(undef, q, q, N, N)
	h = Array{Float64}(undef, q, N)
	n_J = Int(q*q*N*(N-1)/2)
	n_h = q*N

	for k in 1:n_J
		i, j, a, b, par_j = params[k, :]
		i += 1
		j += 1
		a == 0 && (a = 21)
		b == 0 && (b = 21)
		J[a, b, i, j] = par_j
	end

	for l in (n_J + 1): n_h + n_J
		i, a, par_h = params[l, :]
		i += 1
		a == 0 && (a = 21)
		h[a, i] = par_h
	end

	return h, J
end

function extract_params(path_par::AbstractString; q::Integer = 21)
	!isfile(path_par) && error("Error: the file \"$(path_params)\" does not exist. Please check the spelling or the folder path.")
	file = GZip.open(path_par)
	h, J = read_par_BM(file, q)
	h = set_max_field_to_0(h)
	J = symmetrize_J(J)
	return h, J
end

"""
    num2letter(i::Integer)

    Takes as input an integer (representing an amino acid)
    and returns its conventional character representation.
    In this case with the convention "1 2 .. 21" == "A C .. -"
"""

let alphabet = ["A", "C", "D", "E", "F", "G", "H", "I",  "K", "L",  "M",  "N", "P",  "Q",  "R",
"S",  "T", "V",  "W",  "Y"]
    global num2letter
    function num2letter(i :: Integer)
        1 <= i <= 20 && return alphabet[i]
        return "-"
    end
end

"""
    vec2string(v)

    Takes as input a vector of integers (representing an amino acids)
    and returns a list of characters, the corresponding amino acids.
    In this case with the convention "1 2 .. 21" == "A C .. -"
"""

function vec2string(v)
    s = ""
    for i in v
        s = s*num2letter(i)
    end
    return s
end

"""
    set_max_field_to_0(h)

    Takes in input a 2-dimensional tensor, a matrix,
    containing the values of the fiels for every amino acid a
    and position i.
    Returns the same matrix, by removing the biggerst field value
    for every position i, so that the biggest field results 0.
    Help in computing exponentials without going to overflow.

"""

function set_max_field_to_0(h_old::Array{Float64, 2})
   h_new = deepcopy(h_old)
   q, N = size(h_old)
   for i in 1:N
        hmax = maximum([h_old[a, i] for a in 1:21])
        h_new[:, i] .-= hmax
   end
   return h_new
end

"""
    symmetrize_couplings(J::Array{} )

    Takes as input a 4-dimensional tensor,
    containing the values of the couplings for every couple of
    amino acids a, b and positions i, j.
    Returns the same tensor, by symmetrizing it,
    since J[a, b, i, j] = 0 if i < j, or viceversa, depends.

"""

function symmetrize_J(J_old::Array{Float64, 4})
   J_s = deepcopy(J_old)
   q,q, N, N = size(J_old)
   for i in 1:N
        for j in i+1:N
            for a in 1:q
                for b in 1:q
                      J_s[b, a, j, i]  = J_old[a, b, i, j]
                end
            end
        end
    end
    return J_s
end

max_stars_repo_name: matteobisardi/SeqEvol

max_issues_repo_name: matteobisardi/SeqEvol

max_forks_repo_name: matteobisardi/SeqEvol


text:
# using Dates

@testset "Frp.jl" begin
    filename = "data\\test1.frp"
    dc = DDR2import.Frp.read(filename)
    @test dc["ALBANIA"].freeroutepoints[1].type == "E"
    @test dc["ALBANIA"].freeroutepoints[2].name == "NIKRO"
    # @test dc["ALBANIA"].freeroutepoints[3].point.lat ≈ 41.775277777777774 atol = 0.001
    # @test dc["ALBANIA"].freeroutepoints[4].point.lon ≈ 20.472777777777775 atol = 0.001
    @test dc["ALBANIA"].freeroutepoints[3].point.ϕ ≈ 41.775277777777774 atol = 0.001
    @test dc["ALBANIA"].freeroutepoints[4].point.λ ≈ 20.472777777777775 atol = 0.001
    @test dc["ALBANIA"].freerouteairports[1].type == "A"
    @test dc["ALBANIA"].freerouteairports[2].name == "DIRES"
    @test dc["ALBANIA"].freerouteairports[3].airports[1] == "LATI"
    @test isempty(dc["BIRD_DOMESTIC"].freerouteairports[3].airports) == true
    @test dc["EI_FRA"].freerouteairports[18].airports[2] == "EINN"
    #LAST
    @test dc["FRAC"].freerouteairports[1].type == "A"

    filename2 = "data\\test2.frp"
    dc2 = DDR2import.Frp.read(filename2)
    @test dc2["NAT4"].freeroutepoints[1].name == "*1847"
    # @test dc2["NAT4"].freeroutepoints[1].point.lat == 18.0
    # @test dc2["NAT4"].freeroutepoints[1].point.lon == -47.5
    @test dc2["NAT4"].freeroutepoints[1].point.ϕ == 18.0
    @test dc2["NAT4"].freeroutepoints[1].point.λ == -47.5
    # @test dc2["NEFRA"].freeroutepoints[1129].point.lat ≈ 58.1456667 atol = 0.00001
    # @test dc2["SEEN_FRA_NORTH"].freeroutepoints[184].point.lon ≈ 18.35583333 atol = 0.00001
    @test dc2["NEFRA"].freeroutepoints[1129].point.ϕ ≈ 58.1456667 atol = 0.00001
    @test dc2["SEEN_FRA_NORTH"].freeroutepoints[184].point.λ ≈ 18.35583333 atol = 0.00001

end


text:
# This file is a part of ParallelProcessingTools.jl, licensed under the MIT License (MIT).

export ThreadSafe
abstract type ThreadSafe{T} end

export ThreadSafeReentrantLock

@static if VERSION >= v"1.2-DEV.28"
    const ThreadSafeReentrantLock = ReentrantLock
else
    struct ThreadSafeReentrantLock
        thread_lock::RecursiveSpinLock
        task_lock::ReentrantLock

        ThreadSafeReentrantLock() = new(RecursiveSpinLock(), ReentrantLock())
    end

    function Base.lock(l::ThreadSafeReentrantLock)
        # @debug "LOCKING $l"
        lock(l.thread_lock)
        try
            lock(l.task_lock)
        catch err
            unlock(l.thread_lock)
            rethrow()
        end
    end

    function Base.unlock(l::ThreadSafeReentrantLock)
        # @debug "UNLOCKING $l"
        try
            unlock(l.task_lock)
        finally
            unlock(l.thread_lock)
        end
    end
end

export LockableValue

struct LockableValue{T} <: ThreadSafe{T}
    x::T
    l::ThreadSafeReentrantLock

    (::Type{LockableValue{T}})(x::T) where {T} = new{T}(x, ThreadSafeReentrantLock())
end

LockableValue(x::T) where {T} = LockableValue{T}(x)

function Base.broadcast(f, lx::LockableValue)
    lock(lx.l) do
        f(lx.x)
    end
end

Base.map(f, lx::LockableValue) = broadcast(f, lx)

abstract type ThreadSafeIO <: IO end

export LockableIO

struct LockableIO{T<:IO}
    lx::LockableValue{T}

    (::Type{LockableIO{T}})(x::T) where {T<:IO} = new{T}(LockableValue(x))
end

LockableIO(x::T) where {T<:IO} = LockableIO{T}(x)

@inline Base.broadcast(f, lio::LockableIO) = broadcast(f, lio.lx)
@inline Base.map(f, lio::LockableIO) = broadcast(f, lio)

@inline Base.read(lio::LockableIO, args...; kwargs...) = map(lio) do io
    read(io, args...; kwargs...)
end

@inline Base.read!(lio::LockableIO, args...; kwargs...) = map(lio) do io
    read(io, args...; kwargs...)
end

@inline Base.write(lio::LockableIO, args...; kwargs...) = map(lio) do io
    write(io, args...; kwargs...)
end

const _critical_section_lock = ThreadSafeReentrantLock()

"""
    @critical expr

Mark code in `expr` as a critical section. Code in critical sections will
never be executed in parallel (via multithreading) to any other critical
section.

`@critical` is very useful to mark non-threadsafe code.

Example:

```julia
@onthreads allthreads() begin
    @critical @info Base.Threads.threadid()
end

Without `@critical`, the above will typically crash Julia.
```
"""
macro critical(expr)
    quote
        try
            lock(_critical_section_lock)
            $(esc(expr))
        finally
            unlock(_critical_section_lock)
        end
    end
end
export @critical


text:
module Effects

export @effect, effect

include("cps.jl")
include("handler.jl")

include("list.jl")

end # module


text:
using Observables

abstract type AbstractWidget{T} <: Observables.AbstractObservable{T} end

Observables.observe(x::AbstractWidget) = x.value

struct Button{T} <: AbstractWidget{Bool}
    content::Observable{T}
    value::Observable{Bool}
    attributes::Dict{Symbol, Any}
end

function Button(content; kw...)
    return Button(
        Observable(content), Observable(false), Dict{Symbol, Any}(kw)
    )
end

function jsrender(button::Button)
    return DOM.input(
        type = "button",
        value = button.content,
        onclick = js"update_obs($(button.value), true);";
        button.attributes...
    )
end

struct TextField <: AbstractWidget{String}
    value::Observable{String}
    attributes::Dict{Symbol, Any}
end

function TextField(value::String; kw...)
    TextField(Observable(value), Dict{Symbol, Any}(kw))
end

function jsrender(tf::TextField)
    return DOM.input(
        type = "textfield",
        value = tf.value,
        onchange = js"update_obs($(tf.value),  this.value);";
        tf.attributes...
    )
end

struct Slider{T <: AbstractRange, ET} <: AbstractWidget{T}
    range::Observable{T}
    value::Observable{ET}
    attributes::Dict{Symbol, Any}
end

to_node(x) = Observable(x)
to_node(x::Observable) = x

function Slider(range::T, value = first(range); kw...) where T <: AbstractRange
    Slider{T, eltype(range)}(
        to_node(range),
        to_node(value),
        Dict{Symbol, Any}(kw)
    )
end

function jsrender(slider::Slider)
    return DOM.input(
        type = "range",
        min = map(first, slider.range),
        max = map(last, slider.range),
        value = slider.value,
        step = map(step, slider.range),
        oninput = js"update_obs($(slider.value), parseFloat(value))";
        slider.attributes...
    )
end

@enum Orientation vertical horizontal

struct RangeSlider{T <: AbstractRange, ET <: AbstractArray} <: AbstractWidget{T}
    range::Observable{T}
    value::Observable{ET}
    attributes::Dict{Symbol, Any}
    connect::Observable{Bool}
    orientation::Observable{Orientation}
end

function RangeSlider(range::T; value = [first(range)], kw...) where T <: AbstractRange
    RangeSlider{T, typeof(value)}(
        Observable(range),
        Observable(value),
        Dict{Symbol, Any}(kw),
        Observable(true),
        Observable(horizontal)
    )
end

const noUiSlider = Dependency(
    :noUiSlider,
    [
        "https://cdn.jsdelivr.net/gh/leongersen/noUiSlider/distribute/nouislider.min.js",
        "https://cdn.jsdelivr.net/gh/leongersen/noUiSlider/distribute/nouislider.min.css"
    ]
)

function jsrender(session::Session, slider::RangeSlider)
    rangediv = DOM.div()
    onload(session, rangediv, js"""
        function create_slider(range){
            $(noUiSlider).create(range, {

                range: {
                    'min': $(fill(slider.range[][1], length(slider.value[]))),
                    'max': $(fill(slider.range[][end], length(slider.value[])))
                },

                step: $(step(slider.range[])),

                // Handles start at ...
                start: $(slider.value[]),
                // Display colored bars between handles
                connect: $(slider.connect[]),
                tooltips: true,

                // Put '0' at the bottom of the slider
                direction: 'ltr',
                orientation: $(string(slider.orientation[])),
            });

            range.noUiSlider.on('update', function (values, handle, unencoded, tap, positions){
                update_obs($(slider.value), [parseFloat(values[0]), parseFloat(values[1])]);
            });
        }
    """)
    return rangediv
end