text:
using Plots, Statistics, FFTW, SparseArrays, DelimitedFiles, ArgParse, LinearAlgebra#, JLD

#=
#Loading Functions
=#

# function arrange_arrays(data, N)
#     #=
#     This arranges the data as vectors of x,y,z,vx,vy and vz
#     =#
#     x_indices = [2 + 6*i for i in 0:N-1]
#     y_indices = [3 + 6*i for i in 0:N-1]
#     z_indices = [4 + 6*i for i in 0:N-1]
#     vx_indices = [5 + 6*i for i in 0:N-1]
#     vy_indices = [6 + 6*i for i in 0:N-1]
#     vz_indices = [7 + 6*i for i in 0:N-1]
#
#     x, y, z = data[:,x_indices], data[:,y_indices], data[:,z_indices]
#     vx, vy, vz = data[:,vx_indices], data[:,vy_indices], data[:,vz_indices]
#     x, y, z, vx, vy, vz
# end

# function arrange_arrays(idx::Int, N)
#     pos, vel = load_data(i=idx, N=N)
#
#     x = hcat([pos[:,1+3*i] for i in 0:N-1]...)
#     y = hcat([pos[:,2+3*i] for i in 0:N-1]...)
#     z = hcat([pos[:,3+3*i] for i in 0:N-1]...)
#     vx = hcat([vel[:,1+3*i] for i in 0:N-1]...)
#     vy = hcat([vel[:,2+3*i] for i in 0:N-1]...)
#     vz = hcat([vel[:,3+3*i] for i in 0:N-1]...)
#     x, y, z, vx, vy, vz
# end

# function load_data(; i=1, path=PATH, N=num_part)
#     cd(path)
#     files = readdir(path * "/$i")
#     positions = readdlm(path * "/$i/positions.dat")
#     velocities = readdlm(path * "/$i/velocity.dat")
#     Array(reshape(positions, 3N, :)'), Array(reshape(velocities, 3N, :)')
# end

function parseCommandLine(; dir=0)

        # initialize the settings (the description is for the help screen)
        # s = ArgParseSettings(description = "Example 1 for argparse.jl: minimal usage.")
		s = ArgParseSettings(description = "Dynamical Structure Factor code.")

        @add_arg_table! s begin
           "--Filename", "-f"         # another option, with short form
                help = "File name"
                arg_type = String
                default = "Ge00Sb00Te100T823K.cpmd"
			"--pathIn", "-i"         # another option, with short form
                 help = "Directory where input file is stored."
                 arg_type = String
                 default = pwd()
			"--pathOut", "-o"         # another option, with short form
                 help = "Directory where everything will be saved. This directory can be created before or during execution."
                 arg_type = String
                 default = "temp"
           "--n", "-N"
                    help = "number of particles"
					arg_type = Int
                    default = 300
           "--Lchar", "-L"
                help = "box length"
                arg_type = Float32
                default = 42.2976f0
           "--t_max", "-t"
                help = "t maximum"
                arg_type = Int64
                default = 15000
			"--dt", "-s"
                 help = "time step between configurations"
                 arg_type = Float32
                 default = 10.0f0
			 "--dim", "-d"
	 			 help = "dimensions of time step between configurations"
	 			 arg_type = String
	 			 default = "ps"
			"--sparse", "-b"
				help = "Use sparse matrices? Should be faster..."
				arg_type = Int
				default = 0
			"--time_window", "-w"
				help = "List of time window of samples. Should be a  space separated"
				nargs = '*'
				arg_type = Int64
				default = [500, 1000]
	   end

        return parse_args(s) # the result is a Dict{String,Any}
end

function init_dict(path_to_file, pathout, t_MAX, L_char, num_part, title, Δt, dim, t_samp_list)
    # dict_path = PATH_out * "Dicts"
    # isdir(dict_path) || mkdir(dict_path)
    # dict_path = PATH_out * "Dicts/$(parsed_args["n"])"
    # isdir(dict_path) || mkdir(dict_path)
    dict = Dict()
    dict[:PathIn] = path_to_file
    dict[:PathOut] = pathout
	dict[:t_max] = t_MAX
    dict[:L_char] = L_char
	dict[:num_part] = num_part
    dict[:title] = title
	dict[:dt] = Δt
	dict[:dim] = dim
	dict[:t_samp_list] = t_samp_list

    # @save dict[:Path] dict
    return dict
end

function load_data(N, t_max; path = PATH_TO_FILE)
    data1 = readdlm(PATH_TO_FILE)
    x = reshape(data1[:,2],N,:)'
    y = reshape(data1[:,3],N,:)'
    z = reshape(data1[:,4],N,:)'
    vx = reshape(data1[:,5],N,:)'
    vy = reshape(data1[:,6],N,:)'
    vz = reshape(data1[:,7],N,:)'
    # samp = Int(floor(t_max/t_sample))
    # d = [(x[r,:],y[r,:],z[r,:],vx[r,:],vy[r,:],vz[r,:]) for r in Iterators.partition(1:t_max,t_sample)][1:samp]
    @info "Data loaded"
    # t_max = t_max == 0 ? end : t_max
    d = t_max != 0 ? (x[1:t_max,:],y[1:t_max,:],z[1:t_max,:],vx[1:t_max,:],vy[1:t_max,:],vz[1:t_max,:]) : (x,y,z,vx,vy,vz)
end

function partition_data(d, t_sample)
    t_max = size(d[1],1)
    samp = Int(floor(t_max/t_sample))
    data = [(d[1][r,:],d[2][r,:],d[3][r,:],d[4][r,:],d[5][r,:],d[6][r,:]) for r in Iterators.partition(1:t_max,t_sample)][1:samp]
    @info "Data partitioned"
    data, samp * t_sample
end
#=
Tests
=#
function test_data(;Δt = 0.01, steps = 5000, thermal_steps=1000,
			N=300, L = 1, dims=3, arrange = "Crystal", typeHessian="random")

	r, v, M = generateData(N; steps, dims, Δt, arrange, typeHessian)
	((r[1,thermal_steps+1:end,:]/L, r[2,thermal_steps+1:end,:]/L, r[3,thermal_steps+1:end,:]/L,
	v[1,thermal_steps+1:end,:], v[2,thermal_steps+1:end,:], v[3,thermal_steps+1:end,:]), M)
end
#=
#Correlation Functions
=#
function gen_Jays(data, N; kmin = 1, kmax = 10, L = 1, m=1)
    #=
    Generates Jₓ(k,t) = m/N * Σ vₓⁱ(t) exp(-𝚤 k ⋅ rⁱ(t))
    As a matter of sanity check, if data = rand(tmax, 6N+1) where
    mean(rand()) = 1/2, then Jₓ(0,t) = m/2
    =#
    x, y, z, vx, vy, vz = data
    #Transversal to z
    Jxkz = hcat([m/√N*sum(vx .* cos.(2π * k .* z / L), dims=2) for k in kmin:kmax]...)
    JxkzIm = hcat([m/√N*sum(vx .* sin.(2π * k .* z / L), dims=2) for k in kmin:kmax]...)

    Jykz = hcat([m/√N*sum(vy .* cos.(2π * k .* z / L), dims=2) for k in kmin:kmax]...)
    JykzIm = hcat([m/√N*sum(vy .* sin.(2π * k .* z / L), dims=2) for k in kmin:kmax]...)
    #Transversal to x
    Jzkx = hcat([m/√N*sum(vz .* cos.(2π * k .* x / L), dims=2) for k in kmin:kmax]...)
    JzkxIm = hcat([m/√N*sum(vz .* sin.(2π * k .* x / L), dims=2) for k in kmin:kmax]...)

    Jykx = hcat([m/√N*sum(vy .* cos.(2π * k .* x / L), dims=2) for k in kmin:kmax]...)
    JykxIm = hcat([m/√N*sum(vy .* sin.(2π * k .* x / L), dims=2) for k in kmin:kmax]...)
    #Transversal to y
    Jxky = hcat([m/√N*sum(vx .* cos.(2π * k .* y / L), dims=2) for k in kmin:kmax]...)
    JxkyIm = hcat([m/√N*sum(vx .* sin.(2π * k .* y / L), dims=2) for k in kmin:kmax]...)

    Jzky = hcat([m/√N*sum(vz .* cos.(2π * k .* y / L), dims=2) for k in kmin:kmax]...)
    JzkyIm = hcat([m/√N*sum(vz .* sin.(2π * k .* y / L), dims=2) for k in kmin:kmax]...)

    Jxkz, Jykz, Jxky, Jzky, Jykx, Jzkx, JxkzIm, JykzIm, JxkyIm, JzkyIm, JykxIm, JzkxIm
end

