text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiStorageV1VolumeNodeResources <: SwaggerModel
    count::Any # spec type: Union{ Nothing, Int32 } # spec name: count

    function IoK8sApiStorageV1VolumeNodeResources(;count=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1VolumeNodeResources, Symbol("count"), count)
        setfield!(o, Symbol("count"), count)
        o
    end
end # type IoK8sApiStorageV1VolumeNodeResources

const _property_map_IoK8sApiStorageV1VolumeNodeResources = Dict{Symbol,Symbol}(Symbol("count")=>Symbol("count"))
const _property_types_IoK8sApiStorageV1VolumeNodeResources = Dict{Symbol,String}(Symbol("count")=>"Int32")
Base.propertynames(::Type{ IoK8sApiStorageV1VolumeNodeResources }) = collect(keys(_property_map_IoK8sApiStorageV1VolumeNodeResources))
Swagger.property_type(::Type{ IoK8sApiStorageV1VolumeNodeResources }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiStorageV1VolumeNodeResources[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1VolumeNodeResources }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1VolumeNodeResources[property_name]

function check_required(o::IoK8sApiStorageV1VolumeNodeResources)
    true
end

function validate_property(::Type{ IoK8sApiStorageV1VolumeNodeResources }, name::Symbol, val)
end


text:
# ---
# title: 466. Count The Repetitions
# id: problem466
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/count-the-repetitions/description/>
# hidden: true
# ---
#
# Define `S = [s,n]` as the string S which consists of n connected strings s.
# For example, `["abc", 3]` ="abcabcabc".
#
# On the other hand, we define that string s1 can be obtained from string s2 if
# we can remove some characters from s2 such that it becomes s1. For example,
# “abc” can be obtained from “abdbec” based on our definition, but it can not be
# obtained from “acbbe”.
#
# You are given two non-empty strings s1 and s2 (each at most 100 characters
# long) and two integers 0 ≤ n1 ≤ 106 and 1 ≤ n2 ≤ 106. Now consider the strings
# S1 and S2, where `S1=[s1,n1]` and `S2=[s2,n2]`. Find the maximum integer M
# such that `[S2,M]` can be obtained from `S1`.
#
# **Example:**
#
#
#
#     Input:
#     s1="acb", n1=4
#     s2="ab", n2=2
#
#     Return:
#     2
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
# Autogenerated wrapper script for PORTA_jll for x86_64-unknown-freebsd
export libporta, valid, xporta

JLLWrappers.@generate_wrapper_header("PORTA")
JLLWrappers.@declare_library_product(libporta, "libporta.so")
JLLWrappers.@declare_executable_product(valid)
JLLWrappers.@declare_executable_product(xporta)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libporta,
        "lib/libporta.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        valid,
        "bin/valid",
    )

    JLLWrappers.@init_executable_product(
        xporta,
        "bin/xporta",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using LinearMaps
using TensorOperations

"""
TODO:
  * include IterativeSolvers option for als
  * refactor als_gen_eigsolv

Implementation based on the presentation in
Holtz, Sebastian, Thorsten Rohwedder, and Reinhold Schneider. "The alternating linear scheme for tensor optimization in the tensor train format." SIAM Journal on Scientific Computing 34.2 (2012): A683-A713.
"""

function init_H(x_tt::ttvector{T},A_tto::ttoperator{T}) where T<:Number
	d = length(x_tt.ttv_dims)
	H = Array{Array{T}}(undef, d)
	H[d] = ones(T,1,1,1)
	for i = d : -1 : 2
		H[i-1] = zeros(T,A_tto.tto_rks[i],x_tt.ttv_rks[i],x_tt.ttv_rks[i])
		x_vec = x_tt.ttv_vec[i]
		A_vec = A_tto.tto_vec[i]
		update_H!(x_vec,A_vec,H[i],H[i-1])
	end
	return H
end

function update_H!(x_vec::Array{T,3},A_vec::Array{T,4},Hi::Array{T,3},Him::Array{T,3}) where T<:Number
	@tensoropt((ϕ,χ), Him[a,α,β] = conj.(x_vec)[j,α,ϕ]*Hi[z,ϕ,χ]*x_vec[k,β,χ]*A_vec[j,k,a,z]) #size (rim, rim, rAim)
	nothing
end

function init_Hb(x_tt::ttvector{T},b_tt::ttvector{T}) where T<:Number
	d = length(x_tt.ttv_dims)
	H_b = Array{Array{T}}(undef, d)
	H_b[d] = ones(T,1,1)
	for i = d : -1 : 2
		H_b[i-1] = zeros(T,x_tt.ttv_rks[i],b_tt.ttv_rks[i])
		b_vec = b_tt.ttv_vec[i]
		x_vec = x_tt.ttv_vec[i]
		update_Hb!(x_vec,b_vec,H_b[i],H_b[i-1]) #size(rbim, rim)
	end
	return H_b
end

function update_Hb!(x_vec::Array{T,3},b_vec::Array{T,3},H_bi::Array{T,2},H_bim::Array{T,2}) where T<:Number
	@tensoropt((ϕ,χ), H_bim[α,β] = H_bi[ϕ,χ]*b_vec[i,β,χ]*conj.(x_vec)[i,α,ϕ])
	nothing
end

function update_G!(x_vec::Array{T,3},A_vec::Array{T,4},Gi::AbstractArray{T,5},Gip::AbstractArray{T,5}) where T<:Number
	@tensor Gip[j,α,k,β,J] = (conj.(x_vec)[l,ϕ,α]*(Gi[l,ϕ,m,χ,L]*x_vec[m,χ,β]))*A_vec[j,k,L,J]
	nothing
end

function update_Gb!(x_vec::Array{T,3},b_vec::Array{T,3},G_bi::AbstractArray{T,3},G_bip::AbstractArray{T,3}) where T<:Number
	@tensoropt((ϕ,χ), G_bip[i,α,β] = b_vec[i,ϕ,β]*G_bi[j,χ,ϕ]*conj.(x_vec)[j,χ,α])
	nothing
end

#full assemble of matrix K
function K_full(Gi::Array{T,5},Hi::Array{T,3},K_dims::NTuple{3,Int}) where T<:Number
	K = zeros(T,prod(K_dims),prod(K_dims))
	Krshp = reshape(K,(K_dims...,K_dims...))
	@tensor Krshp[a,b,c,d,e,f] = Gi[a,b,d,e,z]*Hi[z,c,f] #size (ni,rim,ri,ni,rim,ri)
	return K
end

function Ksolve(Gi::Array{T,5},G_bi::Array{T,3},Hi::Array{T,3},H_bi::Array{T,2}) where T<:Number
	K_dims = (size(Gi,1),size(Gi,2),size(Hi,2))
	K = K_full(Gi,Hi,K_dims)
	@tensor Pb[i,α1,α2] := G_bi[i,α1,β]*H_bi[α2,β] #size (ni,rim,ri)
	return reshape(K\Pb[:],K_dims)
end

function K_eigmin(Gi::Array{T,5},Hi::Array{T,3},ttv_vec::Array{T,3};it_solver=false,itslv_thresh=256::Int64,maxiter=200::Int64,tol=1e-6::Float64) where T<:Number
	K_dims = (size(Gi,1),size(Gi,2),size(Hi,2))
	if it_solver || prod(K_dims) > itslv_thresh
		H = zeros(T,prod(K_dims))
		function K_matfree(V::AbstractArray{T,1};Gi=Gi::Array{T,5},Hi=Hi::Array{T,3},K_dims=K_dims,H=H::AbstractArray{T,1})
			Hrshp = reshape(H,K_dims)
			@tensoropt((b,c,e,f), Hrshp[a,b,c] = Gi[a,b,d,e,z]*reshape(V,K_dims)[d,e,f]*Hi[z,c,f])
			return H::AbstractArray{T,1}
		end
		r = lobpcg(LinearMap(K_matfree,prod(K_dims);ishermitian = true),false,ttv_vec[:],1;maxiter=maxiter,tol=tol)
		return r.λ[1]::Real, reshape(r.X[:,1],K_dims)::Array{T,3}
	else
		K = K_full(Gi,Hi,K_dims)
		F = eigen(Hermitian(K),1:1)
		return real(F.values[1])::Real,reshape(F.vectors[:,1],K_dims)::Array{T,3}
	end
end

function K_eiggenmin(Gi,Hi,Ki,Li,ttv_vec;it_solver=false,itslv_thresh=2500)
	@tensor begin
		K[a,b,c,d,e,f] := Gi[d,e,a,b,z]*Hi[z,f,c] #size (ni,rim,ri,ni,rim,ri)
		S[a,b,c,d,e,f] := Ki[d,e,a,b,z]*Li[z,f,c] #size (ni,rim,ri,ni,rim,ri)
	end
	if it_solver || prod(size(K)[1:3]) > itslv_thresh
		r = lobpcg(reshape(K,prod(size(K)[1:3]),:),reshape(S,prod(size(S)[1:3]),:),false,ttv_vec[:],1;maxiter=500,tol=1e-8)
		return r.λ[1], reshape(r.X[:,1],size(K)[1:3])
	else
		F = eigen(reshape(K,prod(size(K)[1:3]),:),reshape(S,prod(size(K)[1:3]),:),)
		return real(F.values[1]),reshape(F.vectors[:,1],size(K)[1:3])
	end
end

function left_core_move(x_tt::ttvector{T},V::Array{T,3},i::Int,x_rks) where T<:Number
	rim,ri = x_rks[i],x_rks[i+1]
	ni = x_tt.ttv_dims[i]

	# Prepare core movements
	QV, RV = qr(reshape(permutedims(V, [1 3 2]), ni*ri, :)) #QV: ni*ri x ni*ri; RV ni*ri x rim

	# Apply core movement 3.1
	x_tt.ttv_vec[i] = permutedims(reshape(Matrix(QV)[:, 1:rim], ni, ri, :),[1 3 2])
	x_tt.ttv_ot[i] = 1

	# Apply core movement 3.2
	@tensoropt((b,c,z) , Xim[a,b,c] := x_tt.ttv_vec[i-1][a,b,z]*RV[1:rim,:][c,z]) #size (nim,rim2,rim_new)
	x_tt.ttv_vec[i-1] = Xim
	x_tt.ttv_ot[i-1] = 0
	return x_tt
end

function right_core_move(x_tt::ttvector{T},V::Array{T,3},i::Int,x_rks) where T<:Number
	rim,ri = x_rks[i],x_rks[i+1]
	ni = x_tt.ttv_dims[i]
	QV, RV = qr(reshape(V, ni*rim, :)) #QV: ni*rim x ni*rim; RV ni*rim x ri

	# Apply core movement 3.1
	x_tt.ttv_vec[i] = reshape(Matrix(QV)[:, 1:ri], ni, rim, :)
	x_tt.ttv_ot[i] = -1

	# Apply core movement 3.2
	@tensoropt((b,c,z), Xip[a,b,c] := RV[1:ri,:][b,z]*x_tt.ttv_vec[i+1][a,z,c]) #size (nip,ri,rip)
	x_tt.ttv_vec[i+1] = Xip
	x_tt.ttv_ot[i+1] = 0
	return x_tt
end

"""
Solve Ax=b using the ALS algorithm where A is given as `ttoperator` and `b`, `tt_start` are `ttvector`.
The ranks of the solution is the same as `tt_start`.
`sweep_count` is the number of total sweeps in the ALS.
"""
function als_linsolv(A :: ttoperator{T}, b :: ttvector{T}, tt_start :: ttvector{T} ;sweep_count=2,it_solver=false,r_itsolver=5000) where T<:Number
	# als finds the minimum of the operator J:1/2*<Ax,Ax> - <x,b>
	# input:
	# 	A: the tensor operator in its tensor train format
	#   b: the tensor in its tensor train format
	#	tt_start: start value in its tensor train format
	#	opt_rks: rank vector considered to be optimal enough
	# output:
	#	tt_opt: stationary point of J up to tolerated rank opt_rks
	# 			in its tensor train format

	# Initialize the to be returned tensor in its tensor train format
	tt_opt = orthogonalize(tt_start)
	dims = tt_start.ttv_dims
	d = length(dims)
	# Define the array of ranks of tt_opt [r_0=1,r_1,...,r_d]
	rks = copy(tt_start.ttv_rks)

	# Initialize the arrays of G and G_b
	G = Array{Array{T}}(undef, d)
	G_b = Array{Array{T}}(undef, d)

	# Initialize G[1], G_b[1], H[d] and H_b[d]
	for i in 1:d
		G[i] = zeros(T,dims[i],rks[i],dims[i],rks[i],A.tto_rks[i+1])
		G_b[i] = zeros(dims[i],rks[i],b.ttv_rks[i+1])
	end
	G[1] = reshape(A.tto_vec[1][:,:,1,:], dims[1],1,dims[1], 1, :)
	G_b[1] = reshape(b.ttv_vec[1], dims[1], 1, :)

	#Initialize H and H_b
	H = init_H(tt_opt,A)
	H_b = init_Hb(tt_opt,b)

	nsweeps = 0 #sweeps counter

	while nsweeps < sweep_count
		nsweeps+=1
		# First half sweep
		for i = 1:(d-1)
			println("Forward sweep: core optimization $i out of $d")
			# Define V as solution of K*x=Pb in x
			V = Ksolve(G[i],G_b[i],H[i],H_b[i])
			tt_opt = right_core_move(tt_opt,V,i,rks)
			#update G,G_b
			update_G!(tt_opt.ttv_vec[i],A.tto_vec[i+1],G[i],G[i+1])
			update_Gb!(tt_opt.ttv_vec[i],b.ttv_vec[i+1],G_b[i],G_b[i+1])
		end

		if nsweeps == sweep_count
			return tt_opt
		else
			nsweeps+=1
			# Second half sweep
			for i = d:(-1):2
				println("Backward sweep: core optimization $i out of $d")
				# Define V as solution of K*x=Pb in x
				V = Ksolve(G[i],G_b[i],H[i],H_b[i])
				tt_opt = left_core_move(tt_opt,V,i,rks)
				println(norm(tt_opt.ttv_vec[i-1]))
				update_H!(tt_opt.ttv_vec[i],A.tto_vec[i],H[i],H[i-1])
				update_Hb!(tt_opt.ttv_vec[i],b.ttv_vec[i],H_b[i],H_b[i-1])
			end
		end
	end
	return tt_opt
end

"""
Returns the lowest eigenvalue of A by minimizing the Rayleigh quotient in the ALS algorithm.

The ranks can be increased in the course of the ALS: if `sweep_schedule[k] ≤ i <sweep_schedule[k+1]` is the current number of sweeps then the ranks is given by `rmax_schedule[k]`.
"""
function als_eigsolv(A :: ttoperator{T},
	 tt_start :: ttvector{T} ; #TT initial guess
	 sweep_schedule=[2]::Array{Int64,1}, #Number of sweeps for each bond dimension in rmax_schedule
	 rmax_schedule=[maximum(tt_start.ttv_rks)]::Array{Int64,1}, #bond dimension at each sweep
	 noise_schedule=zeros(length(rmax_schedule))::Array{Float64,1}, #noise at each bond dimension increase
	 it_solver=false::Bool, #linear solver for the microstep
	 itslv_thresh=1024::Int64, #switch from full to iterative
	 maxiter=200::Int64, #maximum of iterations for the iterative solver
	 linsolv_tol=1e-8::Float64) where T <: Number #tolerance of the iterative linear solver
	@assert(length(rmax_schedule)==length(sweep_schedule)==length(noise_schedule),"Sweep schedule error")

	# Initialize the to be returned tensor in its tensor train format
	tt_opt = orthogonalize(tt_start)
	dims = tt_start.ttv_dims
	d = length(dims)
	E = zeros(Float64,2d*(sweep_schedule[end]+1)) #output eigenvalue
	# Define the array of ranks of tt_opt [r_0=1,r_1,...,r_d]
	rks = copy(tt_start.ttv_rks)

	# Initialize the array G
	G = Array{Array{T}}(undef, d)
	for i in 1:d
		G[i] = zeros(T,dims[i],rks[i],dims[i],rks[i],A.tto_rks[i+1])
	end
	G[1] = reshape(A.tto_vec[1][:,:,1,:], dims[1], 1, dims[1], 1, :)

	#Initialize H and H_b
	H = init_H(tt_opt,A)

	nsweeps = 0 #sweeps counter
	i_schedule,i_μit = 1,0
	while i_schedule <= length(sweep_schedule)
		nsweeps+=1
		println("Macro-iteration $nsweeps; bond dimension $(rmax_schedule[i_schedule])")
		if nsweeps == sweep_schedule[i_schedule]
			i_schedule+=1
			if i_schedule > length(sweep_schedule)
				return E[1:i_μit]::Array{Float64,1},tt_opt::ttvector{T}
			else
				tt_opt = tt_up_rks(tt_opt,rmax_schedule[i_schedule];ϵ_wn=noise_schedule[i_schedule])
				tt_opt = orthogonalize(tt_opt)
				H = init_H(tt_opt,A)
				for i in 1:d-1
					Gtemp = zeros(dims[i+1],tt_opt.ttv_rks[i+1],dims[i+1],tt_opt.ttv_rks[i+1],A.tto_rks[i+2])
					Gtemp[1:size(G[i+1],1),1:size(G[i+1],2),1:size(G[i+1],3),1:size(G[i+1],4),1:size(G[i+1],5)] = G[i+1]
					G[i+1] = Gtemp
				end
			end
		end
		# First half sweep
		for i = 1:(d-1)
			println("Forward sweep: core optimization $i out of $(d-1)")
			# Define V as solution of K*x=Pb in x
			i_μit += 1
			E[i_μit],V = K_eigmin(G[i],H[i],tt_opt.ttv_vec[i];it_solver=it_solver,itslv_thresh=itslv_thresh,maxiter=maxiter,tol=linsolv_tol)
			println("Eigenvalue: $(E[i_μit])")
			tt_opt = right_core_move(tt_opt,V,i,tt_opt.ttv_rks)
			#update G
			update_G!(tt_opt.ttv_vec[i],A.tto_vec[i+1],G[i],G[i+1])
		end

		# Second half sweep
		for i = d:(-1):2
			println("Backward sweep: core optimization $(d+1-i) out of $(d-1)")
			# Define V as solution of K*x=Pb in x
			i_μit += 1
			E[i_μit],V = K_eigmin(G[i],H[i],tt_opt.ttv_vec[i];it_solver=it_solver,itslv_thresh=itslv_thresh,maxiter=maxiter,tol=linsolv_tol)
			println("Eigenvalue: $(E[i_μit])")
			tt_opt = left_core_move(tt_opt,V,i,tt_opt.ttv_rks)
			update_H!(tt_opt.ttv_vec[i],A.tto_vec[i],H[i],H[i-1])
		end
	end
	return E[1:i_μit]::Array{Float64,1},tt_opt::ttvector{T}
end

"""
returns the smallest eigenpair Ax = Sx
"""
function als_gen_eigsolv(A :: ttoperator{T}, S::ttoperator{T}, tt_start :: ttvector{T} ; sweep_schedule=[2],rmax_schedule=[maximum(tt_start.ttv_rks)],tol=1e-10,it_solver=false,itslv_thresh=2500) where T<:Number
	# Initialize the to be returned tensor in its tensor train format
	tt_opt = orthogonalize(tt_start)
	dims = tt_start.ttv_dims
	d = length(dims)
	E = zeros(Float64,d*sweep_schedule[end]) #output eigenvalue
	# Define the array of ranks of tt_opt [r_0=1,r_1,...,r_d]
	rks = tt_start.ttv_rks

	# Initialize the arrays of G and K
	G = Array{Array{T}}(undef, d)
	K = Array{Array{T}}(undef, d)

	# Initialize G[1]
	for i in 1:d
		G[i] = zeros(dims[i],rks[i],dims[i],rks[i],A.tto_rks[i+1])
		K[i] = zeros(dims[i],rks[i],dims[i],rks[i],S.tto_rks[i+1])
	end
	G[1] = reshape(A.tto_vec[1][:,:,1,:], dims[1],1,dims[1], 1, :)
	K[1] = reshape(S.tto_vec[1][:,:,1,:], dims[1],1,dims[1], 1, :)

	#Initialize H and H_b
	H = init_H(tt_opt,A)
	L = init_H(tt_opt,S)

	nsweeps = 0 #sweeps counter
	i_schedule,i_μit = 1,0
	while i_schedule <= length(sweep_schedule)
		nsweeps+=1
		if nsweeps == sweep_schedule[i_schedule]
			i_schedule+=1
			if i_schedule > length(sweep_schedule)
				return E[1:i_μit],tt_opt
			else
				tt_opt = tt_up_rks(tt_opt,rmax_schedule[i_schedule])
				for i in 1:d-1
					Htemp = zeros(tt_opt.ttv_rks[i],tt_opt.ttv_rks[i],A.tto_rks[i])
					Ltemp = zeros(tt_opt.ttv_rks[i],tt_opt.ttv_rks[i],S.tto_rks[i])
					Htemp[1:size(H[i],1),1:size(H[i],2),1:size(H[i],3)] = H[i]
					Ltemp[1:size(L[i],1),1:size(L[i],2),1:size(L[i],3)] = L[i]
					H[i] = Htemp
					L[i] = Ltemp
				end
			end
		end

		# First half sweep
		for i = 1:(d-1)
			println("Forward sweep: core optimization $i out of $d")

			# If i is the index of the core matrices do the optimization
			if tt_opt.ttv_ot[i] == 0
				# Define V as solution of K*x=Pb in x
				i_μit += 1
				E[i_μit],V = K_eiggenmin(G[i],H[i],K[i],L[i],tt_opt.ttv_vec[i];it_solver=it_solver,itslv_thresh=itslv_thresh)
				println("Eigenvalue: $(E[i_μit])")
				tt_opt = right_core_move(tt_opt,V,i,rks)
			end

			#update G and K
			update_G!(tt_opt.ttv_vec[i],A.tto_vec[i+1],G[i],G[i+1])
			update_G!(tt_opt.ttv_vec[i],S.tto_vec[i+1],K[i],K[i+1])
		end

		# Second half sweep
		for i = d:(-1):2
			println("Backward sweep: core optimization $i out of $d")
			# Define V as solution of K*x=Pb in x
			i_μit += 1
			E[i_μit],V = K_eiggenmin(G[i],H[i],K[i],L[i],tt_opt.ttv_vec[i];it_solver=it_solver,itslv_thresh=itslv_thresh)
			println("Eigenvalue: $(E[i_μit])")
			tt_opt = left_core_move(tt_opt,V,i,rks)
			update_H!(tt_opt.ttv_vec[i],A.tto_vec[i],H[i],H[i-1])
			update_H!(tt_opt.ttv_vec[i],S.tto_vec[i],L[i],L[i-1])
		end
	end
end

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
using CompScienceMeshes

export meshwaveguidepost

"""
    meshwaveguidepost(width, height, length, post_radius, delta)

Creates a mesh for a waveguide of width (along the x-axis) `width`, height (along
    the y-axis) `height` and length (along the z-axis) `length` by parsing a .geo script
    incorporating these parameters into the GMSH mesher.

A cylindrical post of radius `post_radius` is positioned vertically (i.e.
    parallel to y-axis) midway along the waveguide.

The target edge size is `delta`.
"""

function meshwaveguidepost(width, height, length, post_radius, delta)
    s =
"""

lc = $delta;

//--------------WAVEGUIDE POINTS--------------

Point(1) = {0.0, 0.0, 0.0, lc};
Point(2) = {0.0, $height, 0.0, lc};
Point(3) = {$width, $height, 0.0, lc};
Point(4) = {$width, 0.0, 0.0, lc};

Point(5) = {0.0, 0.0, $length, lc};
Point(6) = {0.0, $height, $length, lc};
Point(7) = {$width, $height, $length, lc};
Point(8) = {$width, 0.0, $length, lc};

//--------------POST POINTS-------------------

Point(9) = {($width)/2,0,($length)/2,lc};

Point(10) = {($width)/2,0,(($length)/2)-(($post_radius)/2),lc};
Point(11) = {(($width)/2)-(($post_radius)/2),0,($length)/2,lc};
Point(12) = {($width)/2,0,(($length)/2)+(($post_radius)/2),lc};
Point(13) = {(($width)/2)+(($post_radius)/2),0,($length)/2,lc};

Point(14) = {($width)/2,$height,($length)/2,lc};

Point(15) = {($width)/2,$height,(($length)/2)-(($post_radius)/2),lc};
Point(16) = {(($width)/2)-(($post_radius)/2),$height,($length)/2,lc};
Point(17) = {($width)/2,$height,(($length)/2)+(($post_radius)/2),lc};
Point(18) = {(($width)/2)+(($post_radius)/2),$height,($length)/2,lc};

//--------------WAVEGUIDE-POST CONTINUOUS SURFACE DEFINITION---------

//Part 1 - waveguide boundaries

Line(1)={1,2};
Line(2)={2,3};
Line(3)={3,4};
Line(4)={4,1};

Line(5)={5,6};
Line(6)={6,7};
Line(7)={7,8};
Line(8)={8,5};

Line(9)={1,5};
Line(10)={2,6};
Line(11)={3,7};
Line(12)={4,8};

//Part 2 - discontinuity boundaries

Circle(13) = {10,9,11};
Circle(14) = {11,9,12};
Circle(15) = {12,9,13};
Circle(16) = {13,9,10};

Circle(17) = {15,14,16};
Circle(18) = {16,14,17};
Circle(19) = {17,14,18};
Circle(20) = {18,14,15};

Line(21) = {10,15};
Line(22) = {11,16};
Line(23) = {12,17};
Line(24) = {13,18};

//Part 3 - creating skeleton

Line Loop(1)={-1,9,5,-10};

Line Loop(2)={-2,10,6,-11};
Line Loop(3)={17,18,19,20};

Line Loop(4)={-3,11,7,-12};

Line Loop(5)={-4,12,8,-9};
Line Loop(6)={-16,-15,-14,-13};

Line Loop(7) = {13,22,-17,-21};
Line Loop(8) = {14,23,-18,-22};
Line Loop(9) = {15,24,-19,-23};
Line Loop(10) = {16,21,-20,-24};

//Part 4 - creating surfaces

Plane Surface(1)={1};
Plane Surface(2)={3,2};
Plane Surface(3)={4};
Plane Surface(4)={6,5};

Ruled Surface(5) = {7};
Ruled Surface(6) = {8};
Ruled Surface(7) = {9};
Ruled Surface(8) = {10};

"""

    fn = tempname()
    io = open(fn, "w")
    try
        print(io, s)
    finally
        close(io)
    end

    # feed the file to gmsh
    fno = tempname()
    run(`gmsh $fn -2 -format msh -o $fno`)
    fdo = open(fno,"r")
    # include(Pkg.dir("CompScienceMeshes","src","gmsh.jl"))
    m = CompScienceMeshes.read_gmsh_mesh(fdo)
    close(fdo)
    rm(fno)
    rm(fn)

    return m

end

# Γ = meshwaveguidepost(22,22,50,3.81,0.8)
# RT = raviartthomas(Γ)
#
# include(Pkg.dir("CompScienceMeshes","examples","plotlyjs_patches.jl"))
# A = rand(numcells(Γ))
# p = patch(Γ, A)
# PlotlyJS.plot([p])


text:
"""
This script generates micromamba artifacts.

Update the parameters (usually just the version) and run `julia src/generate_artifacts.jl`
from the root of the package. Requires ArtifactUtils to be installed.

The latest version and supported platforms can be found at
https://anaconda.org/conda-forge/micromamba
"""

import ArtifactUtils: add_artifact!
import Pkg.Artifacts: Platform

toml = "Artifacts.toml"
version = "0.24.0"
name = "micromamba-$(version)"
platforms = [
    "linux-64" => Platform("x86_64", "linux"),
    "linux-aarch64" => Platform("aarch64", "linux"),
    "linux-ppc64le" => Platform("powerpc64le", "linux"),
    "osx-64" => Platform("x86_64", "macos"),
    "osx-arm64" => Platform("aarch64", "macos"),
    "win-64" => Platform("x86_64", "windows"),
]

@show toml version name

for platform in platforms
    @show platform
    url = "https://micro.mamba.pm/api/micromamba/$(platform[1])/$(version)"
    add_artifact!(toml, name, url; platform=platform[2], lazy=true, clear=true)
end


text:
module Day25

export get_inputs, get_solution1, get_solution2
export transform, find_loop_size

## Input getting
function get_inputs()
    test_input1 = [5764801, 17807724]
    test_output1 = 14897079
    test_input2 = nothing
    test_output2 = nothing
    data = parse.(Int, readlines(joinpath(@__DIR__, "input.txt")))
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

## Solution functions
function transform(subject_number, loop_size)
    value = 1
    for i in 1:loop_size
        value = rem(value * subject_number, 20201227)
    end
    return value
end

function find_loop_size(subject_number, public_key)
    value = 1
    loop_size = 0
    while value != public_key
        value = rem(value * subject_number, 20201227)
        loop_size += 1
    end
    return loop_size
end

# Part 1
function get_solution1(data)
    card_public_key, door_public_key = data
    card_loop_size, door_loop_size = find_loop_size.(7, data)
    return transform(door_public_key, card_loop_size)
end

# Part 2
function get_solution2(data)
    return nothing
end

end


text:
"""
***
# Viscosity of Liquids
***
Data are calculated by solving the equation:

``Log_{10}(\\mu_{G})= A+B/T+C*T+D*T^2``

This function can be used by two methods:
* Get the liquid viscosity by introducing the component's name and the temperature.
* Get the full list of components' names supported.
***
### Example:
`julia> LiquidViscosity("Water",298.15)`

`0.911`

`julia> LiquidViscosity()`

`String["Carbon monoxide", "1,3-Butadiene", "Cyclohexanol", "Styrene", "n-Octane", "Acetic acid", "Hydrogen", "Dichloromethane",
"1-Butanol (n-Butanol)", "1-Octene", "Methyl chloride", "Ethylene glycol", "Acrylic acid", "Ethanol", "1,1,1-Trichloroethane", "Oxygen",
"1,4-Dioxane", "n-Nonane", "Pyridine", "Hydrogen cyanide", "Ammonia", "carbon disulfide", "1-Nonane (n-Nonane)", "Hydrogen chloride",
"1-Decene", "Chlorine", "Acetic anhydride", "1,1-Dichloroethane", "Methyl isobutyl ketone", "Chloroform", "1,1,2-Trichloroethane",
 "p-Xylene", "Butyric acid", "Nitrobenzene", "1-Butene", "Ethylene oxide", "Cyclohexane", "Acetaldehyde", "1,2-Dichloroethane",
 "Methane", "n-Heptane", "n-Hexane", "Nitrogen", "Carbon tetrachloride", "Hydrogen peroxide", "n-Decane", "Ethane", "Water",
 "Diethyl ketone", "Ethyl chloride", "Naphthalene", "Fluorine", "1-Propanol (n-Propanol)", "n-Butane", "o-Xylene", "Toluene",
 "n-Pentane", "Aniline", "Methyl bromide", "m-Xylene", "Carbon dioxide", "Benzene", "Phenol", "i-Butane (iso-Butane)", "Methylamine",
 "Dimethylamine", "Acetone", "Ethyl acetate", "Ethylbenzene", "Cyclopropane", "Glycerol", "Sulfur dioxide", "Ethylene", "Diethyl ether", "Vinyl acetate"]`
"""
function LiquidViscosity(substance::String,T::Float64)
    try
        param=DataLiquidViscosity[substance]

        if T >= param[6] || T <= param[5]
            print_with_color(:red,
            @sprintf("%4.4f K is out of range %4.4f K - %4.4f K.\n",T,
            param[5],param[6]))
            return NaN
        else
            μG=10^(param[1]+param[2]/T+param[3]*T+param[4]*T^2)
        end
    catch
        print_with_color(:red,@sprintf("%s not found in database.\n",substance))
        return NaN
    end
end

function LiquidViscosity()
    println(keys(DataLiquidViscosity))
end


text:
export Fourier

struct Fourier <: AbstractMethod

    kx   :: Vector{Float64}
    kv   :: Vector{Float64}

    function Fourier( meshx, meshv )

        nx  = meshx.len
        dx  = meshx.step
        Lx  = meshx.stop - meshx.start
        kx  = 2π/Lx .* [0:nx÷2-1;-nx÷2:-1]

        nv  = meshv.len
        dv  = meshv.step
        Lv  = meshv.stop - meshv.start
        kv  = 2π/Lv .* [0:nv÷2-1;-nv÷2:-1]

        new( kx, kv)

    end

end

function advection_v!(fᵗ  :: Array{ComplexF64,2},
                      adv :: Fourier,
		      e   :: Vector{ComplexF64},
		      dt  :: Float64 )
    fft!(fᵗ, 1)
    @. fᵗ *= exp(-1im * dt * adv.kv * $transpose(e))
    ifft!(fᵗ, 1)

end

function advection_x!( f   :: Array{ComplexF64,2},
                       adv :: Fourier,
		       e   :: Vector{ComplexF64},
		       v   :: Vector{Float64},
		       dt  :: Float64 )

    ev = exp.(-1im*dt * adv.kx * transpose(v))

    fft!(f,1)
    f .= f .* ev
    dv = v[2]-v[1]
    ρ = dv * vec(sum(f,dims=2))
    for i in 2:length(e)
        e[i] = -1im * ρ[i] ./ adv.kx[i]
    end
    e[1] = 0.0
    ifft!(f,1)
    ifft!(e)
    e .= real(e)
end

max_stars_repo_name: JuliaVlasov/VlasovSolvers

max_issues_repo_name: JuliaVlasov/VlasovSolvers

max_forks_repo_name: JuliaVlasov/VlasovSolvers


text:
# This code calls the first-order eccentric TTV code
# compute_ttv.jl.  Please cite Agol & Deck (2015) if
# you make use of this in your research.

include("compute_ttv.jl")
#include("laplace_coefficients_initialize.jl")
#alpha0 may be wrong - Noah
function test_ttv{T1<:Number,T2<:Number}(jmax::Integer,n1::Integer,n2::Integer,data::Vector{T1}; WriteOutput::Bool = true, num_evals::Integer = 1, profile::Bool = false, alpha0::T2 = abs(data[2]/data[7])^(2/3), b0::Array{T2,2} = TTVFaster.LaplaceCoefficients.initialize(jmax+1,alpha0) )
@assert(jmax>=1)  # Should there be a larger minimum?

@assert(n1>2)
@assert(n2>2)
@assert(length(data)==10)
# Performs a test of the transit_times.jl routine
# Set up planets planar-planet types for the inner and outer planets:
p1=TTVFaster.Planet_plane_hk(data[1],data[2],data[3],data[4],data[ 5])
p2=TTVFaster.Planet_plane_hk(data[6],data[7],data[8],data[9],data[10])
time1 = zeros(eltype(data),n1)
for i in 1:length(time1)
   time1[i] = p1.trans0 + i * p1.period
end
time2 = zeros(eltype(data),n2)
for i in 1:length(time2)
   time2[i] = p2.trans0 + i * p2.period
end
#time1 = collect(p1.trans0 + collect(linspace(0,n1-1,n1)) * p1.period)
#time2 = collect(p2.trans0 + collect(linspace(0,n2-1,n2)) * p2.period)
# Moved alpha0 and b0 to optional parameters, so don't need to recompute every time
# alpha0=abs(p1.period/p2.period)^(2//3)
# Initialize the computation of the Laplace coefficients:
#b0=TTVFaster.LaplaceCoefficients.initialize(jmax+1,alpha0)
# Define arrays to hold the TTVs:
ttv1=Array(T1,n1)
ttv2=Array(T1,n2)
# Define arrays to hold the TTV coefficients and Laplace coefficients:
f1=Array(T1,jmax+2,5)
f2=Array(T1,jmax+2,5)
b=Array(T1,jmax+2,3)
hashsum = 0
for i in 1:num_evals
   # Call the compute_ttv code which implements equation (33)
   TTVFaster.compute_ttv!(jmax,p1,p2,time1,time2,ttv1,ttv2,f1,f2,b,alpha0,b0)
   if profile
     hashsum += hash(ttv1)+hash(ttv2)
   end
end
if profile
  println("# Ignore this: ",hashsum) # This just makes sure optimizer doesn't optimize away important calculations.
end

if WriteOutput
   # Write the mean ephemeris and TTV results to two files:
   writedlm("inner_ttv.txt",[time1 ttv1])
   writedlm("outer_ttv.txt",[time2 ttv2])
end
   return ttv1,ttv2
end

max_stars_repo_name: nwtuchow/TTVFaster

max_issues_repo_name: nwtuchow/TTVFaster

max_forks_repo_name: nwtuchow/TTVFaster


text:
__precompile__()

module FixedLengthStrings

export SString, sstring
export @tiny_string, @short_string, @mid_string, @long_string, @text_string

using StaticArrays

using Compat
import Compat.String

struct SString{n,T <: Union{Char, UInt8}} <: AbstractString
    sv::SVector{n,T}
end

function _string(s::SString)
    replace(
        convert(String, s.sv),
        '\0', ""
    )
end

# coercing constructor:
function sstring(s::String, n=nothing)
    str_length = length(s)

    ( n == nothing ) && ( n = str_length )

    ( str_length > n ) &&
        error("String too large to fit into a fixed length of $n")

    nul_length = n - length(s)

    cur_array = eval(parse("b\"$(s)\""))

    cur_nul = eval(parse("$(eltype(cur_array))(C_NULL)"))

    nul_array = repmat([cur_nul], nul_length)

    append!(cur_array, nul_array)

    sv = SVector{n,Char}(cur_array)

    return SString(sv)
end

sstring(s::AbstractString, n=nothing) = sstring(String(s), n)
sstring(args...; n=nothing) = sstring(string(args...), n)

macro tiny_string(s, flags...) sstring(s, n=4) end
macro short_string(s, flags...) sstring(s, n=16) end
macro mid_string(s, flags...) sstring(s, n=64) end
macro long_string(s, flags...) sstring(s, n=256) end
macro text_string(s, flags...) sstring(s, n=1024) end

import Base: write, endof, getindex, sizeof, search, rsearch, isvalid, next, length, IOBuffer, pointer

if isdefined(Base, :bytestring)
    import Base: bytestring
    bytestring(s::SString) = bytestring(_string(s))
end

endof(s::SString) = endof(_string(s))
next(s::SString, i::Int) = next(_string(s), i)
length(s::SString) = length(_string(s))
getindex(s::SString, i::Int) = getindex(_string(s), i)
getindex(s::SString, i::Integer) = getindex(_string(s), i)
getindex(s::SString, i::Real) = getindex(_string(s), i)
getindex(s::SString, i::UnitRange{Int}) = getindex(_string(s), i)
getindex{T<:Integer}(s::SString, i::UnitRange{T}) = getindex(_string(s), i)
getindex(s::SString, i::AbstractVector) = getindex(_string(s), i)
sizeof(s::SString) = sizeof(_string(s))
search(s::SString, c::Char, i::Integer) = search(_string(s), c, i)
rsearch(s::SString, c::Char, i::Integer) = rsearch(_string(s), c, i)
isvalid(s::SString, i::Integer) = isvalid(_string(s), i)
pointer(s::SString) = pointer(_string(s))
IOBuffer(s::SString) = IOBuffer(_string(s))

import Base.convert
const unsafe_convert = Base.convert

@compat unsafe_convert(T::Union{Type{Ptr{UInt8}},Type{Ptr{Int8}}}, s::SString) = convert(T, _string(s))
unsafe_convert(::Type{Cstring}, s::SString) = unsafe_convert(Cstring, _string(s))

end # module


text:
function localmaxmin(y::Array{Float64,1},nn::Int64)
    dummy_min = Array{Int32}(undef,nn,1)
    dummy_max = Array{Int32}(undef,nn,1)
    dummy_ids = Array{Int32}(undef,nn,1)
    fill!(dummy_min,0)
    fill!(dummy_max,0)
    fill!(dummy_ids,0)

    # ---- Check 1st node
    cmin,cmax,cids = 1,1,1
    if y[2] > y[1] # local minima
        dummy_min[1] = 1
        cmin        += 1
        dummy_ids[1] = 1
        cids        += 1
    else # local maxima
        dummy_max[1] = 1
        cmax        += 1
        dummy_ids[1] = 1
        cids        += 1
    end

    # ---- Check 2:end-1 node
    @inbounds @simd for ii = 2:nn-2
        if y[ii+1] > y[ii] && y[ii-1] > y[ii] # local minima
            dummy_min[cmin] = ii
            cmin    += 1
            dummy_ids[cids] = ii
            cids    += 1
        elseif y[ii+1] < y[ii] && y[ii-1] < y[ii] # local maxima
            dummy_max[cmax] = ii
            cmax    += 1
            dummy_ids[cids] = ii
            cids    += 1
        end
    end

    # ---- Check last node
    if y[nn] < y[nn-1] # local minima
        dummy_min[cmin] = nn
        cmin        += 1
        dummy_ids[cids] = nn
        cids        += 1
    else # local maxima
        dummy_max[cmax] = nn
        cmax        += 1
        dummy_ids[cids] = nn
        cids        += 1
    end

    # ---- Remove zeros
    idx_min = Array{Int32}(undef, cmin-1,1)
    idx_max = Array{Int32}(undef, cmax-1,1)
    idx_ids = Array{Int32}(undef, cids-1,1)
    @inbounds for ii = 1:cmin-1
        idx_min[ii] = dummy_min[ii]
    end
    @inbounds for ii = 1:cmax-1
        idx_max[ii] = dummy_max[ii]
    end
    @inbounds for ii = 1:cids-1
        idx_ids[ii] = dummy_ids[ii]
    end

    return idx_min,idx_max,idx_ids
end

function segdist(idx_ids::Array{Int32,2},x::Vector{Float64},z::Vector{Float64})
    nx        = length(x)
    nsegments = length(idx_ids) - 1
    Lsegments = Array{Float64}(undef,nsegments,1)
    dist2high = Array{Float64}(undef,nx,1)

    @inbounds for ii = 1:nsegments
        id1 = idx_ids[ii];
        id2 = idx_ids[ii+1];

        for jj = id1:id2

            if  z[id2] > z[id1]
                xhigh = x[id2]
            else
                xhigh = x[id1]
            end

            dist2high[jj] = abs(x[jj] - xhigh)

        end
    end

    return dist2high

end ## END segdist FUNCTION

function firstderivative(x,z) # 1st derivative
    dz = (z[2] - z[1])*2
    ∇z = similar(x)

    @avx for i ∈ 2:length(x)-1
        ∇z[i] = (x[i+1] - x[i-1]) / dz
    end

    ∇z[1]   = (x[2] - x[1]) / (dz)
    ∇z[end] = (x[end] - x[end-1]) / (dz)

    return ∇z
end

function secondderivative(x,z) # 2nd derivative
    dz = 2*(z[2] - z[1])^2
    Δz = similar(x)

    @avx for i ∈ 2:length(x)-1
        Δz[i] = (x[i+1] - 2*x[i] + x[i-1]) / dz
    end

    Δz[1]   = Δz[2]
    Δz[end] = Δz[end-1]

    return Δz
end

function erodedvolume(z0,z,x)
    δz = z - z0
    return integrate(x,δz,TrapezoidalEvenFast())
end

# NOTE: thigs should be in-place
function riverlength(x, z, nnod)
    L = fill(0.0,nnod)
    # cummulative length
    @tturbo for i ∈ 2:nnod
        L[i] = L[i-1] + sqrt((x[i]-x[i-1])^2  + (z[i]-z[i-1])^2)
    end
    return L
end

@inline pow32(a::T, b::Real) where T =
    convert(T, convert(Float32,a)^convert(Float32,b))

function uplift!(z1, z2, h)
    @tturbo for i in eachindex(z1)
        z1[i] += h
        z2[i] += h
    end
end


text:
using PyCall

@pyimport brpylib

# file downloaded from https://portal.g-node.org/neo/blackrock/
datafile = "FileSpec2.3001.ns5"

# Open file and extract headers
brns_file = brpylib.NsxFile(datafile)

max_stars_repo_name: paulstey/blackrockio

max_issues_repo_name: paulstey/blackrockio

max_forks_repo_name: paulstey/blackrockio


text:
# Setup:

isbinder() = "jovyan" in split(pwd(), "/")

const REPO = "https://github.com/ablaom/MachineLearningInJulia2020"
using Pkg

if !isbinder()
    Pkg.activate(DIR)
    Pkg.instantiate()
    using CategoricalArrays
    import MLJLinearModels
    import DataFrames
    import CSV
    import DecisionTree
    using MLJ
    import MLJLinearModels
    import MultivariateStats
    import MLJFlux
    import Plots
else
    @info "Skipping package instantiation as binder notebook. "
end
@info "Done loading"


text:
###############################################################################
#
# RT Dimensions
#
###############################################################################
"""
    mutable struct RTDimensions

Struct that stores matrix dimension information

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef mutable struct RTDimensions
    "Number of azimuth angles"
    nAzi  ::Int = 36
    "Number of inclination agles"
    nIncl ::Int = 9
    "Number of canopy layers"
    nLayer::Int = 5
    "Number of canopy layer boundaries nLayer+1"
    nLevel::Int = nLayer+1
    "Number of PAR wave lengths"
    nPAR  ::Int = 10
    "Number of wave lengths"
    nWL   ::Int = 10
    "Number of wave length for excitation"
    nWLE  ::Int = 10
    "Number of wave lengths for SIF"
    nWLF  ::Int = 10
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля
#РЕЗУЛЬТАТ: Робот - в исходном положении, и клетки поля промакированы так: нижний ряд - полностью,
#следующий - весь, за исключением одной последней клетки на Востоке, следующий - за исключением двух
#последних клеток на Востоке, и т.д.

function glavnaia(r::Robot)
    side=Ost
    num_ver=moves!(r,Sud)
    num_hor=moves!(r,Ost)
    count=moves!(r,West)
    num_hor=count-num_hor
    while count>0
        putmarker!(r)
        for _ in 1:count
            move!(r,side)
            putmarker!(r)
        end
        if isborder(r,West)==false
            move!(r,West)
        end
        if isborder(r,Nord)==false
            move!(r,Nord)
            num_ver-=1
        else
            break
        end
        side=inverse(side)
        count-=1
    end
    putmarker!(r)
    moves!(r,West)
    move_back(r,Nord,num_ver)
    move_back(r,Ost,num_hor)
end

function moves!(r::Robot,side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

function move_back(r::Robot,side::HorizonSide,num_steps::Int)
    if (num_steps<0)==true
        side=inverse(side) #обратно
        num_steps*=-1
    end
    for _ in 1:num_steps
        move!(r,side)
    end
end
# Возвращает направление, противоположное заданному
inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2,4))

max_stars_repo_name: Aminulik/mirea-progs

max_issues_repo_name: Aminulik/mirea-progs

max_forks_repo_name: Aminulik/mirea-progs


text:
"""
    gxmake(xmin, iwt, n, nxin, iord)

Define a logarithmically-spaced `x` grid.

# Arguments
- `xmin::Array{Float64,1}`: an input array containitn `n` values
of `x` in ascending order. `xmin[1]` defines the lower end of the grid
and other values define approx positions where the point density
will change according to the values set in `iwt`.
- `iwt::Array{Int32,1}`: input integer weights given in ascending order
and must always be an integer multiple of the previous weight.
- `n::Integer`: number of values specified in `xmin` and `iwt`.
- `nxin::Integer`: Requested number of grid points.
- `iord::Integer`: iord = 2(3) for linear(quadratic) spline interpolation.

# Returns
- `nxout::Integer`: the number of generated grid points.
"""
function gxmake(xmin::Array{Float64,1}, iwt::Array{Int32,1}, n::Integer, nxin::Integer, iord::Integer)

    n = Ref{Int32}(n)
    nxin = Ref{Int32}(nxin)
    iord = Ref{Int32}(iord)
    nxout = Ref{Int32}()

    @ccall gxmake_(xmin::Ref{Float64}, iwt::Ref{Int32}, n::Ref{Int32},
                   nxin::Ref{Int32}, nxout::Ref{Int32}, iord::Ref{Int32})::Nothing

    nxout[]
end

"""
    gqmake(qarr, wgt, n, nqin)

Define a logarithmically-spaced `mu_F^2` grid on which the parton
densities are evolved.

# Arguments
- `qarr::Array{Float64,1}`: input array containing `n` values of `mu^2`
in ascending order. The lower edge of the grid should be above 0.1 GeV^2.
- `wgt::Array{Float64,1}`: relative grip point density in each region
defined by `qarr`.
- `n::Integer`: number of values in `qarr` and `wgt` (n>=2).
- `nqin::Integer`: requested number of grid points.

# Returns
- `nqout::Integer`: number of generated grid points.
"""
function gqmake(qarr::Array{Float64,1}, wgt::Array{Float64,1}, n::Integer, nqin::Integer)

    n = Ref{Int32}(n)
    nqin = Ref{Int32}(nqin)
    nqout = Ref{Int32}()

    @ccall gqmake_(qarr::Ref{Float64}, wgt::Ref{Float64}, n::Ref{Int32},
                   nqin::Ref{Int32}, nqout::Ref{Int32})::Nothing

    nqout[]
end

"""
    ixfrmx(x)

Get grid point index of closest grid point at or
below `x` value.
"""
function ixfrmx(x::Float64)

    x = Ref{Float64}(x)

    ix = @ccall ixfrmx_(x::Ref{Float64})::Int32

    ix[]
end

"""
    xfrmix(ix)

Get `x` value at grid point `ix`.
"""
function xfrmix(ix::Integer)

    ix = Ref{Int32}(ix)

    x = @ccall xfrmix_(ix::Ref{Int32})::Float64

    x[]
end

"""
    iqfrmq(q2)

Get grid point index of the closest grid point at or
below `q2` value.
"""
function iqfrmq(q2::Float64)

    q2 = Ref{Float64}(q2)

    iq = @ccall iqfrmq_(q2::Ref{Float64})::Int32

    iq
end

"""
    qfrmiq(iq)

Get `q2` value at grid point `iq`.
"""
function qfrmiq(iq::Integer)

    iq = Ref{Int32}(iq)

    q2 = @ccall qfrmiq_(iq::Ref{Int32})::Float64

    q2[]
end

"""
    grpars()

Get the current grid definitions.

# Returns
`nx::Integer`: number of points in x grid.
`xmi::Float64`: lower boundary of x grid.
`xma::Float64`: upper boundary of x grid.
`nq::Integer`: number of points in qq grid.
`qmi::Float64`: lower boundary of qq grid.
`qma::Float64`: upper boundary of qq grid.
`iord::Integer`: order of spline interpolation.
"""
function grpars()

    nx = Ref{Int32}()
    xmi = Ref{Float64}()
    xma = Ref{Float64}()
    nq = Ref{Int32}()
    qmi = Ref{Float64}()
    qma = Ref{Float64}()
    iord = Ref{Int32}()

    @ccall grpars_(nx::Ref{Int32}, xmi::Ref{Float64}, xma::Ref{Float64},
                   nq::Ref{Int32}, qmi::Ref{Float64}, qma::Ref{Float64},
                   iord::Ref{Int32})::Nothing

    nx[], xmi[], xma[], nq[], qmi[], qma[], iord[]
end

"""
    gxcopy(n)

Copy the current `x` grid into an array of length `n`.
"""
function gxcopy(n::Integer)

    n = Ref{Int32}(n)
    array = Array{Float64}(undef, n[])
    nx = Ref{Int32}()

    @ccall gxcopy_(array::Ref{Float64}, n::Ref{Int32}, nx::Ref{Int32})::Nothing

    array
end

"""
    gqcopy(n)

Copy the current `mu^2` grid into an array of length `n`.
"""
function gqcopy(n::Integer)

    n = Ref{Int32}(n)
    array = Array{Float64}(undef, n[])
    nx = Ref{Int32}()

    @ccall gqcopy_(array::Ref{Float64}, n::Ref{Int32}, nx::Ref{Int32})::Nothing

    array
end


text:
# define constants here -
const alphavantage_api_url_string = "https://www.alphavantage.co/query"
const path_to_package = dirname(pathof(@__MODULE__))

# load offical packages here -
using DataFrames
using CSV
using HTTP
using JSON
using Dates
using Logging

# load our code here -
include("./base/Network.jl")
include("./base/Utility.jl")
include("./base/User.jl")
include("./base/Error.jl")
include("./base/Handlers.jl")
include("./base/Log.jl")

# STS code includes go here -
include("./sts_functions/STSweekly.jl")
include("./sts_functions/STSDaily.jl")
include("./sts_functions/STSMonthly.jl")
include("./sts_functions/STS_adjusted_daily.jl")

max_stars_repo_name: teamalphavantage/AlphaVantageAPI

max_issues_repo_name: teamalphavantage/AlphaVantageAPI

max_forks_repo_name: teamalphavantage/AlphaVantageAPI


text:
"""Binary Builders
JuliaCon2018 -
Video: https://www.youtube.com/watch?v=2e0PBGSaQaI&ab_channel=TheJuliaProgrammingLanguage
Credits: Elliot `Eli` Saba @staticfloat (thank you so much!) :D

Q.How are we gonna do that?
(I like the way he's seperating topics in)
1. Unpack compiled taballs
1.2. No:
1.2.1. package managers
1.2.2. sudo
1.2.3. Compilation
1.2.4. NO STATE!

2.Build Tarballs
2.1.How: using cross-compilation Environment
2.1.1 Linux based: with:
2.1.1.1. gcc
2.1.1.2. gfortran
2.1.1.3. clang

2.2 Corss compilers with everybody (so that performance gotta be the same )
Note: built for all targets
2.2.1. windows
2.2.2. arm64

2.3 Run inside a Cross Compilation Environment
2.3.1 Buld `all` the things - Anywhere
me: with its own special flavors & Recipies

"""
#=Strategy:
we'll `bake` a little Linux image on any linux or mac machine
Juliacon2018 - Elliot Saba
How: prerequisites: activate current directory and add package BinaryBuilder

1. using BinaryBuilder
2. BinaryBuilder.runshell()

q.

windows -> admin prevelages

either findout chmod and try back on win
    or run in mac

#waiting to be consumed  on mac ...
#1. requires only mac  (Ubuntu references  ) #####requires focal-updates" - Ubuntu Packages # Focal Fossa (Ubuntu) # https://wiki.ubuntu.com/FocalFossa/ReleaseNotes
2. remove
strategy:
either load everything into project toml
=#

    ENV["XDG_DATA_DIRS"] = join(filter(x -> x !== nothing, [
         #   dirname(adwaita_icons_dir), #
         #   dirname(hicolor_icons_dir),
            joinpath(dirname(GTK3_jll.libgdk3_path::String), "..", "share"), #same problem as in Gui.jl
            get(ENV, "XDG_DATA_DIRS", nothing)::Union{String,Nothing},#GTK3_jll not defined
        ]), Sys.iswindows() ? ";" : ":")
#=
=#
#=
TODO: change currrent_pwd() to artifact_hash
# artifact_hash isa notdefined
loaders_cache_hash = artifact_hash(loaders_cache_name, mutable_artifacts_toml)
loaders_dir_hash = artifact_hash(loaders_dir_name, mutable_artifacts_toml)
# Install Ubuntu on WSL2 on Windows 10
# https://aka.ms/wslstore
# install Ubuntu on
https://ubuntu.com/tutorials/install-ubuntu-on-wsl2-on-windows-10#3-download-ubuntu

LoadError: UndefVarError: create_artifact not definedStacktrace:

=#
mutable_artifacts_toml = joinpath(dirname(@__DIR__), "MutableArtifacts.toml")
loaders_cache_name = "gdk-pixbuf-loaders-cache"
#loaders_cache_hash = artifact_hash(loaders_cache_name, mutable_artifacts_toml)
loaders_dir_name = "gdk-pixbuf-loaders-dir"
#loaders_dir_hash = artifact_hash(loaders_dir_name, mutable_artifacts_toml)

# Copy loaders into a directory
loaders_dir_hash = create_artifact() do art_dir
    loaders_dir = mkdir(joinpath(art_dir, "loaders_dir"))
    pixbuf_loaders = joinpath.(gdk_pixbuf_loaders_dir, readdir(gdk_pixbuf_loaders_dir))
    push!(pixbuf_loaders, Librsvg_jll.libpixbufloader_svg)
    cp.(pixbuf_loaders, joinpath.(loaders_dir, basename.(pixbuf_loaders)))
end

loaders_dir = joinpath(artifact_path(loaders_dir_hash), "loaders_dir")
# Pkg removes "execute" permissions on Windows
if Sys.iswindows() && chmod(artifact_path(loaders_dir_hash), 0o755; recursive = true)#grant access to msft

end
using BinaryBuilder
BinaryBuilder.runshell()

max_stars_repo_path: Tutorial/Elliot Saba/BetterBinary.jl

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_path: Tutorial/Elliot Saba/BetterBinary.jl

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_path: Tutorial/Elliot Saba/BetterBinary.jl

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
function calc_napier()
    n_trial =convert(Int,1e8)
    total_count=0
    for i in 1:n_trial
        counter = 0
        accumulated = 0.0
        while true
            accumulated += rand()
            counter += 1
            if accumulated >= 1.0
                break
            end
        end
        total_count+=counter
    end
    total_count/n_trial
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module SerializationCaches

using Serialization
using OrderedCollections

export SerializationCache, fetch!, put!, set_up_cache_path

"""
    SerializationCache{T}(path; in_memory_limit, file_limit, file_gc_ratio)

Return a `SerializationCache` instance that holds elements of type `T`.

A `SerializationCache` (along with [`fetch!`](@ref) and [`put!`](@ref))
implements a simple two-stage cache that is useful for caching objects that take
significantly longer to compute from scratch than to (de)serialize from disk.

Recently fetched objects are cached within an in-memory data structure, while
less recently fetched objects are (de)serialized to/from the filesystem via
Julia's `Serialization` module. Specifically, if adding an object to the
in-memory cache would cause the in-memory cache to grow beyond `in_memory_limit`
elements, the least recently fetched object in the cache is serialized to a
`.jls` file in `path`. Once the number of `.jls` files in `path` exceeds
`file_limit`, less recently deserialized files are deleted to clear space;
the actual ratio of files that are deleted during this "garbage collection"
process is determined by `file_gc_ratio`.

Note that all `.jls` files in `path` at the time of `SerializationCache`
construction are considered to be part of constructed cache.

A valid method for function [`set_up_cache_path(path)`](@ref) must exist for the
given `path`; only `set_up_cache_path(::AbstractString)` currently exists.

See also: [`fetch!`](@ref), [`put!`](@ref), [`set_up_cache_path`](@ref)
"""
struct SerializationCache{T}
    path::Any
    file_keys::OrderedSet{String}
    file_gc_ratio::Float64
    file_limit::Union{Int,Nothing}
    in_memory::OrderedDict{String,T}
    in_memory_limit::Int
    function SerializationCache{T}(path; in_memory_limit, file_limit,
                                   file_gc_ratio) where {T}
        path = set_up_cache_path(path)
        file_keys = OrderedSet{String}(first(splitext(name)) for name in readdir(path)
                                       if endswith(name, ".jls"))
        return new{T}(path, file_keys, file_gc_ratio, file_limit, OrderedDict{String,T}(),
                      in_memory_limit)
    end
end

SerializationCache(args...; kwargs...) = SerializationCache{Any}(args...; kwargs...)

"""
    set_up_cache_path(path::AbstractString)

Return the `SerializationCache`-ready `path`. For input type `AbstractString`,
creates `path` if it is not already a directory. Can be extended to support
other types of (e.g., AWSS3.jl's `S3Path`).
"""
function set_up_cache_path(path::AbstractString)
    path = rstrip(abspath(path), '/')
    isdir(path) || mkpath(path)
    return path
end

"""
    fetch!(f, cache::SerializationCache, key::AbstractString)

Return the object stored at `key` in `cache`. If `key` doesn't exist `cache`,
set `key` to `f()` and return the result.

Note that `key` must be a valid file name.

As part of fetching the requested result, this function performs several
bookkeeping operations to maintain `cache` within its configured limits; see
[`SerializationCache`](@ref) for details.
"""
function fetch!(f, cache::SerializationCache, key::AbstractString)
    if haskey(cache.in_memory, key)
        item = cache.in_memory[key]
        delete!(cache.in_memory, key) # re-prioritize `key`
        cache.in_memory[key] = item
        return item
    end
    if key in cache.file_keys
        item = deserialize(joinpath(cache.path, key * ".jls"))
        delete!(cache.file_keys, key) # re-prioritize `key`
        push!(cache.file_keys, key)
        cache.in_memory_limit == 0 && return item
    else
        item = f()
    end
    return put!(cache, key, item)
end

"""
    put!(cache::SerializationCache, key::AbstractString, item;
         directly_to_file::Bool=false)

Store `item` in `cache` at `key` and return `item`.

If `directly_to_file == true`, `item` is directly serialized to the `cache`'s
filesystem layer, skipping the `cache`'s in-memory layer.

Note that `key` must be a valid file name.
"""
function Base.put!(cache::SerializationCache, key::AbstractString, item;
                   directly_to_file::Bool=false)
    if directly_to_file || cache.in_memory_limit == 0
        _file_layer_put!(cache, key, item)
    else
        haskey(cache.in_memory, key) && delete!(cache.in_memory, key) # re-prioritize 'key'
        cache.in_memory[key] = item
        if length(cache.in_memory) > cache.in_memory_limit
            # too bad `popfirst!` isn't implemented on `OrderedDict`...
            oldest_key = first(keys(cache.in_memory))
            oldest_item = cache.in_memory[oldest_key]
            delete!(cache.in_memory, oldest_key)
            _file_layer_put!(cache, oldest_key, oldest_item)
        end
    end
    return item
end

function _popfirst!(ordered_set::OrderedSet)
    item = first(ordered_set)
    delete!(ordered_set, item)
    return item
end

function _file_layer_put!(cache::SerializationCache, key, item)
    serialize(joinpath(cache.path, key * ".jls"), item)
    if key in cache.file_keys # re-prioritize 'key'
        delete!(cache.file_keys, key)
    end
    push!(cache.file_keys, key)
    file_count = length(cache.file_keys)
    if !isnothing(cache.file_limit) && file_count > cache.file_limit
        delete_count = cache.file_gc_ratio * file_count
        for _ in 1:clamp(floor(Int, delete_count), 1, file_count)
            rm(joinpath(cache.path, _popfirst!(cache.file_keys) * ".jls"))
        end
    end
end

end # module


text:
"""
    @syncthrow ex

Equivalent to just `ex` but let `@sync` include the exception thrown from
`ex` to `CompositeException`.
"""
macro syncthrow(ex)
    @gensym err
    quote
        $Base.@isdefined($(Base.sync_varname)) ||
            $error("`@syncthrow` requires outer `@sync`")
        try
            $ex
        catch $err
            # $Base.@debug "@syncthrow" exception = ($err, catch_backtrace())
            $Base.@sync_add $Thrower($err, $catch_backtrace())
        end
    end |> esc
end

const BacktraceType = try
    error("")
catch
    catch_backtrace()
end |> typeof

struct Thrower <: Exception
    err::Any
    backtrace::BacktraceType
end

Base.wait(t::Thrower) = throw(t)

function Base.showerror(io::IO, t::Thrower)
    print(io, "Thrower: ")
    showerror(io, t.err, t.backtrace)
end

@enum RefState begin
    REF_UNSET
    REF_SET
    REF_CLOSED
end

struct LockedRef{T}
    ref::Base.RefValue{T}
    state::Base.RefValue{RefState}
    cond::Threads.Condition
end

LockedRef{T}() where {T} = LockedRef{T}(Ref{T}(), Ref(REF_UNSET), Threads.Condition())
LockedRef(x) = LockedRef(Ref(x), Ref(REF_SET), Threads.Condition())

function tryfetch(ref::LockedRef)
    lock(ref.cond) do
        while true
            state = ref.state[]
            state == REF_SET && return Some(ref.ref[])
            state == REF_CLOSED && return nothing
            wait(ref.cond)
        end
    end
end

#=
function Base.getindex(ref::LockedRef)
    value = tryfetch(ref)
    value === nothing && error("LockedRef is closed")
    return something(value)
end
=#

function tryset!(ref::LockedRef, value)
    lock(ref.cond) do
        state = ref.state[]
        state == REF_CLOSED && return false
        ref.state[] = REF_SET
        ref.ref[] = value
        state == REF_UNSET && notify(ref.cond)
        return true
    end
end

#=
Base.setindex!(ref::LockedRef, value) =
    tryset!(ref, value) || error("LockedRef is closed")
=#

function Base.close(ref::LockedRef)
    lock(ref.cond) do
        ref.state[] = REF_CLOSED
        notify(ref.cond)
    end
end


text:
function mpccmodel_pointeval_basic(model::MPCCModel, ptevalreq::MPCCPointEvalReq, x0::Vector{S}, pr0::Vector{T}, ps0::Vector{Int64}) where {S <: Real, T <: Real}

    # std entries
    if ( ptevalreq.f )
        f = model.f(x0, pr0, ps0)
    else
        f = missing
    end

    if ( ptevalreq.ce )
        ce = model.ce(x0, pr0, ps0)
    else
        ce = missing
    end

    if ( ptevalreq.ci )
        ci = model.ci(x0, pr0, ps0)
    else
        ci = missing
    end

    if ( ptevalreq.F )
        F = model.F(x0, pr0, ps0)
    else
        F = missing
    end

    # grad entries
    if ( ptevalreq.gradf )
        gradf = model.gradf(x0, pr0, ps0)
    else
        gradf = missing
    end

    if ( ptevalreq.gradce )
        gradce = model.gradce(x0, pr0, ps0)
    else
        gradce = missing
    end

    if ( ptevalreq.gradci )
        gradci = model.gradci(x0, pr0, ps0)
    else
        gradci = missing
    end

    if ( ptevalreq.gradF )
        gradF = model.gradF(x0, pr0, ps0)
    else
        gradF = missing
    end

    # hess entries
    if ( ptevalreq.hessf )
        hessf = model.hessf(x0, pr0, ps0)
    else
        hessf = missing
    end

    if ( ptevalreq.hessce )
        hessce = model.hessce(x0, pr0, ps0)
    else
        hessce = missing
    end

    if ( ptevalreq.hessci )
        hessci = model.hessci(x0, pr0, ps0)
    else
        hessci = missing
    end

    if ( ptevalreq.hessF )
        hessF = model.hessF(x0, pr0, ps0)
    else
        hessF = missing
    end

    # dp entries
    if ( ptevalreq.fdp )
        fdp = model.fdp(x0, pr0, ps0)
    else
        fdp = missing
    end

    if ( ptevalreq.cedp )
        cedp = model.cedp(x0, pr0, ps0)
    else
        cedp = missing
    end

    if ( ptevalreq.cidp )
        cidp = model.cidp(x0, pr0, ps0)
    else
        cidp = missing
    end

    if ( ptevalreq.Fdp )
        Fdp = model.Fdp(x0, pr0, ps0)
    else
        Fdp = missing
    end

    # graddp entries
    if ( ptevalreq.gradfdp )
        gradfdp = model.gradfdp(x0, pr0, ps0)
    else
        gradfdp = missing
    end

    if ( ptevalreq.gradcedp )
        gradcedp = model.gradcedp(x0, pr0, ps0)
    else
        gradcedp = missing
    end

    if ( ptevalreq.gradcidp )
        gradcidp = model.gradcidp(x0, pr0, ps0)
    else
        gradcidp = missing
    end

    if ( ptevalreq.gradFdp )
        gradFdp = model.gradFdp(x0, pr0, ps0)
    else
        gradFdp = missing
    end

    return MPCCPointEval{promote_type(S, T)}(
        f, ce, ci, F,
        gradf, gradce, gradci, gradF,
        hessf, hessce, hessci, hessF,
        fdp, cedp, cidp, Fdp,
        gradfdp, gradcedp, gradcidp, gradFdp
    )

end


text:
EmptyWorkbook() = Workbook(EmptyMSOfficePackage(), Vector{Worksheet}(), false,
    Vector{Relationship}(), SharedStringTable(), Dict{Int, Bool}(), Dict{Int, Bool}(),
    Dict{String, DefinedNameValueTypes}(), Dict{Tuple{Int, String}, DefinedNameValueTypes}(), nothing)

"""
    is_writable(xl::XLSXFile)

Indicates wether this XLSX file can be edited.
This controls if assignment to worksheet cells is allowed.
Writable XLSXFile instances are opened with `XLSX.open_xlsx_template` method.
"""
is_writable(xl::XLSXFile) = xl.is_writable

"""
Lists Worksheet names for this Workbook.
"""
sheetnames(wb::Workbook) = [ s.name for s in wb.sheets ]
@inline sheetnames(xl::XLSXFile) = sheetnames(xl.workbook)

function hassheet(wb::Workbook, sheetname::AbstractString) :: Bool
    for s in wb.sheets
        if s.name == sheetname
            return true
        end
    end
    return false
end

@inline hassheet(xl::XLSXFile, sheetname::AbstractString) = hassheet(xl.workbook, sheetname)

"""
Counts the number of sheets in the Workbook.
"""
@inline sheetcount(wb::Workbook) = length(wb.sheets)
@inline sheetcount(xl::XLSXFile) = sheetcount(xl.workbook)

"""
    isdate1904(wb) :: Bool

Returns true if workbook follows date1904 convention.
"""
@inline isdate1904(wb::Workbook) :: Bool = wb.date1904
@inline isdate1904(xf::XLSXFile) :: Bool = isdate1904(get_workbook(xf))

function getsheet(wb::Workbook, sheetname::String) :: Worksheet
    for ws in wb.sheets
        if ws.name == sheetname
            return ws
        end
    end
    error("$(get_xlsxfile(wb).filepath) does not have a Worksheet named $sheetname.")
end

@inline getsheet(wb::Workbook, sheet_index::Int) :: Worksheet = wb.sheets[sheet_index]
@inline getsheet(xl::XLSXFile, sheetname::String) :: Worksheet = getsheet(xl.workbook, sheetname)
@inline getsheet(xl::XLSXFile, sheet_index::Int) :: Worksheet = getsheet(xl.workbook, sheet_index)

function Base.show(io::IO, xf::XLSXFile)
    wb = xf.workbook
    print(io, "XLSXFile(\"$(basename(xf.filepath))\") ",
              "containing $(sheetcount(wb)) Worksheets\n")
    @printf(io, "%21s %-13s %-13s\n", "sheetname", "size", "range")
    println(io, "-"^(21+1+13+1+13))

    for s in wb.sheets
        name = s.name |> x -> length(x) > 20 ? x[1:20]*"…" : x
        if s.dimension != nothing
            rg = s.dimension
            _size = size(rg) |> x -> string(x[1], "x", x[2])
            @printf(io, "%21s %-13s %-13s\n", name, _size, rg)
        else
            @printf(io, "%21s size unknown\n", name)
        end
    end
end

@inline Base.getindex(xl::XLSXFile, i::Integer) = getsheet(xl, i)

function Base.getindex(xl::XLSXFile, s::AbstractString)
    if hassheet(xl, s)
        return getsheet(xl, s)
    else
        return getdata(xl, s)
    end
end

function getdata(xl::XLSXFile, ref::SheetCellRef)
    @assert hassheet(xl, ref.sheet) "Sheet $(ref.sheet) not found."
    return getdata(getsheet(xl, ref.sheet), ref.cellref)
end

function getdata(xl::XLSXFile, rng::SheetCellRange)
    @assert hassheet(xl, rng.sheet) "Sheet $(rng.sheet) not found."
    return getdata(getsheet(xl, rng.sheet), rng.rng)
end

function getdata(xl::XLSXFile, rng::SheetColumnRange)
    @assert hassheet(xl, rng.sheet) "Sheet $(rng.sheet) not found."
    return getdata(getsheet(xl, rng.sheet), rng.colrng)
end

function getdata(xl::XLSXFile, s::AbstractString)
    if is_valid_sheet_cellname(s)
        return getdata(xl, SheetCellRef(s))
    elseif is_valid_sheet_cellrange(s)
        return getdata(xl, SheetCellRange(s))
    elseif is_valid_sheet_column_range(s)
        return getdata(xl, SheetColumnRange(s))
    elseif is_workbook_defined_name(xl, s)
        v = get_defined_name_value(xl.workbook, s)
        if is_defined_name_value_a_constant(v)
            return v
        elseif is_defined_name_value_a_reference(v)
            return getdata(xl, v)
        else
            error("Unexpected defined name value: $v.")
        end
    end

    error("$s is not a valid sheetname or cell/range reference.")
end

function getcell(xl::XLSXFile, ref::SheetCellRef)
    @assert hassheet(xl, ref.sheet) "Sheet $(ref.sheet) not found."
    return getcell(getsheet(xl, ref.sheet), ref.cellref)
end

getcell(xl::XLSXFile, ref_str::AbstractString) = getcell(xl, SheetCellRef(ref_str))

function getcellrange(xl::XLSXFile, rng::SheetCellRange)
    @assert hassheet(xl, rng.sheet) "Sheet $(rng.sheet) not found."
    return getcellrange(getsheet(xl, rng.sheet), rng.rng)
end

function getcellrange(xl::XLSXFile, rng::SheetColumnRange)
    @assert hassheet(xl, rng.sheet) "Sheet $(rng.sheet) not found."
    return getcellrange(getsheet(xl, rng.sheet), rng.colrng)
end

function getcellrange(xl::XLSXFile, rng_str::AbstractString)
    if is_valid_sheet_cellrange(rng_str)
        return getcellrange(xl, SheetCellRange(rng_str))
    elseif is_valid_sheet_column_range(rng_str)
        return getcellrange(xl, SheetColumnRange(rng_str))
    end

    error("$rng_str is not a valid range reference.")
end

@inline is_workbook_defined_name(wb::Workbook, name::AbstractString) :: Bool = haskey(wb.workbook_names, name)
@inline is_workbook_defined_name(xl::XLSXFile, name::AbstractString) :: Bool = is_workbook_defined_name(get_workbook(xl), name)
@inline is_worksheet_defined_name(ws::Worksheet, name::AbstractString) :: Bool = is_worksheet_defined_name(get_workbook(ws), ws.sheetId, name)
@inline is_worksheet_defined_name(wb::Workbook, sheetId::Int, name::AbstractString) :: Bool = haskey(wb.worksheet_names, (sheetId, name))
@inline is_worksheet_defined_name(wb::Workbook, sheet_name::AbstractString, name::AbstractString) :: Bool = is_worksheet_defined_name(wb, getsheet(wb, sheet_name).sheetId, name)

@inline get_defined_name_value(wb::Workbook, name::AbstractString) :: DefinedNameValueTypes = wb.workbook_names[name]

function get_defined_name_value(ws::Worksheet, name::AbstractString) :: DefinedNameValueTypes
    wb = get_workbook(ws)
    sheetId = ws.sheetId
    return wb.worksheet_names[(sheetId, name)]
end

@inline is_defined_name_value_a_reference(v::DefinedNameValueTypes) = isa(v, SheetCellRef) || isa(v, SheetCellRange)
@inline is_defined_name_value_a_constant(v::DefinedNameValueTypes) = !is_defined_name_value_a_reference(v)


text:
#Basic structure for Unmarked data objects
struct UmData
  y::Array{Union{Int,Missing}}
  site_covs::DataFrame
  obs_covs::DataFrame
  UmData(y,site_covs,obs_covs) = new(y, check_data(site_covs, size(y)[1]),
                                     check_data(obs_covs, prod(size(y))))
end

#Check input data frames for validity
function check_data(dat::Union{DataFrame,Nothing}, exp_rows::Int)
  if isnothing(dat) return DataFrame(_dummy=ones(exp_rows)) end
  if any([any(ismissing.(col)) for col = eachcol(dat)])
    error("Covariate DataFrames cannot contain missing values")
  end
  if nrow(dat) != exp_rows
    error("Supplied DataFrame has incorrect number of rows")
  end
  return dat
end

"Construct Unmarked data objects"
function UmData(y::Union{Array{Int}, Array{Union{Int,Missing}}};
                site_covs::Union{DataFrame,Nothing}=nothing,
                obs_covs::Union{DataFrame,Nothing}=nothing)
  N,J = size(y)
  return UmData(y, check_data(site_covs, N), check_data(obs_covs, N*J))
end

#Subset a UmData object
function _subset_obs(data::UmData, i::Union{Array{Int},UnitRange})
  N,J = size(data.y)
  site_idx = repeat(1:N, inner=J)
  keep = in.(site_idx, (i,))
  return keep
end

function Base.getindex(data::UmData, i::Union{Int,Array{Int},UnitRange})
  if typeof(i) == Int i = i:i end
  new_sc = data.site_covs[i,:]
  new_y = data.y[i,:]
  new_oc = data.obs_covs[_subset_obs(data, i),:]
  UmData(new_y,new_sc,new_oc)
end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
export NullPredictor

"""
    NullPredictor()

A predictor which does no prediction step, i.e., it just returns the input as
its prediction.
"""
struct NullPredictor <: AbstractPredictor end
struct NullPredictorCache <: AbstractPredictorCache end

cache(::NullPredictor, H, x, ẋ, t) = NullPredictorCache()

function predict!(xnext, ::NullPredictor, ::NullPredictorCache, H, x, t, dt, ẋ, Jac)
    xnext .= x
    nothing
end

order(::NullPredictor) = 1


text:
"""
    Free()

Zero external potential everywhere. Useful for modelling free particles.

```jldoctest model
julia> model, R = Free(3), rand(3, 10);

julia> potential(model, R)
0.0

julia> derivative(model, R)
3×10 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
struct Free <: AdiabaticModel
    dofs::Int
end

Free() = Free(1)

NQCModels.ndofs(free::Free) = free.dofs

NQCModels.potential(::Free, ::AbstractMatrix{T}) where {T} = zero(T)
function NQCModels.derivative!(::Free, out::AbstractMatrix, ::AbstractMatrix{T}) where {T}
    out .= 0
    copyto!(out, zero(T))
end


text:
using RbfInterpolationTools
using Documenter

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

makedocs(;
    modules = [RbfInterpolationTools],
    authors = "munozariasjm <munozariasjm@hotmail.com> and contributors",
    repo = "https://github.com/munozariasjm/RbfInterpolationTools.jl",
    sitename = "RbfInterpolationTools.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://munozariasjm.github.io/RbfInterpolationTools.jl",
        assets = String[]
    ),
    pages = [
        "Home" => "index.md",
        "Examples" => "examples.md",
        "Available Kernels" => "kernels.md",
        "Useful References" => "references.md"
    ]
)

deploydocs(;
    repo = "github.com/munozariasjm/RbfInterpolationTools.jl",
    devbranch = "main"
)


text:
export Bloating

"""
    Bloating{N, S<:LazySet{N}} <:LazySet{N}

Type that represents a uniform expansion of a set in a given norm (also
known as *bloating*).

### Fields

- `X` -- set
- `ε` -- (positive) bloating factor
- `p` -- ``p``-norm (should be ``≥ 1``; default: ``2``)

### Notes

If `ε` is positive, then `Bloating(X, ε, p)` is equivalent to the Minkowski sum
of `X` and a ball in the `p`-norm of radius `ε` centered in the origin `O`
(i.e., `X ⊕ Ballp(p, O, ε)`).

The `Bloating` operation preserves convexity: if `X` is convex, then any
bloating of `X` is convex as well.
"""
struct Bloating{N, S<:LazySet{N}} <:LazySet{N}
    X::S
    ε::N
    p::N

    function Bloating(X::S, ε::N, p::N=N(2)) where {N, S<:LazySet{N}}
        @assert p >= one(N) "bloating requires a norm >= 1, but $p was given"

        return new{N, S}(X, ε, p)
    end
end

isoperationtype(::Type{<:Bloating}) = true
isconvextype(::Type{Bloating{N, S}}) where {N, S} = isconvextype(S)

"""
    dim(B::Bloating)

Return the dimension of a bloated set.

### Input

- `B` -- bloated set

### Output

The ambient dimension of the bloated set.
"""
function dim(B::Bloating)
    return dim(B.X)
end

# helper function to compute the bloating ball
function _bloating_ball(B::Bloating{N}) where {N}
    return Ballp(B.p, zeros(N, dim(B)), abs(B.ε))
end

"""
    σ(d::AbstractVector, B::Bloating)

Return the support vector of a bloated set in a given direction.

### Input

- `d` -- direction
- `B` -- bloated set

### Output

The support vector of the bloated set in the given direction.
"""
function σ(d::AbstractVector, B::Bloating)
    @assert !iszero(d) "the support vector in the zero direction is undefined"
    @assert B.ε >= 0 || B.p > 1 "the support vector for negative bloating " *
        "in the 1-norm is not implemented"

    return σ(d, B.X) + sign_cadlag(B.ε) * σ(d, _bloating_ball(B))
end

"""
    ρ(d::AbstractVector, B::Bloating)

Return the support function of a bloated set in a given direction.

### Input

- `d` -- direction
- `B` -- bloated set

### Output

The support function of the bloated set in the given direction.
"""
function ρ(d::AbstractVector, B::Bloating)
    @assert !iszero(d) "the support function in the zero direction is undefined"

    return ρ(d, B.X) + sign_cadlag(B.ε) * ρ(d, _bloating_ball(B))
end

"""
    isbounded(B::Bloating)

Determine whether a bloated set is bounded.

### Input

- `B` -- bloated set

### Output

`true` iff the wrapped set is bounded.
"""
function isbounded(B::Bloating)
    return isbounded(B.X)
end

function isboundedtype(::Type{<:Bloating{N, S}}) where {N, S}
    return isboundedtype(S)
end

"""
    isempty(B::Bloating)

Determine whether a bloated set is empty.

### Input

- `B` -- bloated set

### Output

`true` iff the wrapped set is empty.
"""
function isempty(B::Bloating)
    return isempty(B.X)
end

"""
    an_element(B::Bloating)

Return some element of a bloated set.

### Input

- `B` -- bloated set

### Output

An element in the bloated set.

### Algorithm

The implementation returns the result of `an_element` for the wrapped set.
"""
function an_element(B::Bloating)
    return an_element(B.X)
end

"""
    constraints_list(B::Bloating)

Return the list of constraints of a bloated set.

### Input

- `B` -- bloated set

### Output

The list of constraints of the bloated set.

### Notes

The constraints list is only available for bloating in the `p`-norm for
``p = 1`` or ``p = ∞`` and if `constraints_list` is available for the unbloated
set.

### Algorithm

We call `constraints_list` on the lazy Minkowski sum.
"""
function constraints_list(B::Bloating)
    @assert (B.p == 1 || B.p == Inf) "the constraints list is only available " *
        "for bloating in the 1-norm or in the infinity norm"
    if B.ε < 0
        throw(ArgumentError("computing the constraints list of a negatively " *
                            "bloated set is not supported yet"))
    end

    return constraints_list(MinkowskiSum(B.X, _bloating_ball(B)))
end


text:
# Reproduce Kaplan Violante 2010
# Allen Peters
# November 7, 2019

# t ∈ [1,T]
# vectors n x T

using Plots, LinearAlgebra, Statistics, Random
using BenchmarkTools

# Income process
function sim_income1(N, T_ret, σ_ε, σ_η, σ_z0, κ; y_0 = 0.0)
    Y = zeros(eltype(y_0), n, T_ret)
    ε = sqrt(σ_ε) .* randn(n, T_ret)
    η = sqrt(σ_η) .* randn(n, T_ret)
    z_t1 = sqrt(σ_z0) .* randn(n)

    for t in 1:(T_ret-1)
        z_t = z_t1 + η[:,t]
        Y[:,t] = exp.( κ[t] .+ z_t + ε[:,t] )
        z_t1 = copy(z_t)
    end
    return Y
end

function sim_income2(N, T_ret, σ_ε, σ_η, σ_z0, κ; y_0 = 0.0)
    Y = zeros(eltype(y_0), n, T_ret)
    z_t1 = sqrt(σ_z0) .* randn(n)

    for t in 1:(T_ret-1)
        z_t = z_t1 + sqrt(σ_η) .* randn(n)
        Y[:,t] = exp.( κ[t] .+ z_t + sqrt(σ_ε) .* randn(n) )
        z_t1 = copy(z_t)
    end
    return Y
end

function sim_income3(N, T_ret, σ_ε, σ_η, σ_z0, κ; y_0 = 0.0)
    Y = zeros(eltype(y_0), n, T_ret)
    z_t1 = sqrt(σ_z0) .* randn(n)
    z_t = similar(z_t1)
    for t in 1:(T_ret-1), i in 1:N
        z_t[i] = z_t1[i] + sqrt(σ_η) * randn()
        Y[i,t] = exp.( κ[t] + z_t[i] + sqrt(σ_ε) * randn() )
        z_t1[i] = copy(z_t[i])
    end
    return Y
end

N = 50000
T_ret = 40
T = 70
y_0 = 0.0
κ = zeros(eltype(y_0),N,T)
σ_ε = 1
σ_η = 1
σ_z0 = 1
t = 1:T

@btime Y = sim_income1(N, T_ret, σ_ε, σ_η, σ_z0, κ; y_0 = 0.0)
@btime Y = sim_income2(N, T_ret, σ_ε, σ_η, σ_z0, κ; y_0 = 0.0)
@btime Y = sim_income2flip(N, T_ret, σ_ε, σ_η, σ_z0, κ; y_0 = 0.0)
@time Y = sim_income3(N, T_ret, σ_ε, σ_η, σ_z0, κ; y_0 = 0.0)
@time Y = sim_income3flip(N, T_ret, σ_ε, σ_η, σ_z0, κ; y_0 = 0.0)

Random.seed!(1234)
Y = sim_income2(100, T_ret, σ_ε, σ_η, σ_z0, κ; y_0 = 0.0)

plot!(Y[6,:])

plot()

function sim_income2flip(N, T_ret, σ_ε, σ_η, σ_z0, κ; y_0 = 0.0)
    Y = zeros(eltype(y_0), T_ret,n)
    z_t1 = sqrt(σ_z0) .* randn(n)

    for t in 1:(T_ret-1)
        z_t = z_t1 + sqrt(σ_η) .* randn(n)
        Y[t,:] = exp.( κ[t] .+ z_t + sqrt(σ_ε) .* randn(n) )
        z_t1 = copy(z_t)
    end
    return Y
end

function sim_income3flip(N, T_ret, σ_ε, σ_η, σ_z0, κ; y_0 = 0.0)
    Y = zeros(eltype(y_0), T_ret,n)
    z_t1 = sqrt(σ_z0) .* randn(n)
    z_t = similar(z_t1)
    for i in 1:N, t in 1:(T_ret-1)
        z_t[i] = z_t1[i] + sqrt(σ_η) * randn()
        Y[t,i] = exp.( κ[t] + z_t[i] + sqrt(σ_ε) * randn() )
        z_t1[i] = copy(z_t[i])
    end
    return Y
end

max_stars_repo_path: test/Income tests.jl

max_issues_repo_path: test/Income tests.jl

max_forks_repo_path: test/Income tests.jl


text:
function center!( A )
    colmeans = mean( A, 1 )
    broadcast!( -, A, A, colmeans )
    return colmeans
end

function center( A )
    return broadcast( -, A, mean( A, 1 ) )
end

max_stars_repo_name: tbepler/julia-regression

max_issues_repo_name: tbepler/julia-regression

max_forks_repo_name: tbepler/julia-regression


text:
@test @inferred(ArrayInterface.size(sv5)) === (StaticInt(5),)
@test @inferred(ArrayInterface.size(v5)) === (5,)
@test @inferred(ArrayInterface.size(A)) === (3, 4, 5)
@test @inferred(ArrayInterface.size(Ap)) === (2, 5)
@test @inferred(ArrayInterface.size(A)) === size(A)
@test @inferred(ArrayInterface.size(Ap)) === size(Ap)
@test @inferred(ArrayInterface.size(R)) === (StaticInt(2),)
@test @inferred(ArrayInterface.size(Rnr)) === (StaticInt(4),)
@test @inferred(ArrayInterface.known_length(Rnr)) === 4
@test @inferred(ArrayInterface.size(A2)) === (4, 3, 5)
@test @inferred(ArrayInterface.size(A2r)) === (2, 3, 5)

@test @inferred(ArrayInterface.size(view(rand(4), reshape(1:4, 2, 2)))) == (2, 2)
@test @inferred(ArrayInterface.size(irev)) === (StaticInt(2), StaticInt(3), StaticInt(4))
@test @inferred(ArrayInterface.size(iprod)) === (StaticInt(2), StaticInt(3), StaticInt(4))
@test @inferred(ArrayInterface.size(iflat)) === (static(72),)
@test @inferred(ArrayInterface.size(igen)) === (StaticInt(2), StaticInt(3), StaticInt(4))
@test @inferred(ArrayInterface.size(iacc)) === (StaticInt(2), StaticInt(3), StaticInt(4))
@test @inferred(ArrayInterface.size(ienum)) === (StaticInt(2), StaticInt(3), StaticInt(4))
@test @inferred(ArrayInterface.size(ipairs)) === (StaticInt(2), StaticInt(3), StaticInt(4))
@test @inferred(ArrayInterface.size(izip)) === (StaticInt(2), StaticInt(3), StaticInt(4))
@test @inferred(ArrayInterface.size(zip(S, A_trailingdim))) === (StaticInt(2), StaticInt(3), StaticInt(4), static(1))
@test @inferred(ArrayInterface.size(zip(A_trailingdim, S))) === (StaticInt(2), StaticInt(3), StaticInt(4), static(1))
@test @inferred(ArrayInterface.size(S)) === (StaticInt(2), StaticInt(3), StaticInt(4))
@test @inferred(ArrayInterface.size(Sp)) === (2, 2, StaticInt(3))
@test @inferred(ArrayInterface.size(Sp2)) === (2, StaticInt(3), StaticInt(2))
@test @inferred(ArrayInterface.size(S)) == size(S)
@test @inferred(ArrayInterface.size(Sp)) == size(Sp)
@test @inferred(ArrayInterface.size(parent(Sp2))) === (static(4), static(3), static(2))
@test @inferred(ArrayInterface.size(Sp2)) == size(Sp2)
@test @inferred(ArrayInterface.size(Sp2, StaticInt(1))) === 2
@test @inferred(ArrayInterface.size(Sp2, StaticInt(2))) === StaticInt(3)
@test @inferred(ArrayInterface.size(Sp2, StaticInt(3))) === StaticInt(2)
@test @inferred(ArrayInterface.size(Wrapper(Sp2), StaticInt(3))) === StaticInt(2)
@test @inferred(ArrayInterface.size(Diagonal([1, 2]))) == size(Diagonal([1, 2]))
@test @inferred(ArrayInterface.size(Mp)) === (StaticInt(3), StaticInt(4))
@test @inferred(ArrayInterface.size(Mp2)) === (StaticInt(2), 2)
@test @inferred(ArrayInterface.size(Mp)) == size(Mp)
@test @inferred(ArrayInterface.size(Mp2)) == size(Mp2)

@test @inferred(ArrayInterface.known_size(1)) === ()
@test @inferred(ArrayInterface.known_size(view(rand(4), reshape(1:4, 2, 2)))) == (nothing, nothing)
@test @inferred(ArrayInterface.known_size(A)) === (nothing, nothing, nothing)
@test @inferred(ArrayInterface.known_size(Ap)) === (nothing, nothing)
@test @inferred(ArrayInterface.known_size(Wrapper(Ap))) === (nothing, nothing)
@test @inferred(ArrayInterface.known_size(R)) === (2,)
@test @inferred(ArrayInterface.known_size(Wrapper(R))) === (2,)
@test @inferred(ArrayInterface.known_size(Rnr)) === (4,)
@test @inferred(ArrayInterface.known_size(Rnr, static(1))) === 4
@test @inferred(ArrayInterface.known_size(Ar)) === (nothing, nothing, nothing,)
@test @inferred(ArrayInterface.known_size(Ar, static(1))) === nothing
@test @inferred(ArrayInterface.known_size(Ar, static(4))) === 1
@test @inferred(ArrayInterface.known_size(A2)) === (nothing, nothing, nothing)
@test @inferred(ArrayInterface.known_size(A2r)) === (nothing, nothing, nothing)

@test @inferred(ArrayInterface.known_size(irev)) === (2, 3, 4)
@test @inferred(ArrayInterface.known_size(igen)) === (2, 3, 4)
@test @inferred(ArrayInterface.known_size(iprod)) === (2, 3, 4)
@test @inferred(ArrayInterface.known_size(iflat)) === (72,)
@test @inferred(ArrayInterface.known_size(iacc)) === (2, 3, 4)
@test @inferred(ArrayInterface.known_size(ienum)) === (2, 3, 4)
@test @inferred(ArrayInterface.known_size(izip)) === (2, 3, 4)
@test @inferred(ArrayInterface.known_size(ipairs)) === (2, 3, 4)
@test @inferred(ArrayInterface.known_size(zip(S, A_trailingdim))) === (2, 3, 4, 1)
@test @inferred(ArrayInterface.known_size(zip(A_trailingdim, S))) === (2, 3, 4, 1)
@test @inferred(ArrayInterface.known_length(Iterators.flatten(((x, y) for x in 0:1 for y in 'a':'c')))) === nothing

@test @inferred(ArrayInterface.known_size(S)) === (2, 3, 4)
@test @inferred(ArrayInterface.known_size(Wrapper(S))) === (2, 3, 4)
@test @inferred(ArrayInterface.known_size(Sp)) === (nothing, nothing, 3)
@test @inferred(ArrayInterface.known_size(Wrapper(Sp))) === (nothing, nothing, 3)
@test @inferred(ArrayInterface.known_size(Sp2)) === (nothing, 3, 2)
@test @inferred(ArrayInterface.known_size(Sp2, StaticInt(1))) === nothing
@test @inferred(ArrayInterface.known_size(Sp2, StaticInt(2))) === 3
@test @inferred(ArrayInterface.known_size(Sp2, StaticInt(3))) === 2
@test @inferred(ArrayInterface.known_size(Mp)) === (3, 4)
@test @inferred(ArrayInterface.known_size(Mp2)) === (2, nothing)

@testset "known_length" begin
    @test ArrayInterface.known_length(1:2) === nothing
    @test ArrayInterface.known_length((1,)) == 1
    @test ArrayInterface.known_length((a=1, b=2)) == 2
    @test ArrayInterface.known_length([]) === nothing
    @test ArrayInterface.known_length(CartesianIndex((1, 2, 3))) === 3
    @test @inferred(ArrayInterface.known_length(NDIndex((1, 2, 3)))) === 3

    itr = StaticInt(1):StaticInt(10)
    @inferred(ArrayInterface.known_length((i for i in itr))) == 10
end

@testset "is_lazy_conjugate" begin
    a = rand(ComplexF64, 2)
    @test @inferred(ArrayInterface.is_lazy_conjugate(a)) == false
    b = a'
    @test @inferred(ArrayInterface.is_lazy_conjugate(b)) == true
    c = transpose(b)
    @test @inferred(ArrayInterface.is_lazy_conjugate(c)) == true
    d = c'
    @test @inferred(ArrayInterface.is_lazy_conjugate(d)) == false
    e = permutedims(d)
    @test @inferred(ArrayInterface.is_lazy_conjugate(e)) == false
    @test @inferred(ArrayInterface.is_lazy_conjugate([1, 2, 3]')) == false # We don't care about conj on `<:Real`
end


text:
"""
    orbital_priority(fun, orig_cfg, orbitals)

Generate priorities for the substitution `orbitals`, i.e. the
preferred ordering of the orbitals in configurations excited from
`orig_cfg`. `fun` can optionally transform the labels of substitution
orbitals, in which case they will be ordered just after their parent
orbital in the source configuration; otherwise they will be appended
to the priority list.
"""
function orbital_priority(fun::Function, orig_cfg::Configuration{O₁}, orbitals::Vector{O₂}) where {O₁<:AbstractOrbital,O₂<:AbstractOrbital}
    orbital_priority = Vector{promote_type(O₁,O₂)}()
    append!(orbital_priority, orig_cfg.orbitals)
    for (orb,occ,state) in orig_cfg
        for new_orb in orbitals
            subs_orb = fun(new_orb, orb)
            (new_orb == subs_orb || subs_orb ∈ orbital_priority) && continue
            push!(orbital_priority, subs_orb)
        end
    end
    # This is not exactly correct, since if there is a transformation
    # in effect (`fun` not returning the same orbital), they should
    # not be included in the orbital priority list, whereas if there
    # is no transformation, they should be appended since they have
    # the lowest priority (compared to the original
    # orbitals). However, if there is a transformation, the resulting
    # orbitals will be different from the ingoing ones, and the
    # subsequent excitations will never consider the untransformed
    # substitution orbitals, so there is really no harm done.
    for new_orb in orbitals
        (new_orb ∈ orbital_priority) && continue
        push!(orbital_priority, new_orb)
    end

    Dict(o => i for (i,o) in enumerate(orbital_priority))
end

function single_excitations!(fun::Function,
                             excitations::Vector{<:Configuration},
                             ref_set::Configuration,
                             orbitals::Vector{<:AbstractOrbital},
                             min_occupancy::Vector{Int},
                             max_occupancy::Vector{Int},
                             excite_from::Int)
    orig_cfg = first(excitations)
    orbital_order = orbital_priority(fun, orig_cfg, orbitals)
    for config ∈ excitations[end-excite_from+1:end]
        for (orb_i,(orb,occ,state)) ∈ enumerate(config)
            state != :open && continue
            # If the orbital we propose to excite from is among those
            # from the reference set and already at its minimum
            # occupancy, we continue.
            i = findfirst(isequal(orb), ref_set.orbitals)
            !isnothing(i) && occ == min_occupancy[i] && continue

            for (new_orb_i,new_orb) ∈ enumerate(orbitals)
                subs_orb = fun(new_orb, orb)
                subs_orb == orb && continue
                # If the proposed substitution orbital is among those
                # from the reference set and already present in the
                # configuration to excite, we check if it is already
                # at its maximum occupancy, in which case we continue.
                j = findfirst(isequal(subs_orb), ref_set.orbitals)
                k = findfirst(isequal(subs_orb), config.orbitals)
                !isnothing(j) && !isnothing(k) && config.occupancy[k] == max_occupancy[j] && continue
                # Likewise, if the proposed substitution orbital is
                # already at its maximum, due to the degeneracy, we
                # continue.
                !isnothing(k) && config.occupancy[k] == degeneracy(subs_orb) && continue

                # If we are working with unsorted configurations and
                # the substitution orbital is not already present in
                # config, we need to check if the substitution would
                # generate a configuration that violates the desired
                # orbital priority.
                if isnothing(k) && !config.sorted
                    sp = orbital_order[subs_orb]

                    # Make sure that no preceding orbital has higher
                    # value than subs_orb.
                    any(orbital_order[o] > sp
                        for o in config.orbitals) && continue
                end

                excited_config = replace(config, orb=>subs_orb, append=true)
                excited_config ∉ excitations && push!(excitations, excited_config)
            end
        end
    end
end

# For configurations of spatial orbitals only, the orbitals of the
# reference are valid orbitals to excite to as well.
substitution_orbitals(ref_set::Configuration, orbitals::O...) where {O<:AbstractOrbital} =
    unique(vcat(peel(ref_set).orbitals, orbitals...))

# For spin-orbitals, we do not want to excite to the orbitals of the
# reference set.
substitution_orbitals(ref_set::SpinConfiguration, orbitals::O...) where {O<:SpinOrbital} =
    filter(o -> o ∉ ref_set, O[orbitals...])

function substitutions!(fun, excitations, ref_set, orbitals,
                        min_excitations, max_excitations,
                        min_occupancy, max_occupancy)
    excite_from = 1
    retain = 1
    for i in 1:max_excitations
        i ≤ min_excitations && (retain = length(excitations)+1)
        single_excitations!(fun, excitations, ref_set, orbitals,
                            min_occupancy, max_occupancy, excite_from)
        excite_from = length(excitations)-excite_from
    end
    deleteat!(excitations, 1:retain-1)
end

function substitutions!(fun, excitations::Vector{<:SpinConfiguration},
                        ref_set::SpinConfiguration, orbitals::Vector{<:SpinOrbital},
                        min_excitations, max_excitations,
                        min_occupancy, max_occupancy)
    subs_orbs = Vector{Vector{eltype(orbitals)}}()
    for orb ∈ ref_set.orbitals
        push!(subs_orbs, filter(!isnothing, map(subs_orb -> fun(subs_orb, orb), orbitals)))
    end

    # Either all substitution orbitals are the same for all source
    # orbitals, or all source orbitals have unique substitution
    # orbitals (generated using `fun`); for simplicity, we do not
    # support a mixture.
    all_same = true
    for i ∈ 2:length(subs_orbs)
        d = [isempty(setdiff(subs_orbs[i], subs_orbs[j]))
             for j in 1:i-1]
        i == 2 && !any(d) && (all_same = false)
        all(d) && all_same || !any(d) && !all_same ||
            throw(ArgumentError("All substitution orbitals have to equal or non-overlapping"))
    end

    source_orbitals = findall(iszero, min_occupancy)

    for i ∈ min_excitations:max_excitations
        # Loop over all slots, e.g. all spin-orbitals we want to
        # excite from.
        for srci ∈ combinations(source_orbitals,i)
            all_substitutions = if all_same
                # In this case, in the first slot, any of the i..n
                # substitution orbitals goes, in the next i+1..n, in
                # the second i+2..n, &c.
                combinations(subs_orbs[1], i)
            else
                # In this case, we simply form the direct product of
                # the source-orbital specific substitution orbitals.
                Iterators.product([subs_orbs[src] for src in srci]...)
            end
            for substitutions in all_substitutions
                cfg = ref_set
                for (j,subs_orb) in enumerate(substitutions)
                    cfg = replace(cfg, ref_set.orbitals[srci[j]] => subs_orb)
                end
                push!(excitations, cfg)
            end
        end
    end
    deleteat!(excitations, 1)
end

"""
    excited_configurations([fun::Function, ] cfg::Configuration,
                           orbitals::AbstractOrbital...
                           [; min_excitations=0, max_excitations=:doubles,
                            min_occupancy=[0, 0, ...], max_occupancy=[..., g_i, ...],
                            keep_parity=true])

Generate all excitations from the reference set `cfg` by substituting
at least `min_excitations` and at most `max_excitations` of the
substitution `orbitals`. `min_occupancy` specifies the minimum
occupation number for each of the source orbitals (default `0`) and
equivalently `max_occupancy` specifies the maximum occupation number
(default is the degeneracy for each orbital). `keep_parity` controls
whether the excited configuration has to have the same parity as
`cfg`. Finally, `fun` allows modification of the substitution orbitals
depending on the source orbitals, which is useful for generating
ionized configurations. If `fun` returns `nothing`, that particular
substitution will be rejected.

# Examples

```jldoctest
julia> excited_configurations(c"1s2", o"2s", o"2p")
4-element Array{Configuration{Orbital{Int64}},1}:
 1s²
 1s 2s
 2s²
 2p²

julia> excited_configurations(c"1s2 2p", o"2p")
2-element Array{Configuration{Orbital{Int64}},1}:
 1s² 2p
 2p³

julia> excited_configurations(c"1s2 2p", o"2p", max_occupancy=[2,2])
1-element Array{Configuration{Orbital{Int64}},1}:
 1s² 2p

julia> excited_configurations(first(scs"1s2"), sos"k[s]"...) do dst,src
           if isbound(src)
               # Generate label that indicates src orbital,
               # i.e. the resultant hole
               SpinOrbital(Orbital(Symbol("[\$(src)]"), dst.orb.ℓ), dst.m)
           else
               dst
           end
       end
9-element Array{Configuration{SpinOrbital{#s16,Tuple{Int64,HalfIntegers.Half{Int64}}} where #s16<:Orbital},1}:
 1s₀α 1s₀β
 [1s₀α]s₀α 1s₀β
 [1s₀α]s₀β 1s₀β
 1s₀α [1s₀β]s₀α
 1s₀α [1s₀β]s₀β
 [1s₀α]s₀α [1s₀β]s₀α
 [1s₀α]s₀β [1s₀β]s₀α
 [1s₀α]s₀α [1s₀β]s₀β
 [1s₀α]s₀β [1s₀β]s₀β

julia> excited_configurations((a,b) -> a.m == b.m ? a : nothing,
                              spin_configurations(c"1s"), sos"k[s-d]"..., keep_parity=false)
8-element Array{Configuration{SpinOrbital{#s16,Tuple{Int64,HalfIntegers.Half{Int64}}} where #s16<:Orbital},1}:
 1s₀α
 ks₀α
 kp₀α
 kd₀α
 1s₀β
 ks₀β
 kp₀β
 kd₀β
```
"""
function excited_configurations(fun::Function,
                                ref_set::Configuration{O},
                                orbitals...;
                                min_excitations::Int=zero(Int),
                                max_excitations::Union{Int,Symbol}=:doubles,
                                min_occupancy::Vector{Int}=zeros(Int, length(peel(ref_set))),
                                max_occupancy::Vector{Int}=[degeneracy(first(o)) for o in peel(ref_set)],
                                keep_parity::Bool=true) where {O<:AbstractOrbital}
    if max_excitations isa Symbol
        max_excitations = if max_excitations == :singles
            1
        elseif max_excitations == :doubles
            2
        else
            throw(ArgumentError("Invalid maximum excitations specification $(max_excitations)"))
        end
    elseif max_excitations < 0
        throw(ArgumentError("Invalid maximum excitations specification $(max_excitations)"))
    end

    min_excitations ≥ 0 && min_excitations ≤ max_excitations ||
        throw(ArgumentError("Invalid minimum excitations specification $(min_excitations)"))

    lp = length(peel(ref_set))
    length(min_occupancy) == lp ||
        throw(ArgumentError("Need to specify $(lp) minimum occupancies for active subshells: $(peel(ref_set))"))
    length(max_occupancy) == lp ||
        throw(ArgumentError("Need to specify $(lp) maximum occupancies for active subshells: $(peel(ref_set))"))

    all(min_occupancy .>= 0) || throw(ArgumentError("Invalid minimum occupancy requested"))
    all(min_occupancy .<= max_occupancy .<= [degeneracy(first(o)) for o in peel(ref_set)]) ||
        throw(ArgumentError("Invalid maximum occupancy requested"))

    ref_set_core = core(ref_set)
    ref_set_peel = peel(ref_set)

    orbitals = substitution_orbitals(ref_set, orbitals...)

    Cfg = Configuration{promote_type(O,eltype(orbitals))}
    excitations = Cfg[ref_set_peel]
    substitutions!(fun, excitations, ref_set_peel, orbitals,
                   min_excitations, max_excitations,
                   min_occupancy, max_occupancy)
    keep_parity && filter!(e -> parity(e) == parity(ref_set), excitations)

    Cfg[ref_set_core + e for e in excitations]
end

excited_configurations(fun::Function, ref_set::Vector{<:Configuration}, args...; kwargs...) =
    unique(reduce(vcat, map(r -> excited_configurations(fun, r, args...; kwargs...), ref_set)))

default_substitution(subs_orb,orb) = subs_orb

excited_configurations(ref_set::Configuration, args...; kwargs...) =
    excited_configurations(default_substitution, ref_set, args...; kwargs...)

excited_configurations(ref_set::Vector{<:Configuration}, args...; kwargs...) =
    excited_configurations(default_substitution, ref_set, args...; kwargs...)

ion_continuum(neutral::Configuration{<:Orbital{<:Integer}},
              continuum_orbitals::Vector{Orbital{Symbol}},
              max_excitations=:singles) =
                  excited_configurations(neutral, continuum_orbitals...;
                                         max_excitations=max_excitations, keep_parity=false)

export excited_configurations, ion_continuum


text:
const UNITSQUARE = [0.0, 1.0, 0.0, 1.0]
const NULLPAIR = (0.0, 0.0)

_min(a) = minimum(filter(!isnan, a))
_max(a) = maximum(filter(!isnan, a))

function hasnan(a)
    for el ∈ a
        (el === NaN || el === missing) && return true
    end
    return false
end

lookup(item::Real, dict) = Int(item)
lookup(item::T, dict::Dict{T,V}) where {T} where V = dict[item]

# Fetch example from filename and return it as String
_example(name) = read(joinpath(dirname(@__FILE__), "../examples/docstrings", "$name.jl"), String)


text:
module LifeValData

using XLSX, DataFrames, JuliaDB, DelimitedFiles

using LifeValStructures

export Initialise, Read_Data

# set some global constants that need to be worked into functions at some point
function Initialise(filePath, fileName)
    main = DataFrame(XLSX.readtable(fileName, "MAIN")...)
    runsettings = DataFrame(XLSX.readtable(fileName, main.RS[1])...)
    spcode = DataFrame(XLSX.readtable(fileName, main.SP[1])...)

        #writes these as lastrun SETTINGS
    XLSX.writetable(filePath * "lastrun.xlsx", MAIN =( DataFrames._columns(main), DataFrames._names(main) ),
        TEMP1 =( DataFrames._columns(runsettings), DataFrames._names(runsettings) ),
        TEMP2 =( DataFrames._columns(spcode), DataFrames._names(spcode) ), overwrite=true)

    XLSX.openxlsx(filePath * "lastrun.xlsx", mode="rw") do xf
        sheet = xf[2]
        XLSX.rename!(sheet, main.RS[1])
        sheet = xf[3]
        XLSX.rename!(sheet, main.SP[1])
    end

    #read bases - this should happen before reading policies for e.g. impact of unit pricing on policies
    qxname = Vector{String}()
    wxname = Vector{String}()
    vtname = Vector{String}()
    expname = Vector{String}()
    chname = Vector{String}()
    fname = Vector{String}()
    fbname = Vector{String}()

    #qxLookup = unique(spcode, :qxbasis)
    for runtype in eachrow(runsettings)
        qxn = filePath * runtype.qxname
        wxn = filePath * runtype.wxname
        vtn = filePath * runtype.vtname
        expn = filePath * runtype.expname
        chn = filePath * runtype.chname
        fn = filePath * runtype.fname
        fbn = runtype.fbasis

        push!(qxname, qxn)
        push!(wxname, wxn)
        push!(vtname, vtn)
        push!(expname, expn)
        push!(chname, chn)
        push!(fname, fn)
        push!(fbname, fbn)
    end

    names = DataFrame([qxname, wxname, vtname, expname, chname, fname, fbname])
    rename!(names, [:qxname,:wxname,:vtname,:expname,:chname,:fname,:fbname])

    outpn = main.outv[1]

    return [runsettings, spcode, names, outpn]
end

function Read_Data(runsettings, spcode, names, filePath)
    #read policyData - []
    policyData = Vector{IndexedTable}(undef, nrow(runsettings))
    fB = DataFrame(XLSX.readtable(names.fname[1], names.fbname[1])...) #temporary, assume all f has same #s
    fData = Array{LifeValStructures.fTables}(undef, nrow(runsettings), nrow(fB))
    for runtype in eachrow(runsettings)
        policyData[runtype.runno] = load(filePath * runtype.policyfile)

        #fund data
        if runtype.fbasis != ""
            fB = DataFrame(XLSX.readtable(names.fname[runtype.runno], runtype.fbasis)...)
            for fbItem in eachrow(fB)
                fData[runtype.runno, fbItem.ID]=LifeValStructures.fTables(runtype.runno, fbItem.ID, DataFrame(fbItem))
            end
        end
    end

    #find all (unique) spcodes
    spcodes = []
    for runno in eachrow(spcode)
        spcod = runno.spc
        #if spcod !in spcodes
            push!(spcodes, spcod)
        #end
    end
    spcodes = unique(spcodes)

    qxData = Array{LifeValStructures.qxTables}(undef, nrow(runsettings), 100)
    qxpData = Array{LifeValStructures.qxTables}(undef, nrow(runsettings), 100)
    wxData = Array{LifeValStructures.wxTables}(undef, nrow(runsettings), 100)
    wxpData = Array{LifeValStructures.wxTables}(undef, nrow(runsettings), 100)
    vtData = Array{LifeValStructures.vtTables}(undef, nrow(runsettings), 100)
    vtpData = Array{LifeValStructures.vtTables}(undef, nrow(runsettings), 100)
    expData = Array{LifeValStructures.expTables}(undef, nrow(runsettings), 100)
    exppData = Array{LifeValStructures.expTables}(undef, nrow(runsettings), 100)
    chData = Array{LifeValStructures.chTables}(undef, nrow(runsettings), 100)

    for runno in eachrow(spcode)
        runtemp = runno.runno
        spcod = runno.spc

        qxB = DataFrame(XLSX.readtable(names.qxname[runtemp], runno.qxbasis)...)
        qxData[runtemp,spcod]=LifeValStructures.qxTables(runtemp,spcod,qxB)

        if runno.qxbasisp != ""
            qxpB = DataFrame(XLSX.readtable(names.qxname[runtemp], runno.qxbasisp)...)
            qxpData[runtemp,spcod]=LifeValStructures.qxTables(runtemp,spcod,qxpB)
        end

        wxB = DataFrame(XLSX.readtable(names.wxname[runtemp], runno.wxbasis)...)
        wxData[runtemp,spcod]=LifeValStructures.wxTables(runtemp,spcod,wxB)

        if runno.wxbasisp != ""
            wxpB = DataFrame(XLSX.readtable(names.wxname[runtemp], runno.wxbasisp)...)
            wxpData[runtemp,spcod]=LifeValStructures.wxTables(runtemp,spcod,wxpB)
        end

        vtB = DataFrame(XLSX.readtable(names.vtname[runtemp], runno.vtbasis)...)
        vtData[runtemp,spcod]=LifeValStructures.vtTables(runtemp,spcod,vtB)

        if runno.vtbasisp != ""
            vtpB = DataFrame(XLSX.readtable(names.vtname[runtemp], runno.vtbasisp)...)
            vtpData[runtemp,spcod]=LifeValStructures.vtTables(runtemp,spcod,vtpB)
        end

        expB = DataFrame(XLSX.readtable(names.expname[runtemp], runno.ebasis)...)
        expItem = filter(row -> row.SPCODE == runno.spc, expB) #don't need every spcode on spcode level
        expData[runtemp,spcod]=LifeValStructures.expTables(runtemp,spcod,expItem)

        if runno.ebasisp != ""
            exppB = DataFrame(XLSX.readtable(names.expname[runtemp], runno.ebasisp)...)
            exppItem = filter(row -> row.SPCODE == runno.spc, exppB) #don't need every spcode on spcode level
            exppData[runtemp,spcod]=LifeValStructures.expTables(runtemp,spcod,exppItem)
        end

        if runno.chbasis != ""
            chB = DataFrame(XLSX.readtable(names.chname[runtemp], runno.chbasis)...)
            chData[runtemp,spcod]=LifeValStructures.chTables(runtemp,spcod,chB)
        end
    end

    basisData = [qxData, qxpData, wxData, wxpData, vtData, vtpData, expData, exppData, chData]
    #names!(basisData, [:qxData,:qxpData,:wxData,:wxpData,:vtData,:vtpData,:expData,:exppData,:chData,:fData])

    return [policyData, basisData, fData]
end

end