function gen_JiJi(data, N, tmax; kmin = 1, kmax = 10, L = 1, m=1)
    #=
    Generates ⟨ JxJx ⟩(k,τ) =  Σₜ Jₓ(k,t+τ)Jₓ(k,t)/Z where Z = tmax-τ
    As a matter of sanity check, if data = rand(tmax, 6N+1) where
    mean(rand()) = 1/2, then ⟨ JxJx ⟩(k,τ) = m²/4
    =#
    # Jxkz, Jykz, Jxky, Jzky, Jykx, Jzkx = gen_Jays(data, N; kmin = kmin, kmax = kmax, L = L, m=m)
    Jxkz, Jykz, Jxky, Jzky, Jykx, Jzkx, JxkzIm, JykzIm, JxkyIm, JzkyIm,
        JykxIm, JzkxIm = gen_Jays(data, N; kmin = kmin, kmax = kmax, L = L, m=m)

    JxJxkz = hcat([sum([Jxkz[τ+t,:] .* Jxkz[t,:] .+ JxkzIm[τ+t,:] .* JxkzIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    JxJxky = hcat([sum([Jxky[τ+t,:] .* Jxky[t,:] .+ JxkyIm[τ+t,:] .* JxkyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JxJxkzIm = hcat([sum([Jxkz[τ+t,:] .* JxkzIm[t,:] .+ JxkzIm[τ+t,:] .* Jxkz[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    JxJxkyIm = hcat([sum([Jxky[τ+t,:] .* JxkyIm[t,:] .+ JxkyIm[τ+t,:] .* Jxky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JyJykz = hcat([sum([Jykz[τ+t,:] .* Jykz[t,:] .+ JykzIm[τ+t,:] .* JykzIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    JyJykx = hcat([sum([Jykx[τ+t,:] .* Jykx[t,:] .+ JykxIm[τ+t,:] .* JykxIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JyJykzIm = hcat([sum([Jykz[τ+t,:] .* JykzIm[t,:] .+ JykzIm[τ+t,:] .* Jykz[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    JyJykxIm = hcat([sum([Jykx[τ+t,:] .* JykxIm[t,:] .+ JykxIm[τ+t,:] .* Jykx[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JzJzkx = hcat([sum([Jzkx[τ+t,:] .* Jzkx[t,:] .+ JzkxIm[τ+t,:] .* JzkxIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    JzJzky = hcat([sum([Jzky[τ+t,:] .* Jzky[t,:] .+ JzkyIm[τ+t,:] .* JzkyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JzJzkxIm = hcat([sum([Jzkx[τ+t,:] .* JzkxIm[t,:] .+ JzkxIm[τ+t,:] .* Jzkx[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
    JzJzkyIm = hcat([sum([Jzky[τ+t,:] .* JzkyIm[t,:] .+ JzkyIm[τ+t,:] .* Jzky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JxJxkz, JxJxky, JyJykz, JyJykx, JzJzkx, JzJzky, JxJxkzIm, JxJxkyIm, JyJykzIm, JyJykxIm, JzJzkxIm, JzJzkyIm
end

function gen_JiJi(data, N, tmax, β::Int; kmin = 1, kmax = 10, L = 1, m=1)
    #=
    Generates ⟨ JxJx ⟩(k,τ) =  Σₜ Jₓ(k,t+τ)Jₓ(k,t)/Z where Z = tmax-τ
    As a matter of sanity check, if data = rand(tmax, 6N+1) where
    mean(rand()) = 1/2, then ⟨ JxJx ⟩(k,τ) = m²/4
    =#
    # Jxkz, Jykz, Jxky, Jzky, Jykx, Jzkx = gen_Jays(data, N; kmin = kmin, kmax = kmax, L = L, m=m)
    Jxkz, Jykz, Jxky, Jzky, Jykx, Jzkx, JxkzIm, JykzIm, JxkyIm, JzkyIm,
        JykxIm, JzkxIm = gen_Jays(data, N; kmin = kmin, kmax = kmax, L = L, m=m)

	JxJxkz = hcat([diag((Jxkz' * sparse_with_diag(i,tmax) * Jxkz .+ JxkzIm' * sparse_with_diag(i,tmax) * JxkzIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JxJxkz = hcat([sum([Jxkz[τ+t,:] .* Jxkz[t,:] .+ JxkzIm[τ+t,:] .* JxkzIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	JxJxky = hcat([diag((Jxky' * sparse_with_diag(i,tmax) * Jxky .+ JxkyIm' * sparse_with_diag(i,tmax) * JxkyIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JxJxky = hcat([sum([Jxky[τ+t,:] .* Jxky[t,:] .+ JxkyIm[τ+t,:] .* JxkyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JxJxkzIm = hcat([diag((Jxkz' * sparse_with_diag(i,tmax) * JxkzIm .+ JxkzIm' * sparse_with_diag(i,tmax) * Jxkz)/(tmax-i)) for i in 0:tmax-1]...)
    # JxJxkzIm = hcat([sum([Jxkz[τ+t,:] .* JxkzIm[t,:] .+ JxkzIm[τ+t,:] .* Jxkz[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	JxJxkyIm = hcat([diag((Jxky' * sparse_with_diag(i,tmax) * JxkyIm .+ JxkyIm' * sparse_with_diag(i,tmax) * Jxky)/(tmax-i)) for i in 0:tmax-1]...)
    # JxJxkyIm = hcat([sum([Jxky[τ+t,:] .* JxkyIm[t,:] .+ JxkyIm[τ+t,:] .* Jxky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JyJykz = hcat([diag((Jykz' * sparse_with_diag(i,tmax) * Jykz .+ JykzIm' * sparse_with_diag(i,tmax) * JykzIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JyJykz = hcat([sum([Jykz[τ+t,:] .* Jykz[t,:] .+ JykzIm[τ+t,:] .* JykzIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	JyJykx = hcat([diag((Jykx' * sparse_with_diag(i,tmax) * Jykx .+ JykxIm' * sparse_with_diag(i,tmax) * JykxIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JyJykx = hcat([sum([Jykx[τ+t,:] .* Jykx[t,:] .+ JykxIm[τ+t,:] .* JykxIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JyJykzIm = hcat([diag((Jykz' * sparse_with_diag(i,tmax) * JykzIm .+ JykzIm' * sparse_with_diag(i,tmax) * Jykz)/(tmax-i)) for i in 0:tmax-1]...)
    # JyJykzIm = hcat([sum([Jykz[τ+t,:] .* JykzIm[t,:] .+ JykzIm[τ+t,:] .* Jykz[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	JyJykxIm = hcat([diag((Jykx' * sparse_with_diag(i,tmax) * JykxIm .+ JykxIm' * sparse_with_diag(i,tmax) * Jykx)/(tmax-i)) for i in 0:tmax-1]...)
    # JyJykxIm = hcat([sum([Jykx[τ+t,:] .* JykxIm[t,:] .+ JykxIm[τ+t,:] .* Jykx[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JzJzkx = hcat([diag((Jzkx' * sparse_with_diag(i,tmax) * Jzkx .+ JzkxIm' * sparse_with_diag(i,tmax) * JzkxIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JzJzkx = hcat([sum([Jzkx[τ+t,:] .* Jzkx[t,:] .+ JzkxIm[τ+t,:] .* JzkxIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	JzJzky = hcat([diag((Jzky' * sparse_with_diag(i,tmax) * Jzky .+ JzkyIm' * sparse_with_diag(i,tmax) * JzkyIm)/(tmax-i)) for i in 0:tmax-1]...)
    # JzJzky = hcat([sum([Jzky[τ+t,:] .* Jzky[t,:] .+ JzkyIm[τ+t,:] .* JzkyIm[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

	JzJzkxIm = hcat([diag((Jzkx' * sparse_with_diag(i,tmax) * JzkxIm .+ JzkxIm' * sparse_with_diag(i,tmax) * Jzkx)/(tmax-i)) for i in 0:tmax-1]...)
    # JzJzkxIm = hcat([sum([Jzkx[τ+t,:] .* JzkxIm[t,:] .+ JzkxIm[τ+t,:] .* Jzkx[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)
	JzJzkyIm = hcat([diag((Jzky' * sparse_with_diag(i,tmax) * JzkyIm .+ JzkyIm' * sparse_with_diag(i,tmax) * Jzky)/(tmax-i)) for i in 0:tmax-1]...)
    # JzJzkyIm = hcat([sum([Jzky[τ+t,:] .* JzkyIm[t,:] .+ JzkyIm[τ+t,:] .* Jzky[t,:] for t in 1:tmax-τ])/(tmax-τ) for τ in 0:tmax-1]...)

    JxJxkz, JxJxky, JyJykz, JyJykx, JzJzkx, JzJzky, JxJxkzIm, JxJxkyIm, JyJykzIm, JyJykxIm, JzJzkxIm, JzJzkyIm
end

function sparse_with_diag(idx, tmax)
    A = blockdiag(sparse([i for i in 1+idx:tmax],[i for i in 1:tmax-idx],[1 for i in 1:tmax-idx]), sparse(0I, 0, idx))
    A
end

function get_Corr(data, N, tmax; kmin = 1, kmax = 10, L = 1, m=1)
	if use_sparse == 0
    	CT = gen_JiJi(data, N, tmax; kmin = kmin, kmax = kmax, L = L, m=m)
		return mean(CT[1:6]), mean(CT[7:12])
	elseif use_sparse == 1
		CT = gen_JiJi(data, N, tmax,1; kmin = kmin, kmax = kmax, L = L, m=m)
		return  mean(CT[1:6]), mean(CT[7:12])
	end
end

function get_Corr_over_all_samples(t_samp, N, t_max, d; kmin = 1, kmax = 10,
			L = 1, m=1, in_Fourier=false)
    tmax = t_samp
    samp = Int(floor(t_max/t_samp))
    # d = load_data(N,t_samp,t_max)
    data = d[1]
    CT, CTIm = get_Corr(data, N, tmax; kmin = kmin, kmax = kmax, L = L, m=m)
    CT_av, CTIm_av = reshape(CT,1,kmax-kmin+1,tmax), reshape(CTIm,1,kmax-kmin+1,tmax)
    for i in 2:samp
        @info "Sample $i"
        data = d[i]
        CT, CTIm = get_Corr(data, N, tmax; kmin = kmin, kmax = kmax, L = L, m=m)

        CT_av = vcat(CT_av, reshape(CT,1,kmax-kmin+1,tmax))
        CTIm_av = vcat(CTIm_av, reshape(CTIm,1,kmax-kmin+1,tmax))
    end
	# return CT_av, CTIm_av
	do_the_fou(CT_av, CTIm_av, tmax; in_Fourier, L)
end

function do_the_fou(CT_av, CTIm_av, tmax; kmin = 1, kmax = 10,
			L = 1, m=1, in_Fourier=true)
	if in_Fourier
        CT_av, CTIm_av = reshape(mean(CT_av,dims=1),kmax-kmin+1,tmax), reshape(mean(CTIm_av,dims=1),kmax-kmin+1,tmax)
        CT_FT2 = abs.(hcat(fft.([CT_av[i,:] for i in 1:kmax-kmin+1])...))
        CTIm_FT = abs.(hcat(fft.([CTIm_av[i,:] for i in 1:kmax-kmin+1])...))
        return CT_FT2, CTIm_FT
    else
        return reshape(mean(CT_av,dims=1),kmax-kmin+1,tmax), reshape(mean(CTIm_av,dims=1),kmax-kmin+1,tmax)
    end
end

# function get_Corr_over_all_samples_in_fourier(t_samp, N, t_max, d; kmin = 1, kmax = 10, L = 1, m=1)
#     tmax = t_samp
#     samp = Int(floor(t_max/t_samp))
#     # d = load_data(N,t_samp,t_max)
#     data = d[1]
#     CT, CTIm = get_Corr(data, N, tmax; kmin = kmin, kmax = kmax, L = L, m=m)
#     CT_av, CTIm_av = reshape(CT,1,kmax-kmin+1,tmax), reshape(CTIm,1,kmax-kmin+1,tmax)
#     for i in 2:samp
#         @info "Sample $i"
#         data = d[i]
#         CT, CTIm = get_Corr(data, N, tmax; kmin = kmin, kmax = kmax, L = L, m=m)
#
#         CT_av = vcat(CT_av, reshape(CT,1,kmax-kmin+1,tmax))
#         CTIm_av = vcat(CTIm_av, reshape(CTIm,1,kmax-kmin+1,tmax))
#     end
# 	apply_fourier(CT_av, samp, kmax-kmin+1), apply_fourier(CTIm_av, samp, kmax-kmin+1)
# end

# function apply_fourier(corr, samp, knum)
#     ff = fft.([corr[i,j,:] for i in 1:samp, j in 1:knum])
#     ff_ar = cat([hcat(ff[i,:]...) for i in 1:samp]..., dims=3)
#     ff_ar_m = reshape(mean(ff_ar, dims=3),:,knum)
#     fp = abs.(ff_ar_m)
#
#     fp
# end

# function get_ω_sample(N, t_max, data, t_samp_list; kmin = 1, kmax = 10, L = 1,
#         m=1, Δs=10) #ωlim = 30 data_samp_freq=10
# 	data_samp_period=dict[:dt]
#     tsamp = Int(floor(t_max/Δs))
#     t_samp_list = t_samp_list == 0 ? [i*tsamp for i in 1:Δs] : t_samp_list
#     n_samp = Int.(floor.(t_max ./ t_samp_list))
#     ω_list = zeros(size(t_samp_list,1),kmax-kmin + 1)
# 	CT_arr = zeros(maximum(t_samp_list),kmax-kmin + 1,2,size(t_samp_list,1))
#
#     for j in 1:size(t_samp_list,1)
#         ts = t_samp_list[j]
#         d_part, t_m = partition_data(data, ts)
#         @info "$j : t_samp set to $ts"
#         if ts * n_samp[j] != t_m
#             @warn "t_sample x n_sample does NOT equal t_m"
#             break
#         end
#         ωlim = Int(floor(ts/2))
#         CT, _ = get_Corr_over_all_samples_in_fourier(ts, N, t_max, d_part; kmin = kmin, kmax = kmax, L = L, m=m)
#         ω_list[j,:] = vcat([findall(x->x == maximum(CT[1:ωlim,i]), CT[1:ωlim,i]) for i in 1:size(CT,2)]...)
#
# 		freq = hcat([[l for l in 0:size(CT,1)-1] ./ (ts * data_samp_period) for i in 1:size(CT,2)]...)
# 		CT_arr[1:size(CT,1),:,:,j] = cat(CT,freq, dims=3)
#     end
#     freq_list = hcat([(ω_list[i,:] .- 1) ./(t_samp_list[i] * data_samp_period) for i in 1:size(t_samp_list,1)]...)
#     ω_list, freq_list, arrange_CT_arr(CT_arr; len=Int(floor(minimum(t_samp_list)/2)), Δk = kmax-kmin + 1)
# end

function get_ω_sample(N, t_max, data, t_samp_list, c::Bool; kmin = 1, kmax = 10, L = 1,
        m=1, Δs=10) #ωlim = 30 data_samp_period=10
	data_samp_period=dict[:dt]
    tsamp = Int(floor(t_max/Δs))
    t_samp_list = t_samp_list == 0 ? [1000,2000,5000] : t_samp_list
    n_samp = Int.(floor.(t_max ./ t_samp_list))
    ω_list = zeros(size(t_samp_list,1),kmax-kmin + 1)
	CT_arr = zeros(maximum(t_samp_list),kmax-kmin + 1,2,size(t_samp_list,1))

    for j in 1:size(t_samp_list,1)
        ts = t_samp_list[j]
        d_part, t_m = partition_data(data, ts)
        @info "$j : t_samp set to $ts"
        if ts * n_samp[j] != t_m
            @warn "t_sample x n_sample does NOT equal t_m"
            break
        end
        ωlim = Int(floor(ts/2))
        CT, _ = get_Corr_over_all_samples(ts, N, t_max, d_part; kmin = kmin, kmax = kmax, L = L, m=m, in_Fourier = c)
        ω_list[j,:] = vcat([findall(x->x == maximum(CT[1:ωlim,i]), CT[1:ωlim,i]) for i in 1:size(CT,2)]...)

		freq = hcat([[l for l in 0:size(CT,1)-1] ./ (ts * data_samp_period) for i in 1:size(CT,2)]...)
		CT_arr[1:size(CT,1),:,:,j] = cat(CT,freq, dims=3)
    end
    freq_list = hcat([(ω_list[i,:] .- 1) ./(t_samp_list[i] * data_samp_period) for i in 1:size(t_samp_list,1)]...)
    # ω_list, freq_list, arrange_CT_arr(CT_arr; len=Int(floor(minimum(t_samp_list)/2)), Δk = kmax-kmin + 1)
	CT_arrv2 = CT_arr
	CT_arr = arrange_CT_arr(CT_arr; len=Int(floor(minimum(t_samp_list)/2)), Δk = kmax-kmin + 1)
	CT_arr = envelope.([CT_arr[i,:,:] for i in 1:kmax-kmin + 1])
	corr = permutedims(cat(CT_arr...,dims=3),(3,1,2))
	ω_list, freq_list, corr, CT_arrv2
	# ω_list, freq_list, CT_arr
end

function arrange_CT_arr(CT_arr; len=30, Δk=10)
        t_len=size(CT_arr,4)
        l = [i*len for i in 1:t_len]
        new_CT = zeros(Δk,sum(l),2)
        for k in 1:Δk
                x = vcat([CT_arr[1:i*len,k,2,i] for i in 1:t_len]...)
                y = vcat([CT_arr[1:i*len,k,1,i] for i in 1:t_len]...)
                idx = sortperm(x)
                new_CT[k,:,:] = hcat([x[idx],y[idx]]...)
        end
        new_CT
end

function envelope(ar)
    arr = ar[:,1]
    arr2 = ar[:,2]

    l = zeros(size(arr))
    ind = 1
    for i in 1:size(l,1)
        if ind <= size(arr2,1)
            ind = findall(x->x == maximum(arr2[arr[ind] .== arr]), arr2)[1]
            l[i] = ind
            ind = maximum(indexin(arr2[arr[ind] .== arr], arr2)) + 1
        end
    end
    ll = Array{Int,1}(l[l .>0])
    ar[ll,:]
end

function get_freq_max_2(CT_arr)
        freq_max = zeros(10,2)
        for k in 1:10
                idx = sortperm(CT_arr[k,:,2])[end]
                freq_max[k,:] = [k, CT_arr[k,idx,1]]
        end
        freq_max
end

function get_freqs(dict; test=false, ω = 0.0002, Δt = 10.0, δ = 0.5,
			t_samp_list=0)
    #dict for t_MAX, L_char,, num_part
    t_MAX = dict[:t_max]
    L_char = dict[:L_char]
    num_part = dict[:num_part]

    if test
        d, M = test_data(steps = t_MAX + 2000, thermal_steps=2000, L = L_char,
				Δt = dict[:dt], N=num_part, arrange = "Crystal", typeHessian="random")
    else
        d = load_data(num_part, t_MAX)
    end

    get_ω_sample(num_part, t_MAX, d, t_samp_list, true; L=L_char)
	# get_ω_sample(num_part, t_MAX, d, [500,1000], true; L=L_char)
end

function saveData(ω_list, freq_list, CT, CT2, t_samp_list)
	fig_name = split(filename,".")[1]
	writedlm(dict[:PathOut] * "/" * fig_name * "_modes.dat", ω_list')
	writedlm(dict[:PathOut] * "/" * fig_name * "_freq.dat", freq_list)
	# writedlm(dict[:PathOut] * "/" * fig_name * "_CT2.dat", CT_arrv2)
	# writedlm(dict[:PathOut] * "/" * fig_name * "_CT.dat", CT_arr)

	isdir(dict[:PathOut] * "/" * fig_name * "_CorrelationDataTimeWindowsMixed") || mkdir(dict[:PathOut] * "/" * fig_name * "_CorrelationDataTimeWindowsMixed")
	for k in 1:size(CT,1)
	    writedlm(dict[:PathOut] * "/" * fig_name * "_CorrelationDataTimeWindowsMixed/CT_k_$k.dat", CT[k,:,:])
	end

	isdir(dict[:PathOut] * "/" * fig_name * "_CorrelationData") || mkdir(dict[:PathOut] * "/" * fig_name * "_CorrelationData")
	for k in 1:size(CT2,2), s in 1:size(CT2,4)
		tw = t_samp_list[s]
	    tw2 = Int(floor(tw/2))
	    writedlm(dict[:PathOut] * "/" * fig_name * "_CorrelationData/CT_k_$(k)_TW_$tw.dat", [CT2[1:tw2,k,2,s] CT2[1:tw2,k,1,s]])
	end
end

function savePlots(ω_list, freq_list, CT, CT2, t_samp_list)
	dim = dict[:dim]
	m = reshape(mean(freq_list, dims=2),:)
	st = reshape(std(freq_list, dims=2),:)
	fig = plot(1:10, m, ribbon=st, fillalpha=0.2, label="mean", lw=4,
	    	frame=:box, xlabel="k (wave number)", ylabel="freq (1/$dim)", legend=:topleft,
	    	title=dict[:title], margin = 10Plots.mm, ms=5, markershapes = [:circle], markerstrokewidth=0)

	figs_CT = []
    for k in 1:size(CT,1)
        f_CT = plot(CT[k,1:200,1],CT[k,1:200,2], frame=:box, xlabel="freq (1/$dim)", ylabel="CT", label="k = $k", margin = 5Plots.mm, ms=4, markershapes = [:circle], markerstrokewidth=0)
        push!(figs_CT, f_CT)
    end
    figs_freq = []
    figs_ω = []
    figs_CT2 = []
    for s in 1:size(ω_list,1)
		tw = t_samp_list[s]
        f_ω = plot(ω_list[s,:], frame=:box, xlabel="k", ylabel="ω max", label="Time window = $tw", margin = 5Plots.mm, ms=4, markershapes = [:circle], markerstrokewidth=0)
        push!(figs_ω, f_ω)

        f_freq = plot(freq_list[:,s], frame=:box, xlabel="k", ylabel="freq max (1/$dim)", label="Time window = $tw", margin = 5Plots.mm, ms=4, markershapes = [:circle], markerstrokewidth=0)
        push!(figs_freq, f_freq)

        for k in 1:size(CT2,2)
			tw2 = Int(floor(tw/2))
            f_CT2 = plot(CT2[1:tw2,k,2,s], CT2[1:tw2,k,1,s], frame=:box, xlabel="freq (1/$dim)", ylabel="CT", label="k = $k, Time Window $tw", margin = 5Plots.mm, ms=4, markershapes = [:circle], markerstrokewidth=0)
            push!(figs_CT2, f_CT2)
        end
    end

	fig_name = split(filename,".")[1]
	savefig(fig, dict[:PathOut] * "/" * fig_name * "_freqMean.png")
	savefig(plot(figs_CT..., size=(1500,900)), dict[:PathOut] * "/" * fig_name * "_CT.png")
	savefig(plot(figs_freq..., size=(1500,900)), dict[:PathOut] * "/" * fig_name * "_freq.png")
	savefig(plot(figs_ω..., size=(1500,900)), dict[:PathOut] * "/" * fig_name * "_modes.png")
	savefig(plot(figs_CT2..., size=(2500,1900)), dict[:PathOut] * "/" * fig_name * "_CT2.png")
end

function main(;t_samp_list=0)
	# dim = dict[:dim]
	@info "Starting..."
	ω_list, freq_list, CT_arr, CT_arrv2 = get_freqs(dict; t_samp_list)
	@info "Computation done"
	# m = reshape(mean(freq_list, dims=2),:)
	# st = reshape(std(freq_list, dims=2),:)
	# fig = plot(1:10, m, ribbon=st, fillalpha=0.2, label="mean", lw=4,
	#     	frame=:box, xlabel="k (wave number)", ylabel="freq (1/$dim)", legend=:topleft,
	#     	title=dict[:title], margin = 10Plots.mm, ms=5, markershapes = [:circle], markerstrokewidth=0)

	# fig_name = split(filename,".")[1]
	# savefig(fig, dict[:PathOut] * "/" * fig_name * ".png")
	@info "Generating Plots..."
	savePlots(ω_list, freq_list, CT_arr, CT_arrv2, t_samp_list)
	# writedlm(dict[:PathOut] * "/" * fig_name * "_modes.dat", ω_list')
	# writedlm(dict[:PathOut] * "/" * fig_name * "_freq.dat", freq_list)
	# writedlm(dict[:PathOut] * "/" * fig_name * "_Corr.dat", CT_arrv2)
	@info "Saving Data..."
	saveData(ω_list, freq_list, CT_arr, CT_arrv2, t_samp_list)
	@info "Done!"
	# writedlm(dict[:PathOut] * "/" * fig_name * "_freq_method2.dat",
	# 		get_freq_max_2(CT_arr))
end

max_stars_repo_name: jquetzalcoatl/DynamicStructureFactor

max_issues_repo_name: jquetzalcoatl/DynamicStructureFactor

max_forks_repo_name: jquetzalcoatl/DynamicStructureFactor


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

mutable struct DeploymentOperation <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    operationId::Any # spec type: Union{ Nothing, String } # spec name: operationId
    properties::Any # spec type: Union{ Nothing, DeploymentOperationProperties } # spec name: properties

    function DeploymentOperation(;id=nothing, operationId=nothing, properties=nothing)
        o = new()
        validate_property(DeploymentOperation, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(DeploymentOperation, Symbol("operationId"), operationId)
        setfield!(o, Symbol("operationId"), operationId)
        validate_property(DeploymentOperation, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type DeploymentOperation

const _property_map_DeploymentOperation = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("operationId")=>Symbol("operationId"), Symbol("properties")=>Symbol("properties"))
const _property_types_DeploymentOperation = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("operationId")=>"String", Symbol("properties")=>"DeploymentOperationProperties")
Base.propertynames(::Type{ DeploymentOperation }) = collect(keys(_property_map_DeploymentOperation))
Swagger.property_type(::Type{ DeploymentOperation }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_DeploymentOperation[name]))}
Swagger.field_name(::Type{ DeploymentOperation }, property_name::Symbol) =  _property_map_DeploymentOperation[property_name]

function check_required(o::DeploymentOperation)
    true
end

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

ext: jl


text:
## 08.jl
# Typoglycemia

function typoglycemia(words)
  s = []
  for word in words
    if length(word) > 4
      middle = join(shuffle(split(word[2:end-1], "")), "")
      push!(s, string(word[1]) * middle * string(word[end]))
    else
      push!(s, word)
    end
  end
  return join(s, " ")
end

sentence = "I couldn't believe that I could actually understand
 what I was reading : the phenomenal power of the human mind ."

words = split(sentence, " ")
println(typoglycemia(words))

# I cdon'ult bvlieee that I cluod actullay udnaetsdrn
# what I was rdiaeng : the pnnaoehmel peowr of the huamn mind .


text:
export data

abstract type TensorStorage end

#
# Generic ITensor storage functions
#

storage_randn!(S::TensorStorage) = randn!(data(S))
storage_norm(S::TensorStorage) = norm(data(S))
storage_conj(S::T) where {T<:TensorStorage}= T(conj(data(S)))

#
# This high level part of contract is generic to different
# storage types
#
# TODO: make this storage_contract!(), where C is pre-allocated.
#       This will allow for in-place multiplication
# TODO: optimize the contraction logic so C doesn't get permuted?
#
function storage_contract(Astore::TensorStorage,
                          Ais::IndexSet,
                          Bstore::TensorStorage,
                          Bis::IndexSet)
  if length(Ais)==0
    Cis = Bis
    Cstore = storage_scalar(Astore)*Bstore
  elseif length(Bis)==0
    Cis = Ais
    Cstore = storage_scalar(Bstore)*Astore
  else
    #TODO: check for special case when Ais and Bis are disjoint sets
    #I think we should do this analysis outside of storage_contract, at the ITensor level
    #(since it is universal for any storage type and just analyzes in indices)
    (Alabels,Blabels) = compute_contraction_labels(Ais,Bis)
    if is_outer(Alabels,Blabels)
      Cstore,Cis = storage_outer(Astore,Ais,Bstore,Bis)
    else
      (Cis,Clabels) = contract_inds(Ais,Alabels,Bis,Blabels)
      Cstore = _contract(Cis,Clabels,Astore,Ais,Alabels,Bstore,Bis,Blabels)
    end
  end
  return (Cis,Cstore)
end

# Generic outer function that handles proper
# storage promotion
# TODO: should this handle promotion with storage
# type switching?
# TODO: we should combine all of the storage_add!
# outer wrappers into a single call that promotes
# based on the storage type, i.e. promote_type(Dense,Diag) -> Dense
function storage_add!(Bstore::BT,
                      Bis::IndexSet,
                      Astore::AT,
                      Ais::IndexSet,
                      x::Number = 1) where {BT<:TensorStorage,AT<:TensorStorage}
  NT = promote_type(AT,BT)
  if NT == BT
    _add!(Bstore,Bis,Astore,Ais, x)
    return Bstore
  end
  Nstore = storage_convert(NT,Bstore,Bis)
  _add!(Nstore,Bis,Astore,Ais, x)
  return Nstore
end


text:
function check(result::Some)::(Union{Nothing,T} where {T<:Any})

    # ok, so check, do we have an error object?
    # Yes: log the error if we have a logger, then throw the error.
    # No: return the result.value

    # Error case -
    if (isa(something(result), Exception) == true)

        # get the error object -
        error_object = result.value

        # get the error message as a String -
        error_message = sprint(showerror, error_object, catch_backtrace())
        @error(error_message)

        # throw -
        throw(result.value)
    end

    # default -
    return result.value
end


text:
using JuLIP, JuLIP.Potentials
using JuLIPMaterials: Vec3
import JuLIPMaterials
const CLE = JuLIPMaterials.CLE

using JuLIPMaterials.CLE: grad, onb3D, QSB
using JuLIPMaterials.Testing
using Einsum, LinearAlgebra, Printf, ForwardDiff

println("---------------------------------------------------------------")
println(" Testing the 3D Anisotropic Dislocation Solution Implementation")
println("---------------------------------------------------------------")

# Random Lam\'e Parameters
λ = 1.0 + rand()
μ = 1.0 + rand()
Ciso = CLE.isotropic_moduli(λ, μ)

# TODO: Should probably not be in this set of tests
println("Test basis orientation implemented correctly: ")
maxerr = 0.0
# Check vertical
a = [0.0,0.0,1.0]
b,c = onb3D(JuLIPMaterials.Vec3(a))
maxerr = max( maxerr, abs( (a×b)⋅c - 1.0 ) )
# Check random directions
for n = 1:10
   a = randvec3()
   a /= norm(a)
   b,c = onb3D(JuLIPMaterials.Vec3(a))
   global maxerr = max( maxerr, abs( (a×b)⋅c - 1.0 ) )
end
print("maxerr = $maxerr: ")
println(@test maxerr < 1e-12)

println("Test implementation of Q,S,B matrices using BBS formulae (3.6.20-21): ")
maxerr_1 = 0.0
maxerr_2 = 0.0
# Check vertical
a = [0.0,0.0,1.0]
b,c = onb3D(JuLIPMaterials.Vec3(a))
Q,S,B = QSB(Ciso,b,c,30)
maxerr_1 = max( maxerr_1, norm(4*π*B*Q+S*S + I, Inf) )
maxerr_2 = max( maxerr_2, norm(Q*S'+S*Q, Inf) )
# Check random directions
for n = 1:10
   a = randvec3()
   a /= norm(a)
   b,c = onb3D(JuLIPMaterials.Vec3(a))
   Q,S,B = QSB(Ciso,b,c,30)
   global maxerr_1 = max( maxerr_1, norm(4*π*B*Q+S*S + I, Inf) )
   global maxerr_2 = max( maxerr_2, norm(Q*S'+S*Q, Inf) )
end
print("maxerr_1 = $maxerr_1: ")
println(@test maxerr_1 < 1e-12)
println("maxerr_2 = $maxerr_2")
println(@test maxerr_2 < 1e-12)

# Test explicit edge formula vs anistropic one
println("Test agreement between anisotopic and isotropic edge dislocation: ")
b = Vec3(1.0,0.0,0.0)
t = Vec3(0.0,0.0,1.0)
u0 = CLE.Dislocation(b, t, Ciso, Nquad = 40)
uedge = CLE.IsoEdgeDislocation3D(λ, μ, b[1])
maxerr = 0.0
maxerr_g = 0.0
for n = 1:10
   x = randvec3()
   global maxerr = max( maxerr, norm(u0(x) - uedge(x), Inf) )
   global maxerr_g = max(maxerr_g, norm(grad(u0, x) - grad(uedge, x), Inf) )
end
print("maxerr = $maxerr, maxerr_g = $maxerr_g: ")
print(@test maxerr < 1e-12)
println(@test maxerr_g < 1e-12)

# Test explicit screw formula vs anistropic one
println("Test agreement between anisotopic and isotropic screw dislocation: ")
Ciso = CLE.isotropic_moduli(λ, μ)
b = Vec3(0.0,0.0,1.0)
t = Vec3(0.0,0.0,1.0)
u0 = CLE.Dislocation(b, t, Ciso, Nquad = 40)
uscrew = CLE.IsoScrewDislocation3D(b[3])
maxerr = 0.0
maxerr_g = 0.0
for n = 1:10
   x = randvec3()
   global maxerr = max( maxerr, norm(u0(x) - uscrew(x), Inf) )
   global maxerr_g = max(maxerr_g, norm(grad(u0, x) - grad(uscrew, x), Inf) )
end
print("maxerr = $maxerr, maxerr_g = $maxerr_g: ")
print(@test maxerr < 1e-12)
println(@test maxerr_g < 1e-12)

# Generate random elastic moduli
Crand = randmoduli()
# Generate a random dislocation
b = randvec3();
b /= norm(b);
t = randvec3();
t /= norm(t);

# Test gradient implementation matches displacement implementation
for (Disl, id, C) in [ (CLE.IsoEdgeDislocation3D(λ, μ, 1.0), "IsoEdgeDislocation3D", Ciso),
         (CLE.IsoScrewDislocation3D(1.0), "IsoScrewDislocation3D", Ciso),
         (CLE.Dislocation(b,t,Crand, Nquad = 40), "Dislocation(30)", Crand) ]
   print("u = $id : test that ∇u is consistent with u: ")
   maxerr = 0.0
   for n = 1:10
      a, x = randvec3(), randvec3()
      u = x_ -> Disl(x_)
      ∂u = x_ -> grad(Disl, x_) * a
      ∂uad = x_ -> ForwardDiff.jacobian(Disl, x_) * a
      maxerr = max( maxerr, norm(∂u(x) - ∂uad(x), Inf) )
   end
   print("maxerr = $maxerr: ")
   println(@test maxerr < 1e-12)
end

# Test PDE is solved in all implementations
for (Disl, id, C) in [(CLE.IsoEdgeDislocation3D(λ, μ, 1.0), "IsoEdgeDislocation3D", Ciso),
         (CLE.IsoScrewDislocation3D(1.0), "IsoScrewDislocation3D", Ciso),
         (CLE.Dislocation(b,t,Crand, Nquad = 40), "Dislocation(30)", Crand) ]
   println("u = $id: test that u satisfies the PDE: ")
   maxerr = 0.0
   for n = 1:10
      x = Vec3(randvec3())
      u = x_ -> Disl(x_)
      maxerr = max( norm(cleforce(x, u, C), Inf), maxerr )
   end
   print("maxerr = $maxerr: ")
   println(@test maxerr < 1e-10)
end

# Test Burgers vector for edge dislocation implementation
Disl = CLE.IsoEdgeDislocation3D(λ, μ, 1.0)
println("u = IsoEdgeDislocation3D: test Burgers vector: ")
err = 0.0
# Integrate around loop
n = 30;
II = zeros(3)
DuTau = zeros(3)
for ω in range(0.0, step=pi/n, length=2*n)
   x = [cos(ω),sin(ω),0.0]
   tau = [-sin(ω),cos(ω),0.0]
   global II += grad(Disl,x) * tau
end
II = II*pi/n
maxerr = norm(II-[1.0,0.0,0.0])
print("maxerr = $maxerr: ")
println(@test maxerr < 1e-12)

# Test Burgers vector for screw dislocation implementation
Disl = CLE.IsoScrewDislocation3D(1.0)
println("u = IsoScrewDislocation3D: test Burgers vector: ")
err = 0.0
# Integrate around loop
n = 30;
II = zeros(3)
DuTau = zeros(3)
for ω in range(0.0, step=pi/n, length=2*n)
   x = [cos(ω),sin(ω),0.0]
   tau = [-sin(ω),cos(ω),0.0]
   global II += grad(Disl,x) * tau
end
II = II*pi/n
maxerr = norm(II-[0.0,0.0,1.0])
print("maxerr = $maxerr: ")
println(@test maxerr < 1e-12)

# Test Burgers vector for arbitrary anisotropic implementation
Disl = CLE.Dislocation(b,[0.0,0.0,1.0],Crand, Nquad = 30)
println("u = Dislocation(30): test Burgers vector: ")
err = 0.0
# Integrate around loop
n = 30;
II = zeros(3)
DuTau = zeros(3)
for ω in range(0.0, step=pi/n, length=2*n)
   x = [cos(ω),sin(ω),0.0]
   tau = [-sin(ω),cos(ω),0.0]
   global II += grad(Disl,x) * tau
end
II = II*pi/n
maxerr = norm(II-b)
print("maxerr = $maxerr: ")
println(@test maxerr < 1e-12)

# Test convergence with increasing quadrature points
println("Convergence of u with random C (please test this visually!): ")
println(" nquad |    err    |   err_g")
println("-------|-----------|-----------")
C = Crand
xtest = [ randvec3() for n = 1:10 ]
u0 = CLE.Dislocation(b,t,C, Nquad = 40)
for nquad in [2, 4, 6, 8, 10, 12, 14]
   u = CLE.Dislocation(b,t,C, Nquad = nquad)
   err = maximum( norm(u0(x) - u(x), Inf)   for x in xtest )
   err_g = maximum( norm(grad(u0, x) - grad(u, x), Inf)   for x in xtest )
   @printf("   %2d  | %.3e | %.3e  \n", nquad, err, err_g)
end


text:
# -- PRIVATE FUNCTIONS NOT EXPORTED -------------------------------------------------------------- #
function _extract_section(file_buffer_array::Array{String,1}, start_section_marker::String,
    end_section_marker::String)::Array{String,1}

    # initialize -
    section_buffer = String[]

    # find the SECTION START AND END -
    section_line_start = 1
    section_line_end = 1
    for (index,line) in enumerate(file_buffer_array)

        if (occursin(start_section_marker,line) == true)
            section_line_start = index
        elseif (occursin(end_section_marker,line) == true)
            section_line_end = index
        end
    end

    for line_index = (section_line_start+1):(section_line_end-1)
        line_item = file_buffer_array[line_index]
        push!(section_buffer,line_item)
    end

    # return -
    return section_buffer
end

function _extract_molecular_symbol_array(reaction_table::DataFrame)::Array{String,1}

    # strategy: compile a list of strings, in the order in which they are found from the L and R phrases.
    # Next, cut around the +, ane finally look for any * (which denote st coeff)

    # initialize -
    reaction_phrase_array = Array{String,1}()
    molecular_symbol_array = Array{String,1}()

    # what is the size of the reaction_table?
    (number_of_reaction, number_of_fields) = size(reaction_table)

    # process the L and R phrases -
    [push!(reaction_phrase_array,line) for line in reaction_table[!,:left_phrase]]
    [push!(reaction_phrase_array,line) for line in reaction_table[!,:right_phrase]]

    # ok, so now let's cut around the +'s -
    for phrase in reaction_phrase_array

        # compnents -
        if (occursin("+",phrase) == true)

            # ok, so we have some +'s
            component_array = split(phrase,"+")
            for component in component_array
                inner_component_array = split(component,"*")
                push!(molecular_symbol_array, string(last(inner_component_array)))
            end
        else
            # we have a single spcies, check for * -
            if (occursin("*",phrase) == true)

                # split -
                component_array = split(phrase,"*")

                # the symbol will be the last element -
                push!(molecular_symbol_array, string(last(component_array)))
            else
               # grab -
               push!(molecular_symbol_array, phrase)
            end
        end
    end

    # we want the array to be unique -
    unique!(molecular_symbol_array)

    # return the unique array -
    # NOTE: this is in random order
    return molecular_symbol_array
end

# ------------------------------------------------------------------------------------------------ #

# -- PUBLIC FUNCTIONS THAT ARE EXPORTED ---------------------------------------------------------- #
"""
    parse_vff_sequence_section(buffer::Array{String,1})::VLResult

Parse contents of the sequence section in the model file and extract amino acid and nucleotide sequences for protein and DNA, respectively. The sequence section is written between the tags #TXTL-SEQUENCE::START and #TXTL-SEQUENCE::STOP.

Input arguments:
`buffer::Array{String,1}` - 1-D array holding information extracted from the model file.

Output arguments:
`VLResult::VLResult` - concrete type holding protein and DNA sequence information.

"""
function parse_vff_sequence_section(buffer::Array{String,1})::VLResult

    # initialize -
    operation_type_array = Array{Symbol,1}()
    molecular_symbol_array = Array{String,1}()
    molecular_species_type_array = Array{Symbol,1}()
    sequence_array = Array{String,1}()

    # what are the alphabets for the sequences?
    dna_sequence_alphabet = ["a","t","g","c"]
    protein_sequence_alphabet = ["G","A","L","M","F","W","K","Q","E","S","P","V","I","C","Y","H","R","N","D","T"]

    try

        # extract the sequence section -
        sequence_section_buffer = _extract_section(buffer, "#TXTL-SEQUENCE::START", "#TXTL-SEQUENCE::STOP")
        if (isempty(sequence_section_buffer) == true)
            @warn "Hmmm. No SEQUENCE section was found. That's ok. We'll skip to the next section ..."
            return VLResult(nothing)
        end

        # make the buffer flat, and then write to tmp_file -
        flat_buffer = ""
        for line in sequence_section_buffer
            flat_buffer *= line

            if (occursin(";",line) == true)
                flat_buffer *= "\n"
            end
        end

        # remove the ; -
        flat_buffer = replace(flat_buffer,";"=>"")

        # ok, so now lets load the tmp file in CSV, and put into a
        df_tmp = CSV.read(IOBuffer(flat_buffer), DataFrame; header=false)

        # let's create a DataFrame w/the each peice of data in a col -
        [push!(molecular_symbol_array,line) for line in df_tmp[!,:Column1]]
        [push!(sequence_array,line) for line in df_tmp[!,:Column2]]

        # ok, so we need to do some additional logic for this next bit -
        # what type is the molecular symbol?
        for sequence_record in sequence_array

            # grab the first element -
            first_sequence_element = string(first(sequence_record))

            # if this in the dna alphabet?
            if (in(first_sequence_element, dna_sequence_alphabet) == true)

                push!(molecular_species_type_array,:DNA)

            elseif (in(first_sequence_element, protein_sequence_alphabet) == true)

                push!(molecular_species_type_array,:PROTEIN)
            else

                # we don't know what this is ... put type as missing for now
                push!(molecular_species_type_array, :UNKOWN)
            end
        end

        # Finally, add everything to new data frame -
        new_sequence_df = DataFrame(
            sequence_type=molecular_species_type_array,
            molecular_symbol=molecular_symbol_array,
            sequence=sequence_array)

        # return the new_df -
        return VLResult(new_sequence_df)
    catch error
        return VLResult(error)
    end
end

"""
    parse_vff_metabolic_section(buffer::Array{String,1})::VLResult

Parse contents of the metabolic section in the model file and extract all biochemical reactions and molecular species involved in the metabolic network. The metabolic section is written between the tags #METABOLISM::START and #METABOLISM::STOP.

Input arguments:
`buffer::Array{String,1}` - single array holding information extracted from the model file.

Output arguments:
`VLResult::VLResult` - concrete type holding information on reactions and species involved in metabolism.
"""
function parse_vff_metabolic_section(buffer::Array{String,1})::VLResult

    try
        # extract the metabolic section -
        metabolic_section_buffer = _extract_section(buffer, "#METABOLISM::START", "#METABOLISM::STOP")
        if (isempty(metabolic_section_buffer) == true)
            @warn "Hmmm. No METABOLISM section was found. That's ok. We'll skip to the next section ..."
            return VLResult(nothing)
        end

        # make the buffer flat, and then write to tmp_file -
        flat_buffer = ""
        for line in metabolic_section_buffer
            flat_buffer *= line
            flat_buffer *= "\n"
        end

        # ok, so now lets load the tmp file in CSV, and put into a
        df_tmp = CSV.read(IOBuffer(flat_buffer),DataFrame; header=false)

        # rename the cols -
        rename!(df_tmp,:Column1 => :reaction_tag)
        rename!(df_tmp,:Column2 => :ec_number_array)
        rename!(df_tmp,:Column3 => :left_phrase)
        rename!(df_tmp,:Column4 => :right_phrase)
        rename!(df_tmp,:Column5 => :is_reaction_reversible)

        # Next, let's get the list of molecular symbols -
        molecular_symbol_array = _extract_molecular_symbol_array(df_tmp)

        # we need to put the symbols in order -
        molecular_symbol_array = reorder_molecular_symbol_array(molecular_symbol_array; callback = nothing)

        # put the reactions in order -
        reaction_tag_array = df_tmp[!,:reaction_tag]
        reaction_tag_array = reorder_reaction_symbol_array(reaction_tag_array; callback = nothing)

        # ok, so let's build a NamedTuple and return to the caller =
        results_tuple = (reaction_table=df_tmp, molecular_symbol_array=molecular_symbol_array, reaction_tag_array=reaction_tag_array)

        # return the data frame -
        return VLResult(results_tuple)
    catch error
        return VLResult(error)
    end
end

"""
    parse_vff_grn_section(buffer::Array{String,1})::VLResult

Parse contents of the GRN section in the model file and extract information contained on gene regulation. The GRN section is written between the tags #GRN::START and #GRN::STOP.

Input arguments:
`buffer::Array{String,1}` - 1-D array holding information extracted from the model file.

Output arguments:
`VLResult::VLResult` - concrete type holding information on the gene regulatory network.
"""
function parse_vff_grn_section(buffer::Array{String,1})::VLResult

    # initialize -
    original_record_buffer_dictionary = Dict{Int64,Any}()

    try

        # get the grn section -
        grn_section_buffer = _extract_section(buffer, "#GRN::START", "#GRN::STOP")
        if (isempty(grn_section_buffer) == true)
            @warn "Hmmm. No GRN section was found. That's ok. We'll skip to the next section ..."
            return VLResult(nothing)
        end

        # initial tokenize by spaces -
        for (index,grn_record) in enumerate(grn_section_buffer)

            # split -
            fragment_array = split(grn_record, " ")

            # grab -
            original_record_buffer_dictionary[index] = string.(fragment_array)
        end

        # scan the original record buffer -> to produce the cannonical_reduced_array
        scan_result = minerva_scanner(original_record_buffer_dictionary, minerva_grn_scan_function)
        cannonical_reduced_array = check(scan_result)

        # for now - the cannonical_reduced_array -
        return VLResult(cannonical_reduced_array)
    catch error
        return VLResult(error)
    end
end

"""
    parse_vff_species_bounds_section(buffer::Array{String,1}, metabolic_results_tuple::NamedTuple)::VLResult

Parse contents of the model file and extract information on species bounds. The species bounds section is written between the tags #SPECIES_BOUNDS::START and #SPECIES_BOUNDS::STOP.

Input arguments:
`buffer::Array{String,1}` - 1-D array holding information extracted from the model file.

Output arguments:
`VLResult::VLResult` - concrete type holding information on species bounds.
"""
function parse_vff_species_bounds_section(buffer::Array{String,1})::VLResult

    # initialize -
    original_record_buffer_dictionary = Dict{Int64,Any}()

    try

        # extract the species bounds section -
        species_bounds_section_buffer = _extract_section(buffer, "#SPECIES_BOUNDS::START", "#SPECIES_BOUNDS::STOP")
        if (isempty(species_bounds_section_buffer) == true)
            @warn "Hmmm. No SPECIES_BOUNDS section was found. That's ok. We'll skip to the next section ..."
            return VLResult(nothing)
        end

        # initial tokenize by spaces -
        for (index, species_bound_record) in enumerate(species_bounds_section_buffer)

            # split -
            fragment_array = split(species_bound_record, " ")

            # grab -
            original_record_buffer_dictionary[index] = string.(fragment_array)
        end

        # scan the original record buffer -> to produce the cannonical_reduced_array
        scan_result = minerva_scanner(original_record_buffer_dictionary, minerva_bound_type_assignment_scan_function)
        cannonical_reduced_array = check(scan_result)

        # parse the cannonical records -

        # for now - the cannonical_reduced_array -
        return VLResult(cannonical_reduced_array)
    catch error
        return VLResult(error)
    end
end

"""
    parse_vff_bio_types_section(buffer::Array{String,1})::VLResult

Parse contents of the model file and tokenize biological types. The bio-types section is written between the tags #BIO-TYPE-PREFIXES::START and #BIO-TYPE-PREFIXES::STOP.

Input arguments:
`buffer::Array{String,1}` - single array holding information extracted from the model file.

Output arguments:
`VLResult::VLResult` - concrete type holding information on biological types.
"""
function parse_vff_bio_types_section(buffer::Array{String,1})::VLResult

    # initialize -
    original_record_buffer_dictionary = Dict{Int64,Any}()

    try

        # load the types section buffer -
        types_section_buffer = _extract_section(buffer, "BIO-TYPE-PREFIXES::START", "BIO-TYPE-PREFIXES::STOP")
        if (isempty(types_section_buffer) == true)
            @warn "Hmmm. No BIO-TYPES-PREFIXES section was found. That's ok. We'll skip to the next section ..."
            return VLResult(nothing)
        end

        # ok, lets see how this tokenizes -
        for (index,type_record) in enumerate(types_section_buffer)

            # tokenize -
            token_array = tokenize(type_record)

            # grab -
            original_record_buffer_dictionary[index] = token_array
        end

         # scan the original record buffer -> to produce the cannonical_reduced_array
         scan_result = minerva_scanner(original_record_buffer_dictionary, biological_type_assignment_scan_function)
         cannonical_reduced_array = check(scan_result)

         # for now - the cannonical_reduced_array -
         return VLResult(cannonical_reduced_array)
    catch error
        return VLResult(error)
    end
end

"""
    parse_vff_model_document(model::VLAbstractModelObject)::VLResult

Parse contents of the model file, extract information from various sections of the file and construct a rudimentary data dictionary.

Input arguments:
`model::VLAbstractModelObject` - abstract Julia data object holding information for generating model code.

Output arguments:
`VLResult::VLResult` - concrete type holding an intermediate data dictionary construction.
"""
function parse_vff_model_document(model::VLAbstractModelObject)::VLResult

    # initialize -
    intermediate_representation_dictionary = Dict{String,Any}()

    try

         # get the path to the vff model -
        vff_model_file_path = model.path_to_model_file

        # load the vff buffer -
        vff_file_buffer = read_model_document(vff_model_file_path)

        # -- TYPES SECTION ------------------------------------------------------------------------------- #
        types_parse_result = parse_vff_bio_types_section(vff_file_buffer)

        # ------------------------------------------------------------------------------------------------ #

        # -- SEQ SECTION --------------------------------------------------------------------------------- #
        result = parse_vff_sequence_section(vff_file_buffer)
        if (isa(result.value, Exception) == true)
            throw(result.value)
        end
        sequence_data_table = result.value
        # ------------------------------------------------------------------------------------------------ #

        # -- METABOLISM SECTION -------------------------------------------------------------------------- #
        result = parse_vff_metabolic_section(vff_file_buffer)
        if (isa(result.value,Exception) == true)
            throw(result.value)
        end
        metabolic_section_results_tuple = result.value
        # ------------------------------------------------------------------------------------------------ #

        # -- GRN SECTION --------------------------------------------------------------------------------- #
        grn_section_result = parse_vff_grn_section(vff_file_buffer)
        grn_scanned_dictionary = check(grn_section_result)
        # ------------------------------------------------------------------------------------------------ #

        # -- SPECIES BOUNDS SECTION ---------------------------------------------------------------------- #
        species_bound_result = parse_vff_species_bounds_section(vff_file_buffer)
        species_scanner_dictionary = check(species_bound_result)
        # ------------------------------------------------------------------------------------------------ #

        # ok, put stuff in the IR dictionary -
        intermediate_representation_dictionary[ir_master_reaction_table_key] = metabolic_section_results_tuple.reaction_table
        intermediate_representation_dictionary[ir_list_of_molecular_species_key] = metabolic_section_results_tuple.molecular_symbol_array
        intermediate_representation_dictionary[ir_list_of_reaction_tags_key] = metabolic_section_results_tuple.reaction_tag_array
        intermediate_representation_dictionary[ir_master_species_bounds_table_key] = species_scanner_dictionary
        intermediate_representation_dictionary[ir_sequence_section_table_key] = sequence_data_table

        # return -
        return VLResult(intermediate_representation_dictionary)
    catch error
        throw(error)
    end
end
# ------------------------------------------------------------------------------------------------ #

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
tests = ["direct.jl",
         "indirect.jl",
         "mpb_linear.jl",
         "options.jl"]

println("Running tests:")

for curtest in tests
    println(" Test: $(curtest)")
    include(curtest)
end

import SCS
include(joinpath(Pkg.dir("MathProgBase"),"test","conicinterface.jl"))
coniclineartest(SCS.SCSSolver(), duals=true, tol=1e-2)
conicSOCtest(SCS.SCSSolver(), duals=true, tol=1e-2)
conicEXPtest(SCS.SCSSolver(), duals=true, tol=1e-2)
conicSDPtest(SCS.SCSSolver(), duals=true, tol=1e-2)


text:
module Source

using SparseArrays

"""
Determine the indices and values of the source discretization stencil that
spans `p` cells.

Input arguments:
x : List of coordinates
xs : Source location
p : Polynomial degree.
s : Number of smoothness conditions.

Return values:
stencil : An array of the stencil indices that maps to the `x` coordinates.
coef : An array of stencil coefficients (values).
alpha : The signed distance from the source location to the nearest grid point
"""
function source_smoothness(x, xs, p, s)

        idx = argmin(abs.(x .- xs))

        # Determine signed distance from source location to nearest point.
        alpha = x[idx] - xs

        shift::Int64 = 0
        if alpha < 0
                shift = 0
        else
                shift = -1
        end

        if p == 0
                coef = 1
                stencil = idx
                return
        end

        if mod(p,1) == 1
                error("Polynomial degree must be odd")
        end

        left = round((p-1)/2) + round(s/2);
        right = p + s - left - 1;

        idx_l::Int64 = idx - left + shift;
        idx_r::Int64 = idx + right + shift + 1;

        adjust::Int64 = 0;
        if idx_l < 1
                adjust = 1 - idx_l;
        end
        if idx_r > length(x)
                adjust = length(x) - idx_r;
        end

        stencil = collect((adjust + idx_l):(adjust + idx_r));

        # Moment conditions
        e = ones(size(stencil))
        xl = x[stencil] .- x[idx];
        M = zeros(p+1, size(xl, 1))
        for i in 1:p+1
        for j in 1:size(xl, 1)
                M[i,j] = xl[j] ^ (i - 1)
        end
        end
        b1 = (xs - x[idx]).^(0:p);

        if s > 0
                # Smoothness conditions
                S = zeros(s,size(xl, 1))
                for i in 1:s
                for j in 1:size(xl, 1)
                        S[i,j] = (-1)^j * xl[j] ^ (i - 1)
                end
                end
                A = [M; S]
                b = [b1; zeros(s,1)]
        else
                A = M
                b = b1
        end
        coef = A\b;
return stencil, coef, alpha
end

function source_discretize_1d(r1_s, p, s, r1, h1)
        d = spzeros(length(r1))
        stencil1, coef1, alpha1 = source_smoothness(r1, r1_s, p, s)
        for i=1:length(stencil1)
                d[stencil1[i]] = coef1[i] / h1
        end

        return d
end
"""
d = source_disc(xs,ys, h, gr, J)
Discretize point source using nearest grid point

Arguments:
xs, ys : Source location
h : Grid spacing in parameter space
gr : grid object that contains X,Y meshgrids.
J : Jacobian

Returns
A sparse vector that is zero everywhere except for one index, which is set to
a weight d = 1/(cell-area).
"""
function source_discretize_2d(r1_s, r2_s, p, s, r1, r2, h1, h2)
        d = spzeros(length(r1)*length(r2))
        stencil1, coef1, alpha1 = source_smoothness(r1, r1_s, p, s)
        stencil2, coef2, alpha2 = source_smoothness(r2, r2_s, p, s)

        n2 = length(r2)
        reci_area = 1 / (h1 * h2)
        for i=1:length(stencil1)
        for j=1:length(stencil2)
          d[stencil2[j] + n2 * (stencil1[i] - 1)] = coef1[i] * coef2[j] * reci_area
        end
        end

        return d
end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
# floatfuncs2arg = Dict{Symbol,Any}(
# :copysign => :todo,                      # floatfuncs,operators
# :flipsign => :todo,                      # floatfuncs,operators
# )

# TODO:

# eval
# copysign
# flipsign
# signbit
# maxintfloat
# isinteger
# num2hex
# hex2num
# round
# $(Expr(:$, :f)): Not a symbol
# _signif_og: Not exported
# signif
# isapprox
# ≉
# rtoldefault: Not exported
# fma_libm: Not exported
# fma_llvm: Not exported
# fma


text:
#Type definition for the aerodynamic database, including a function to
# add new contributions per contribution. (TBD: Including a remove function)
export AbstractDataBase

export AeroDataBase

abstract type AbstractDataBase end

mutable struct AeroDataBase <: AbstractDataBase
  Variable::String
  ListOfInputs::Array{String}
  ListOfContributions::Array{Any}
  ListOfInputValues::Array{Any}
  ListOfInputOrder::Array{Any}
  AddContribution::Function

  function AeroDataBase(varkey::String)
    this = new()
    this.Variable = varkey
    this.ListOfInputs = []
    this.ListOfContributions = []
    this.ListOfInputValues = []
    this.ListOfInputOrder = []

    #Add method to add contributions to this specific database
    this.AddContribution = function (Contr::Array{Float64},varargs::Tuple)
      #check if correct number of dimensions have been assigned
      nvarargs = length(varargs)
      ndims(Contr) != nvarargs ? error("Incorrect number of inputs") : nothing

      #Add Contr to ListOfContributions
      this.ListOfContributions = [this.ListOfContributions..., Contr]
      NewInput = Vector(nvarargs)
      NewInputOrder = Vector(nvarargs)

      for i=1:length(varargs)
        if IsContainedInStringArray(this.ListOfInputs,varargs[i][1])
          NewInput[i] = varargs[i][2:end]
          NewInputOrder[i] = IdInStringArray(this.ListOfInputs,varargs[i][1])
        else
          this.ListOfInputs = [this.ListOfInputs...,varargs[i][1]]
          NewInput[i] = varargs[i][2:end]
          NewInputOrder[i] = length(this.ListOfInputs)
        end
      end
      this.ListOfInputValues = push!(this.ListOfInputValues, NewInput)
      this.ListOfInputOrder = push!(this.ListOfInputOrder, NewInputOrder)
      return nothing
    end #this.AddContribution
    return this
  end #Constructor
end #typedef

#Example
#=
Cm1 = rand(4,7)
Cm2 = rand(5,7)
x = ["alpha", 2,3,4,5]
y = ["beta",1,5,2,7,9,4,9]
z = ["sigma",1,2,3,4,5]
Cmdb = AeroDataBase("Cm")
Cmdb.AddContribution(Cm1,(x,y))
Cmdb.AddContribution(Cm2,(z,y))
=#


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module GaussianMixtureTest

using Distributions, StatsBase, StatsFuns, Yeppp

#import Yeppp: add!, exp!, log!
import StatsBase: RealArray, RealVector, RealArray, IntegerArray, IntegerVector, IntegerMatrix, IntUnitRange

export gmm, gmmrepeat, asymptoticdistribution, kstest

include("estimate.jl")
include("arithmetic.jl")
include("test.jl")
end # module


text:
import Kuber, HTTP, DataStructures, JSON, Dates

include("client.jl")
include("resource.jl")
include("controllers/Controllers.jl")

"""
    objectListWatcher(client::K8sClient)

A objectListWatcher will keep a cache of the current state of the object CRDs in the cluster by watching
the resource and periodically listing it.
"""
mutable struct ListWatcher
    client::K8sClient
    chan::Channel
    resource::Resource
    objects::Dict{String,Any}
end

function listwatch(listwatcher::ListWatcher)
    @info "listing objects"
    objects = list_objects(listwatcher)

    # check if object in in current cache otherwise send to chan
    for objectid in keys(objects)
        object = objects[objectid]
        existing = get(listwatcher.objects, object["metadata"]["uid"], Dict{String, Any}())
        if length(existing) == 0
            @info "sending add/update to queue"
            put!(listwatcher.chan, object)
            listwatcher.objects[objectid] = object
        end
    end
    # check if any current objects no longer exist, if so send to chan
    for objectid in keys(listwatcher.objects)
        object = listwatcher.objects[objectid]
        existing = get(objects, object["metadata"]["uid"], Dict{String, Any}())
        if length(existing) == 0
            @info "sending deletion to queue"
            put!(listwatcher.chan, object)
            delete!(listwatcher.objects, objectid)
        end
    end
end

function loop(listwatcher::ListWatcher, interval::Dates.Second)
    @info "starting loop..."
    while true
        listwatch(listwatcher)
        sleep(Dates.value(interval))
    end
end

function list_objects(listwatcher::ListWatcher)::Dict{String, Any}
    resp = list(listwatcher.client, listwatcher.resource)
    str = String(resp.body)
    jobj = JSON.Parser.parse(str)
    currentObjects = Dict{String, Any}()
    for item in jobj["items"]
        uid = item["metadata"]["uid"]
        currentObjects[uid] = item
    end
    return currentObjects
end

function watch_object(listwatcher::ListWatcher)
    # TODO
end

# TODO: make this a more generic crd informer
"""
    Informer(listwatcher::ListWatcher)

An Informer will inform on what changes have occurred with the Resources, and apply the appropriate
handler functions.
"""
mutable struct Informer
    listwatcher::ListWatcher
    chan::Channel
    resource::Resource
    controller
end

function run(informer::Informer)
    @info "running async"
    @async loop(informer.listwatcher, Dates.Second(2))
    while true
        @info "taking from chan"
        data = take!(informer.chan)
        @debug "received data: " data
        dispatch(informer.listwatcher, informer.resource, informer.controller, data)
    end
end

function dispatch(listwatcher::ListWatcher, resource::Resource, controller, object)
    found = get(listwatcher.client, resource, object["metadata"]["name"], object["metadata"]["namespace"])
    # object is missing
    if length(found) == 0
        @info "deleting object: " object
        delete(controller, object)
        return
    end
    # object is present, check if ids match
    if object["metadata"]["resourceVersion"] != found["metadata"]["resourceVersion"]
        @info "updating object: " object
        update(controller, object)
        return
    end
    # resource versions and names match, its an add
    @info "adding object: " object
    add(controller, object)
    return
end

max_stars_repo_name: pbarker/kubemeta

max_issues_repo_name: pbarker/kubemeta

max_forks_repo_name: pbarker/kubemeta


text:
export CuboidRoom

"""
`cuboid room object`

Returns a `CuboidRoom` object containing the geometrical information of the room.

# Usage

* CuboidRoom(Lx::Float64,Ly::Float64,Lz::Float64,ξ::Array{Float64,1}, env::FDTDEnv)
* CuboidRoom(Nx::Int64,Ny::Int64,Nz::Int64,ξ::Array{Float64,1}, env::FDTDEnv)
  * `ξ::Array{Float64}`: must be 6 element Array containg acoustic impedance of the 6 walls

"""
immutable CuboidRoom <: AbstractGeometry
	Lx::Float64   #x dimension
	Ly::Float64   #y dimension
	Lz::Float64   #z dimension

	Nx::Int64   #x-axis samples
	Ny::Int64   #y-axis samples
	Nz::Int64   #z-axis samples

	ξ::Array{Float64,1}   #Impedance
	G::Array{UInt8,3}     #NodeProperty
	env::FDTDEnv          #FDTD acoustic env.
	function CuboidRoom(Lx,Ly,Lz,Nx,Ny,Nz,ξ,env)
		if(any([Lx;Ly;Lz].< 0)) error("room dimensions L should be positive") end
		if(length(ξ)!= 6) error("length(ξ) must be 6") end
		g = ones(Bool,Nx+2,Ny+2,Nz+2)
		g[1,:,:] = g[end,:,:] = false
		g[:,1,:] = g[:,end,:] = false
		g[:,:,1] = g[:,:,end] = false
		G = AcFDTD.getNodeProperty(g)
		new(Lx,Ly,Lz,Nx,Ny,Nz,ξ,G,env)
	end

end

function CuboidRoom(Lx::Float64,Ly::Float64,Lz::Float64,
		    ξ::Array{Float64,1}, env::FDTDEnv)

	Nx = round(Int64,Lx/env.X)
	Ny = round(Int64,Ly/env.X)
	Nz = round(Int64,Lz/env.X)

	CuboidRoom(env.X*Nx,env.X*Ny,env.X*Nz,Nx,Ny,Nz,ξ,env)
end

function CuboidRoom(Nx::Int64,Ny::Int64,Nz::Int64,
		    ξ::Array{Float64,1}, env::FDTDEnv)
	Lx = Nx*env.X
	Ly = Ny*env.X
	Lz = Nz*env.X

	CuboidRoom(Lx,Ly,Lz,Nx,Ny,Nz,ξ,env)
end

fun_name(f::CuboidRoom) = "cuboid room"
fun_dim(f::CuboidRoom) = @sprintf("Lx x Ly x Lz = %.2f x %.2f x %.2f m^3", f.Lx, f.Ly, f.Lz)
fun_Nxyz(f::CuboidRoom) = @sprintf("Nx x Ny x Nz = %d x %d x %d", f.Nx, f.Ny, f.Nz)
fun_ξ(f::CuboidRoom) =
@sprintf("frequency independent   \n\n |  ξx1 |  ξx2 |  ξy1 |  ξy2 |  ξz1 | ξz2  | \n | % 02.1f| % 02.1f| % 02.1f| % 02.1f| % 02.1f| % 02.1f|\n" , f.ξ[1],f.ξ[2],f.ξ[3],f.ξ[4],f.ξ[5],f.ξ[6])
fun_X(f::CuboidRoom) = @sprintf(" %.2f cm", f.env.X*100)
fun_Fs(f::CuboidRoom) = @sprintf(" %.2f kHz", f.env.Fs/1000)

max_stars_repo_name: nantonel/JuFdtd

max_issues_repo_name: nantonel/JuFdtd

max_forks_repo_name: nantonel/JuFdtd


text:
using IterTools, LinearAlgebra
import LinearAlgebra: AbstractMatrix
import Base: +, *, ==, size, eachindex, print
import Distances: pairwise
export CrossKernel, Kernel, cov, xcov, EQ, RQ, Linear, Poly, Noise, Wiener, WienerVelocity,
    Exponential, ConstantKernel, isstationary, ZeroKernel
using IterTools: product

############################# Define CrossKernels and Kernels ##############################

abstract type CrossKernel end
abstract type Kernel <: CrossKernel end

@inline isstationary(::Type{<:CrossKernel}) = false
@inline isstationary(x::CrossKernel) = isstationary(typeof(x))

# Some fallback definitions.
size(::CrossKernel, N::Int) = (N ∈ (1, 2)) ? Inf : 1
size(k::CrossKernel) = (size(k, 1), size(k, 2))

eachindex(k::Kernel, N::Int) = eachindex(k)
size(k::Kernel) = (length(k), length(k))
size(k::Kernel, dim::Int) = size(k)[dim]
length(k::Kernel) = Inf

"""
    AbstractMatrix(k::Kernel)

Convert `k` into an `AbstractMatrix`, if such a representation exists.
"""
function AbstractMatrix(k::Kernel)
    @assert isfinite(size(k, 1))
    return pairwise(k, eachindex(k, 1))
end

"""
    AbstractMatrix(k::CrossKernel)

Convert `k` into an `AbstractMatrix`, if such a representation exists.
"""
function AbstractMatrix(k::CrossKernel)
    @assert isfinite(size(k, 1))
    @assert isfinite(size(k, 2))
    return pairwise(k, eachindex(k, 1), eachindex(k, 2))
end

###################### `map` and `pairwise` fallback implementations #######################

# Unary map / _map
_map_fallback(k::CrossKernel, X::AV) = [k(x, x) for x in X]
_map_fallback(k::Kernel, X::AV) = [k(x) for x in X]
_map(k::CrossKernel, X::AV) = _map_fallback(k, X)
map(k::CrossKernel, X::BlockData) = BlockVector([map(k, x) for x in blocks(X)])
map(k::CrossKernel, X::AV) = _map(k, X)

# Binary map / _map
_map_fallback(k::CrossKernel, X::AV, X′::AV) = [k(x, x′) for (x, x′) in zip(X, X′)]
_map(k::CrossKernel, X::AV, X′::AV) = _map_fallback(k, X, X′)
function map(k::CrossKernel, X::BlockData, X′::BlockData)
    return BlockVector([map(k, x, x′) for (x, x′) in zip(blocks(X), blocks(X′))])
end
map(k::CrossKernel, X::AV, X′::AV) = _map(k, X, X′)

# Unary pairwise / _pairwise
_pairwise(k::Kernel, X::AV) = _pairwise(k, X, X)
pairwise(k::Kernel, X::AV) = LazyPDMat(_pairwise(k, X))
function pairwise(k::Kernel, X::BlockData)
    Σ = BlockMatrix([pairwise(k, x, x′) for x in blocks(X), x′ in blocks(X)])
    return LazyPDMat(Symmetric(Σ))
end
pairwise(k::CrossKernel, X::BlockData) = pairwise(k, X, X)
pairwise(k::CrossKernel, X::AV) = _pairwise(k, X)

# Binary pairwise / _pairwise
function _pairwise_fallback(k::CrossKernel, X::AV, X′::AV)
    return [k(X[p], X′[q]) for p in eachindex(X), q in eachindex(X′)]
end
_pairwise(k::CrossKernel, X::AV, X′::AV) = _pairwise_fallback(k, X, X′)
_pairwise(k::CrossKernel, X::AV) = _pairwise(k, X, X)

pairwise(k::CrossKernel, X::AV, X′::AV) = _pairwise(k, X, X′)
function pairwise(k::CrossKernel, X::BlockData, X′::BlockData)
    return BlockMatrix([pairwise(k, x, x′) for x in blocks(X), x′ in blocks(X′)])
end
pairwise(k::CrossKernel, X::BlockData, X′::AV) = pairwise(k, X, BlockData([X′]))
pairwise(k::CrossKernel, X::AV, X′::BlockData) = pairwise(k, BlockData([X]), X′)

# Sugar for `eachindex` things.
for op in [:map, :pairwise]
    @eval begin
        $op(k::CrossKernel, ::Colon) = $op(k, eachindex(k))
        $op(k::CrossKernel, ::Colon, ::Colon) = $op(k, :)
        $op(k::CrossKernel, ::Colon, X′::AV) = $op(k, eachindex(k, 1), X′)
        $op(k::CrossKernel, X::AV, ::Colon) = $op(k, X, eachindex(k, 2))
    end
end

# Optimisation for Toeplitz covariance matrices.
function pairwise(k::Kernel, x::StepRangeLen{<:Real})
    if isstationary(k)
        return LazyPDMat(SymmetricToeplitz(map(k, x, Fill(x[1], length(x)))))
    else
        return LazyPDMat(_pairwise(k, x))
    end
end
function pairwise(k::CrossKernel, x::StepRangeLen{<:Real}, x′::StepRangeLen{<:Real})
    if isstationary(k) && x.step == x′.step
        return Toeplitz(
            map(k, x, Fill(x′[1], length(x))),
            map(k, Fill(x[1], length(x′)), x′),
        )
    else
        return _pairwise(k, x, x′)
    end
end

################################ Define some basic kernels #################################

# An error that Kernels may throw if `eachindex` is undefined.
function eachindex_err(::T) where T<:Kernel
    throw(ArgumentError("`eachindex` undefined for kernel of type $T."))
end

"""
    ZeroKernel <: Kernel

A rank 0 `Kernel` that always returns zero.
"""
struct ZeroKernel{T<:Real} <: Kernel end
(::ZeroKernel{T})(x, x′) where T = zero(T)
(::ZeroKernel{T})(x) where T = zero(T)
isstationary(::Type{<:ZeroKernel}) = true
@inline _map(::ZeroKernel{T}, X::AV, X′::AV) where T = Zeros{T}(length(X))
@inline function _pairwise(::ZeroKernel{T}, X::AV, X′::AV) where T
    return Zeros{T}(length(X), length(X′))
end
@inline ==(::ZeroKernel{<:Any}, ::ZeroKernel{<:Any}) = true
@inline eachindex(k::ZeroKernel) = eachindex_err(k)

const ZK = ZeroKernel{Float64}
zero(k::Kernel) = length(k) < Inf ? FiniteZeroKernel(eachindex(k)) : ZK()

"""
    ConstantKernel{T<:Real} <: Kernel

A rank 1 constant `Kernel`. Useful for consistency when creating composite Kernels,
but (almost certainly) shouldn't be used as a base `Kernel`.
"""
struct ConstantKernel{T<:Real} <: Kernel
    c::T
end
(k::ConstantKernel)(x::T, x′::T) where T = k.c
(k::ConstantKernel)(x) = k.c
isstationary(::Type{<:ConstantKernel}) = true
_map(k::ConstantKernel, X::AV, ::AV) = Fill(k.c, length(X))
_pairwise(k::ConstantKernel, X::AV, X′::AV) = Fill(k.c, length(X), length(X′))
==(k::ConstantKernel, k′::ConstantKernel) = k.c == k′.c
@inline eachindex(k::ConstantKernel) = eachindex_err(k)

# ConstantKernel-specific optimisations.
+(k::ConstantKernel, k′::ConstantKernel) = ConstantKernel(k.c + k′.c)
*(k::ConstantKernel, k′::ConstantKernel) = ConstantKernel(k.c * k′.c)

"""
    EQ <: Kernel

The standardised Exponentiated Quadratic kernel with no free parameters.
"""
struct EQ <: Kernel end
isstationary(::Type{<:EQ}) = true
(::EQ)(x, x′) = exp(-0.5 * sqeuclidean(x, x′))
(::EQ)(x::Real, x′::Real) = exp(-0.5 * (x - x′)^2)
(::EQ)(x::T) where T = one(Float64)
_pairwise(::EQ, X::ColsAreObs) = exp.(-0.5 .* pairwise(SqEuclidean(), X.X))
_pairwise(::EQ, X::ColsAreObs, X′::ColsAreObs) = exp.(-0.5 .* pairwise(SqEuclidean(), X.X, X′.X))
@inline eachindex(k::EQ) = eachindex_err(k)

"""
    PerEQ{Tp<:Real}

The usual periodic kernel derived by mapping the input domain onto a circle.
"""
struct PerEQ{Tp<:Real} <: Kernel
    p::Tp
end
isstationary(::Type{<:PerEQ}) = true
(k::PerEQ)(x::Real, x′::Real) = exp(-2 * sin(π * abs(x - x′) / k.p)^2)
(k::PerEQ)(x::Real) = one(typeof(x))
@inline eachindex(k::PerEQ) = eachindex_err(k)

"""
    Exponential <: Kernel

The standardised Exponential kernel.
"""
struct Exponential <: Kernel end
isstationary(::Type{<:Exponential}) = true
(::Exponential)(x::Real, x′::Real) = exp(-abs(x - x′))
(::Exponential)(x) = one(Float64)
@inline eachindex(k::Exponential) = eachindex_err(k)

# """
#     RQ{T<:Real} <: Kernel

# The standardised Rational Quadratic. `RQ(α)` creates an `RQ` `Kernel{Stationary}` whose
# kurtosis is `α`.
# """
# struct RQ{T<:Real} <: Kernel
#     α::T
# end
# @inline (k::RQ)(x::Real, y::Real) = (1 + 0.5 * abs2(x - y) / k.α)^(-k.α)
# ==(a::RQ, b::RQ) = a.α == b.α
# isstationary(::Type{<:RQ}) = true
# show(io::IO, k::RQ) = show(io, "RQ($(k.α))")

"""
    Linear{T<:Real} <: Kernel

Standardised linear kernel. `Linear(c)` creates a `Linear` `Kernel{NonStationary}` whose
intercept is `c`.
"""
struct Linear{T<:Union{Real, Vector{<:Real}}} <: Kernel
    c::T
end
==(a::Linear, b::Linear) = a.c == b.c
(k::Linear)(x, x′) = dot(x .- k.c, x′ .- k.c)
(k::Linear)(x) = sum(abs2, x .- k.c)
@inline eachindex(k::Linear) = eachindex_err(k)

@inline _pairwise(k::Linear, x::AV) = _pairwise(k, ColsAreObs(x'))
@inline function _pairwise(k::Linear, x::AV, x′::AV)
    return _pairwise(k, ColsAreObs(x'), ColsAreObs(x′'))
end

function _pairwise(k::Linear, D::ColsAreObs)
    Δ = D.X .- k.c
    return Δ' * Δ
end
_pairwise(k::Linear, X::ColsAreObs, X′::ColsAreObs) = (X.X .- k.c)' * (X′.X .- k.c)

# """
#     Poly{Tσ<:Real} <: Kernel

# Standardised Polynomial kernel. `Poly(p, σ)` creates a `Poly`.
# """
# struct Poly{Tσ<:Real} <: Kernel
#     p::Int
#     σ::Tσ
# end
# @inline (k::Poly)(x::Real, x′::Real) = (x * x′ + k.σ)^k.p
# show(io::IO, k::Poly) = show(io, "Poly($(k.p))")

"""
    Noise{T<:Real} <: Kernel

A white-noise kernel with a single scalar parameter.
"""
struct Noise{T<:Real} <: Kernel
    σ²::T
end
isstationary(::Type{<:Noise}) = true
==(a::Noise, b::Noise) = a.σ² == b.σ²
(k::Noise)(x, x′) = x === x′ || x == x′ ? k.σ² : zero(k.σ²)
(k::Noise)(x) = k.σ²
_pairwise(k::Noise, X::AV) = Diagonal(Fill(k.σ², length(X)))
function _pairwise(k::Noise, X::AV, X′::AV)
    if X === X′
        return _pairwise(k, X)
    else
        return [view(X, p) == view(X′, q) ? k.σ² : 0
            for p in eachindex(X), q in eachindex(X′)]
    end
end
@inline eachindex(k::Noise) = eachindex_err(k)

# """
#     Wiener <: Kernel

# The standardised stationary Wiener-process kernel.
# """
# struct Wiener <: Kernel end
# @inline (::Wiener)(x::Real, x′::Real) = min(x, x′)
# cov(::Wiener, X::AM, X′::AM) =
# show(io::IO, ::Wiener) = show(io, "Wiener")

# """
#     WienerVelocity <: Kernel

# The standardised WienerVelocity kernel.
# """
# struct WienerVelocity <: Kernel end
# @inline (::WienerVelocity)(x::Real, x′::Real) =
#     min(x, x′)^3 / 3 + abs(x - x′) * min(x, x′)^2 / 2
# show(io::IO, ::WienerVelocity) = show(io, "WienerVelocity")

"""
    EmpiricalKernel <: Kernel

A finite-dimensional kernel defined in terms of a PSD matrix `Σ`.
"""
struct EmpiricalKernel{T<:LazyPDMat} <: Kernel
    Σ::T
end
@inline (k::EmpiricalKernel)(q::Int, q′::Int) = k.Σ[q, q′]
@inline (k::EmpiricalKernel)(q::Int) = k.Σ[q, q]
@inline length(k::EmpiricalKernel) = size(k.Σ, 1)
eachindex(k::EmpiricalKernel) = eachindex(k.Σ, 1)

function _pairwise(k::EmpiricalKernel, X::AV)
    if X == eachindex(k)
        return k.Σ
    else
        return k.Σ[X, X]
    end
end

function _pairwise(k::EmpiricalKernel, X::AV, X′::AV)
    if X == eachindex(k) && X′ == eachindex(k)
        return k.Σ
    else
        return k.Σ[X, X′]
    end
end
AbstractMatrix(k::EmpiricalKernel) = k.Σ

+(x::ZeroKernel, x′::ZeroKernel) = zero(x)
function +(k::CrossKernel, k′::CrossKernel)
    @assert size(k) == size(k′)
    if iszero(k)
        return k′
    elseif iszero(k′)
        return k
    else
        return CompositeCrossKernel(+, k, k′)
    end
end
function +(k::Kernel, k′::Kernel)
    @assert size(k) == size(k′)
    if iszero(k)
        return k′
    elseif iszero(k′)
        return k
    else
        return CompositeKernel(+, k, k′)
    end
end
function *(k::Kernel, k′::Kernel)
    @assert size(k) == size(k′)
    return iszero(k) || iszero(k′) ? zero(k) : CompositeKernel(*, k, k′)
end
function *(k::CrossKernel, k′::CrossKernel)
    @assert size(k) == size(k′)
    return iszero(k) || iszero(k′) ? zero(k) : CompositeCrossKernel(*, k, k′)
end


text:
module STMMRunner

using DocStringExtensions
using LinearAlgebra: norm
using Random: rand

include("config.jl")
include("utils.jl")
include("MSTM/MSTM.jl")

end


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

using ArraysOfArrays
using Test

using StaticArrays

@testset "functions" begin
    function gen_nested()
        A11 = [1 2; 3 4]
        A21 = [4 5 6; 7 8 9]
        A12 = [10 11; 12 13; 14 15]
        A22 = [16 17; 18 19]
        hcat([A11, A21], [A12, A22])
    end

    @testset "flatview and nestedview" begin
        A = [(@SArray randn(3, 2, 4)) for i in 1:2, j in 1:2]
        @test @inferred(nestedview(flatview(A), Val(3))) == A

        B = rand(3, 2, 4)
        @test @inferred(nestedview(flatview(B), SVector{3})) == @inferred(nestedview(B, Val(1)))
    end

    @testset "innersize" begin
        @test @inferred(innersize([[1, 2, 3], [4, 5, 6]])) == (3,)
        @test @inferred(innersize([[]])) == (0,)
        @test_throws DimensionMismatch @inferred(innersize([[1, 2, 3], [4, 5]]))
    end

    @testset "deepgetindex" begin
        A = gen_nested()
        @test @inferred(deepgetindex(A, 1, 2)) === A[1, 2]
        @test @inferred(deepgetindex(A, 1, 2, 2, 1)) == A[1, 2][2, 1]
        @test @inferred(deepgetindex(A, 1, 2, 1:2, 1)) == A[1, 2][1:2, 1]
        @test @inferred(deepgetindex(A, 1, 1:2, 2, 1)) == [A[1, 1][2, 1], A[1, 2][2, 1]]
        @test @inferred(deepgetindex(A, 1, 1:2, 1:2, 2)) == [A[1, 1][1:2, 2], A[1, 2][1:2, 2]]
        @test_throws MethodError deepgetindex(A, 1, 2, 2, 1, 2)

        B = rand(3, 4, 5)
        @test @inferred(deepgetindex(B, 2, 3, 4)) === getindex(B, 2, 3, 4)
        @test_throws MethodError deepgetindex(B, 6)
        @test_throws MethodError deepgetindex(B, 2, 3)
    end

    @testset "deepsetindex!" begin
        B12 = [20 21 22; 23 24 25]
        X = [[41, 42], [43, 44]]

        A = gen_nested()
        @test @inferred(deepsetindex!(A, B12, 1, 2)) === A
        @test A[1, 2] == B12

        A = gen_nested()
        @test @inferred(deepsetindex!(A, 42, 1, 2, 2, 1)) === A
        @test A[1, 2][2, 1] == 42

        A = gen_nested()
        @test @inferred(deepsetindex!(A, X, 1, 1:2, 1:2, 2)) === A
        @test deepgetindex(A, 1, 1:2, 1:2, 2) == X
    end

    @testset "deepview" begin
        A = gen_nested()
        @test @inferred(deepview(A, 1, 2)) == view(A, 1, 2)
        @test @inferred(deepview(A, 1, 2, 2, 1)) == view(A[1, 2], 2, 1)
        @test @inferred(deepview(A, 1, 2, 1:2, 1)) == view(A[1, 2], 1:2, 1)
        @test_throws ArgumentError deepview(A, 1, 1:2, 2, 1)
        @test @inferred(deepview(A, 1, 1:2, 1:2, 2)) == [A[1, 1][1:2, 2], A[1, 2][1:2, 2]]
        @test_throws MethodError deepview(A, 1, 2, 2, 1, 2)

        B = rand(3, 4, 5)
        @test @inferred(deepview(B, 2, 3, 4)) === view(B, 2, 3, 4)
        @test_throws MethodError deepview(B, 6)
        @test_throws MethodError deepview(B, 2, 3)
    end

    @testset "abstract_nestedarray_type" begin
        @test @inferred(abstract_nestedarray_type(Int, Val(()))) == Int
        @test @inferred(abstract_nestedarray_type(Int, Val((2,)))) == AbstractArray{Int, 2}
        @test @inferred(abstract_nestedarray_type(Float32, Val((2,3,4)))) ==
            AbstractArray{<:AbstractArray{<:AbstractArray{Float32, 4}, 3}, 2}
    end
end


text:
using DataStructures
lines = readlines("inputs/11-15/12.txt")

function edges(lines)
  e = DefaultDict(() -> String[])
  for l in lines
    a,b = split(l,'-')
    push!(e[a], b)
    push!(e[b], a)
  end
  e
end

small(x) = x != uppercase(x)

function cp(adj, start, goal, v=Set([]), twice_ = false)
  (start == goal) && return 1
  t = 0
  @show start, goal
  for ne in adj[start]
    twice = twice_
    if small(ne) && ne ∈ v
      if twice || ne ∈ ["start","end"]
        continue
      else
        twice = true
      end
    end
    push!(v, start)
    t += cp(adj, ne, goal, v, twice)
  end
  t
end

function day12()
  A = edges(lines)
  cp(A, "start", "end", Set([]), false)
end


text:
# Form request body
"""
    Form(dict::Dict)

A type representing a request body using the multipart/form-data encoding.
The key-value pairs in the Dict argument will constitute the name and value of each multipart boundary chunk.
Files and other large data arguments can be provided as values as IO arguments: either an `IOStream` such as returned via `open(file)`,
an `IOBuffer` for in-memory data. For complete control over a multipart chunk's details, an
`HTTP.Multipart` type is provided to support setting the `Content-Type`, `filename`, and `Content-Transfer-Encoding` if desired. See `?HTTP.Multipart` for more details.
"""
mutable struct Form <: IO
    data::Vector{IO}
    index::Int
    boundary::String
end

Form(f::Form) = f
Base.eof(f::Form) = f.index > length(f.data)
Base.isopen(f::Form) = false
Base.close(f::Form) = nothing
Base.length(f::Form) = sum(x->isa(x, IOStream) ? filesize(x) - position(x) : bytesavailable(x), f.data)
function Base.position(f::Form)
    index = f.index
    foreach(mark, f.data)
    return index
end
function Base.seek(f::Form, pos)
    f.index = pos
    foreach(reset, f.data)
    return
end

Base.readavailable(f::Form) = read(f)
function Base.read(f::Form)
    result = UInt8[]
    for io in f.data
        append!(result, read(io))
    end
    f.index = length(f.data) + 1
    return result
end

function Base.read(f::Form, n::Integer)
    nb = 0
    result = UInt8[]
    while nb < n
        d = read(f.data[f.index], n - nb)
        nb += length(d)
        append!(result, d)
        eof(f.data[f.index]) && (f.index += 1)
        f.index > length(f.data) && break
    end
    return result
end

function Form(d; boundary=string(rand(UInt128), base=16))
    # https://www.w3.org/Protocols/rfc1341/7_2_Multipart.html
    bcharsnospace = raw"\w'\(\)\+,-\./:=\?"
    boundary_re = Regex("^[$bcharsnospace ]{0,69}[$bcharsnospace]\$")
    @require match(boundary_re, boundary) !== nothing
    @require eltype(d) <: Pair
    data = IO[]
    io = IOBuffer()
    len = length(d)
    for (i, (k, v)) in enumerate(d)
        write(io, (i == 1 ? "" : "\r\n") * "--" * boundary * "\r\n")
        write(io, "Content-Disposition: form-data; name=\"$k\"")
        if isa(v, IO)
            writemultipartheader(io, v)
            seekstart(io)
            push!(data, io)
            push!(data, v)
            io = IOBuffer()
        else
            write(io, "\r\n\r\n")
            write(io, v)
        end
        i == len && write(io, "\r\n--" * boundary * "--" * "\r\n")
    end
    seekstart(io)
    push!(data, io)
    return Form(data, 1, boundary)
end

function writemultipartheader(io::IOBuffer, i::IOStream)
    write(io, "; filename=\"$(basename(i.name[7:end-1]))\"\r\n")
    write(io, "Content-Type: $(HTTP.sniff(i))\r\n\r\n")
    return
end
function writemultipartheader(io::IOBuffer, i::IO)
    write(io, "\r\n\r\n")
    return
end

"""
    Multipart(filename::String, data::IO, content_type=HTTP.sniff(data), content_transfer_encoding="")

A type to represent a single multipart upload chunk for a file. This type would be used as the value in a
key-value pair of a Dict passed to an http request, like `HTTP.post(url; body=Dict("key"=>HTTP.Multipart("MyFile.txt", open("MyFile.txt"))))`.
The `data` argument must be an `IO` type such as `IOStream`, or `IOBuffer`.
The `content_type` and `content_transfer_encoding` arguments allow the manual setting of these multipart headers. `Content-Type` will default to the result
of the `HTTP.sniff(data)` mimetype detection algorithm, whereas `Content-Transfer-Encoding` will be left out if not specified.

Filename SHOULD be included when the Multipart represents the contents of a file
[RFC7578 4.2](https://tools.ietf.org/html/rfc7578#section-4.2)

Content-Disposition set to "form-data" MUST be included with each Multipart.
An additional "name" parameter MUST be included
An optional "filename" parameter SHOULD be included if the contents of a file are sent
This will be formatted such as:
  Content-Disposition: form-data; name="user"; filename="myfile.txt"
[RFC7578 4.2](https://tools.ietf.org/html/rfc7578#section-4.2)

Content-Type for each Multipart is optional, but SHOULD be included if the contents
of a file are sent.
[RFC7578 4.4](https://tools.ietf.org/html/rfc7578#section-4.4)

Content-Transfer-Encoding for each Multipart is deprecated
[RFC7578 4.7](https://tools.ietf.org/html/rfc7578#section-4.7)

Other Content- header fields MUST be ignored
[RFC7578 4.8](https://tools.ietf.org/html/rfc7578#section-4.8)
"""
mutable struct Multipart{T <: IO} <: IO
    filename::Union{String, Nothing}
    data::T
    contenttype::String
    contenttransferencoding::String
    name::String
end

function Multipart(f::Union{AbstractString, Nothing}, data::T, ct::AbstractString="", cte::AbstractString="", name::AbstractString="") where {T<:IO}
    f = f !== nothing ? String(f) : nothing
    return Multipart{T}(f, data, String(ct), String(cte), String(name))
end

function Base.show(io::IO, m::Multipart{T}) where {T}
    items = ["data=::$T", "contenttype=\"$(m.contenttype)\"", "contenttransferencoding=\"$(m.contenttransferencoding)\")"]
    m.filename === nothing || pushfirst!(items, "filename=\"$(m.filename)\"")
    print(io, "HTTP.Multipart($(join(items, ", ")))")
end

Base.bytesavailable(m::Multipart{T}) where {T} = isa(m.data, IOStream) ? filesize(m.data) - position(m.data) : bytesavailable(m.data)
Base.eof(m::Multipart{T}) where {T} = eof(m.data)
Base.read(m::Multipart{T}, n::Integer) where {T} = read(m.data, n)
Base.read(m::Multipart{T}) where {T} = read(m.data)
Base.mark(m::Multipart{T}) where {T} = mark(m.data)
Base.reset(m::Multipart{T}) where {T} = reset(m.data)
Base.seekstart(m::Multipart{T}) where {T} = seekstart(m.data)

function writemultipartheader(io::IOBuffer, i::Multipart)
    if i.filename === nothing
        write(io, "\r\n")
    else
        write(io, "; filename=\"$(i.filename)\"\r\n")
    end
    contenttype = i.contenttype == "" ? HTTP.sniff(i.data) : i.contenttype
    write(io, "Content-Type: $(contenttype)\r\n")
    write(io, i.contenttransferencoding == "" ? "\r\n" : "Content-Transfer-Encoding: $(i.contenttransferencoding)\r\n\r\n")
    return
end

content_type(f::Form) = "Content-Type" =>
                        "multipart/form-data; boundary=$(f.boundary)"

# Deprecation can be removed in HTTP 0.10.0 (or 1.0.0).
@deprecate post(url, f::Form; kw...) post(url, [], f; kw...)


text:
@testset "Hebisch Problems" begin
    (x, ) = @variables x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Waldek*Hebisch - email*May*2013=#

    #= ::Subsection:: =#
    #=Problem #1=#

    @test_int [(x^6 - x^5 + x^4 - x^3 + 1)*exp(x), x, 25, 871*ℯ^x - 870*ℯ^x*x + 435*ℯ^x*x^2 - 145*ℯ^x*x^3 + 36*ℯ^x*x^4 - 7*ℯ^x*x^5 + ℯ^x*x^6]

    #= ::Subsection:: =#
    #=Problem #2=#

    @test_int [(2 - x^2)*exp(x/(x^2 + 2))/(x^3 + 2*x), x, -5, ExpIntegralEi(x/(2 + x^2))]

    @test_int [(2 + 2*x + 3*x^2 - x^3 + 2*x^4)*exp(x/(2 + x^2))/(x^3 + 2*x), x, -5, ℯ^(x/(2 + x^2))*(2 + x^2) + ExpIntegralEi(x/(2 + x^2))]

    #= ::Subsection:: =#
    #=Problem #3=#

    @test_int [(exp(x) + 1)*(exp(exp(x) + x)/(exp(x) + x)), x, 2, ExpIntegralEi(ℯ^x + x)]

    #= ::Subsection:: =#
    #=Problem #4=#

    @test_int [(x^3 - x^2 - 3*x + 1)*(exp(1/(x^2 - 1))/(x^3 - x^2 - x + 1)), x, -6, ℯ^(1/(-1 + x^2))*(1 + x)]

    #= ::Subsection:: =#
    #=Problem #5=#

    @test_int [(log(x)^2 - 1)*exp(1 + 1/log(x))/log(x)^2, x, 1, x*ℯ^(1 + 1/log(x))]

    @test_int [((x + 1)*log(x)^2  - 1)*exp(x + 1/log(x))/log(x)^2, x, -2, ℯ^(x + 1/log(x))*x]
    end

max_stars_repo_path: test/0 Independent test suites/Hebisch Problems.jl


text:
struct BasicReg{L,R} <: RegressionModel
    nobs::Int
    formula::FormulaTerm{L,R}
    coef::Vector{Float64}
    coefnames::Vector{String}
    yname::String
    tss::Float64
    rss::Float64
    residuals::Union{Vector{Float64}, Nothing}
    function BasicReg(nobs, formula::FormulaTerm{L,R}, coef, coefnames, yname, tss, rss, residuals) where {L,R}
        # @assert rss >= 0 "Residual sum of squares must be greater than 0"
        # @assert tss >= 0 "Total sum of squares must be greater than 0"
        # @assert nobs >= 0 "Observations must be greater than 0"
        # @assert length(coef) == length(coefnames) "Number of coefficients must be same as number of coefficient names"
        new{L,R}(nobs, formula, coef, coefnames, yname, tss, rss, residuals)
    end
    BasicReg(x::Int, f::FormulaTerm{L,R}) where {L, R} = new{L,R}(x, f)
end

# These are created to minimize the amount of allocations Julia does
# Julia typically allocates a vector for each loop, which when using
# so many loops, can create real garbage collection problems
# As it turns out, doing sum(abs2, resp - mean(resp)) also does
# an allocation, which could mean allocating a huge amount
# caluclating rss was even worse, so these functions are only
# meant to be used internally but do not allocate if passed a view
function calc_tss(resp)
    out = 0.0
    m = mean(resp)
    @simd for x in resp
        out += (x - m) ^ 2
    end
    out
end
function fast_pred(pred, coef, i)
    out = 0.0
    @simd for j in 1:length(coef)
        @inbounds out += pred[i, j] * coef[j]
    end
    out
end
function calc_rss(resp, pred, coef)
    out = 0.0
    @simd for i in 1:length(resp)
        @inbounds out += (resp[i] - fast_pred(pred, coef, i)) ^ 2
    end
    out
end
"""
    function BasicReg(
        resp::AbstractVector{Float64},
        pred::AbstractMatrix{Float64},
        yname::String,
        xnames::Vector{String},
        f::FormulaTerm{L,R};
        save_residuals::Bool=false,
        minobs::Int=1
    )::BasicReg{L,R} where {L,R}

## Arguments
- resp::AbstractVector{Float64}: The "Y" or response in a linear regression
- pred::AbstractMatrix{Float64}: The "X" matrix in a linear regression
- yname::String: The name of the response variable
- xnames::Vector{String}: The names of the prediction variables
- f::FormulaTerm{L,R}: A StatsModels.jl formula, saved in the resulting struct
- save_residuals::Bool=false: Whether or not to save the vector of residuals from
    the regression. Note for large numbers of regressions this can significantly slow
    down the speed
- minobs::Int=1: The minimum length of the response vector for the regression to
    run. The regression will also not run if the length of the response vector is
    less than or equal to the number of columns in the prediction matrix.

BasicReg is an intentionally simplistic linear regression. It also attempts to produce
a minimum number of allocations if views of vectors are passed.
"""
function BasicReg(
    resp::AbstractVector{Float64},
    pred::AbstractMatrix{Float64},
    yname::String,
    xnames::Vector{String},
    f::FormulaTerm{L,R};
    save_residuals::Bool=false,
    minobs=1
)::BasicReg{L,R} where {L,R}
    if length(resp) <= size(pred, 2) || length(resp) < minobs
        return BasicReg(length(resp), f)
    end

    coef = cholesky!(Symmetric(pred' * pred)) \ (pred' * resp)

    BasicReg(
        length(resp),
        f,
        coef,
        xnames,
        yname,
        calc_tss(resp),
        calc_rss(resp, pred, coef),
        save_residuals ? resp - pred * coef : nothing
    )
end

function create_pred_matrix(data::IterateTimelineTable, sch)
    id = first(data)[1]
    create_pred_matrix(
        parent(data)[id, internal_termvars(sch.rhs), AllowMissing{true}],
        sch
    )
end

function create_pred_matrix(data::TimelineTable{true}, sch)
    # This allows missing to make sure all dates are in the data

    mat = modelcols(sch.rhs, data)
    DataMatrix(
        mat,
        data_dates(data),
        calendar(data)
    )
end

"""
    quick_reg(
        data::TimelineTable{false},
        f::FormulaTerm;
        minobs::Real=0.8,
        save_residuals::Bool=false
    )

    quick_reg(
        data::IterateTimelineTable,
        f::FormulaTerm;
        minobs::Real=0.8,
        save_residuals::Bool=false
    )

Calculates a linear regression for the supplied data based on the formula (formula from StatsModels.jl).
Unless the formula explicitly excludes the intercept (i.e., `@formula(y ~ 0 + x)`), an intercept is added.

If `data` is of the type `IterateTimelineTable`, then the formula is applied to each `TimelineTable` in an
optimized way and returns a `Vector{BasicReg}`.

## Arguments
- `minobs::Real`: The minimum number of observations to return a completed regression. If less than 1,
    the value is used as a percentage relative to the total number of business days in the time period.
    Therefore, the default of 0.8 corresponds to at least 80% of the business days over the time period have values.
- `save_residuals::Bool=false`: Whether to save the residuals into `BasicReg`, This can have significant performance implications.

"""
function quick_reg(
    data::TimelineTable{false},# only works if no missing data for multiplication
    f::FormulaTerm;
    minobs::V=0.8,
    save_residuals::Bool=false
) where {V<:Real}

    if !StatsModels.omitsintercept(f) & !StatsModels.hasintercept(f)
        f = FormulaTerm(f.lhs, InterceptTerm{true}() + f.rhs)
    end
    sch = apply_schema(f, schema(f, data))
    select!(data, internal_termvars(sch))

    resp, pred = modelcols(sch, data)
    yname, xnames = coefnames(sch)
    BasicReg(
        resp,
        pred,
        yname,
        xnames,
        f;
        save_residuals,
        minobs=adjust_minobs(minobs, calendar(data), norm_dates(data))
    )
end

function fill_vector_reg!(
    out_vector::Vector{BasicReg{L,R}},
    iter_data::IterateTimelineTable{T},
    cache::DataMatrix,
    sch::FormulaTerm,
    f::FormulaTerm{L,R};
    save_residuals::Bool=false,
    minobs=0.8
) where {T, L, R}
    @assert validate_iterator(iter_data, out_vector) "Length of out_vector does not match the number of indexes in the iterator"
    cols = internal_termvars(sch)
    yname, xnames = coefnames(sch)

    Threads.@threads for (u_id, iter_indexes) in iter_data
        data = parent(iter_data)[u_id, cols, AllowMissing{true}]
        resp = datavector_modelcols(int_lhs(sch), data)
        for iter in iter_indexes
            cur_dates = dates_min_max(data_dates(data), iter_dates(iter))
            if nnz(data_missing_bdays(data)) == 0
                x = view(resp, cur_dates)
                y = view(cache, cur_dates)
            else
                new_mssngs = get_missing_bdays(data, cur_dates)
                x = view(resp, cur_dates, new_mssngs)
                y = view(cache, cur_dates, new_mssngs)
            end
            @inbounds out_vector[iter_index(iter)] = BasicReg(
                x,
                y,
                yname,
                xnames,
                f;
                save_residuals,
                minobs=adjust_minobs(minobs, calendar(parent(data)), iter_dates(iter))
            )
        end
    end
    out_vector
end

function quick_reg(
    data::IterateTimelineTable,
    f::FormulaTerm;
    minobs::V=0.8,
    save_residuals::Bool=false
) where {V<:Real}
    if !StatsModels.omitsintercept(f) & !StatsModels.hasintercept(f)
        f = FormulaTerm(f.lhs, InterceptTerm{true}() + f.rhs)
    end
    sch = apply_schema(f, schema(f, parent(data)))
    cache = create_pred_matrix(data, sch)
    out = fill(BasicReg(0, f), total_length(data))
    fill_vector_reg!(
        out,
        data,
        cache,
        sch,
        f;
        save_residuals,
        minobs
    )
end

StatsBase.predict(mod::BasicReg, x) = x * coef(mod)

StatsBase.coef(x::BasicReg) = isdefined(x, :coef) ? x.coef : missing
StatsBase.coefnames(x::BasicReg) = isdefined(x, :coef) ? x.coefnames : missing
StatsBase.responsename(x::BasicReg) = isdefined(x, :coef) ? x.yname : missing
StatsBase.nobs(x::BasicReg) = x.nobs
StatsBase.dof_residual(x::BasicReg) = isdefined(x, :coef) ? nobs(x) - length(coef(x)) : missing
StatsBase.r2(x::BasicReg) = 1 - (rss(x) / deviance(x))
StatsBase.adjr2(x::BasicReg) = 1 - rss(x) / deviance(x) * (nobs(x) - 1) / dof_residual(x)
StatsBase.islinear(x::BasicReg) = true
StatsBase.deviance(x::BasicReg) = isdefined(x, :coef) ? x.tss : missing
StatsBase.rss(x::BasicReg) = isdefined(x, :coef) ? x.rss : missing
function StatsBase.residuals(x::BasicReg)
    if !isdefined(x, :coef)
        return missing
    end
    if x.residuals === nothing
        @error("To access residuals, run `cache_reg` with the option `save_residuals=true`")
    else
        x.residuals
    end
end

function rhs_str(nms, vals; intercept = "(Intercept)")
    out = String[]
    for (nm, val) in zip(nms, vals)
        val_str = string(round(val, digits=3))
        if nm == intercept
            push!(out, val_str)
        else
            push!(out, val_str * "*" * nm)
        end
    end
    join(out, " + ")
end

function Base.show(io::IO, rr::BasicReg)
    if !isdefined(rr, :coef)
        print(io, "Obs: $(nobs(rr)), $(rr.formula)")
    else
        print(io, "Obs: $(nobs(rr)), $(responsename(rr)) ~ $(rhs_str(coefnames(rr), coef(rr))), AdjR2: ", round(adjr2(rr) * 100, digits=3), "%")
        # mat = hcat(
        #     coefnames(rr),
        #     string.(round.(coef(rr), digits=3))
        # )
        # PrettyTables.pretty_table(io, mat)
    end
end


text:
#!/usr/bin/env julia

#include("../src/MsgPackRpcClient.jl")
#using MsgPackRpcClient

tests = ["sync_call", "async_call", "arguments"]

println("Running tests:")

#@sync begin

addprocs(1)
println("INFO: Start test_server.rb")
@spawn remotecall(2, run(`ruby test_server.rb`))

for t in tests
  test_fn = "$t.jl"
  println(" * $test_fn")
  run(`julia $test_fn`)
end

println("INFO: Stop test_server.rb")
try
  run(`pkill ruby`)
catch
end
println("INFO: Finished.")

#end

max_stars_repo_name: keithseahus/msgpack-rpc-julia

max_issues_repo_name: keithseahus/msgpack-rpc-julia

max_forks_repo_name: keithseahus/msgpack-rpc-julia