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

using Markdown
using InteractiveUtils

# ╔═╡ 723b4419-4413-431a-9a2b-5a501cbd2607
begin
	import Pkg;
    Pkg.activate(joinpath(homedir(),"Code","RvSpectMLEcoSystem","NeidSolarScripts"))
	begin
	using CSV, JLD2, FileIO
	using DataFrames, Query, DataValues
	using Dates
	using NeidSolarScripts
	using Polynomials
	using Statistics, NaNMath
	using PlutoUI
	using Plots, ColorSchemes
end
end

# ╔═╡ 67b3208b-428b-4f48-94e5-3cb11beee67d
begin
	path = "/mnt/data_simons/NEID/pyrohelio"
	files = readdir(path)
end

# ╔═╡ b2c1942e-7d8a-4b5e-afba-d42aff7f6471
begin
	df = DataFrame()
	for f in files
		if !occursin("neid_ljpyrohelio_chv0",f) continue end
		df_tmp = CSV.read(joinpath(path,f), header=[:time,:val1,:val2], DataFrame) |> @map({time=DateTime(_.time[1:23]),flux=_.val2}) |> DataFrame
		append!(df,df_tmp)
	end
	size(df)
end

# ╔═╡ b677c0ee-0c58-4a48-8d19-d1c54f7771bb
days_to_keep = 360

# ╔═╡ 90172c7e-3339-49da-bd3c-27a7fb0dd3f0
df_by_date = df |> @take(days_to_keep*24*60*60) |> @mutate( date=Date(_.time #=-Dates.Hour(13)=#) ) |> @groupby(_.date) |> @map({ date=_.date, time=_.time, flux=_.flux }) |> DataFrame;

# ╔═╡ 003ff73f-aad3-40bf-a673-7c28091e51e4
begin
	datestr = "20210312"
	manifest = CSV.read(joinpath(path,datestr,"manifest.csv"), DataFrame)
end

# ╔═╡ 8b662f93-b270-4ce4-ace0-ac3f3bb1dc20
first(df_by_date[1,:date])

# ╔═╡ 1276658b-97fc-4859-84e2-8df7cd43d1f6
m = match(r"neidL1_(\d{8}T\d{6})\.fits",manifest.Filename[3])

# ╔═╡ 1cffcbe4-fa66-4cbb-835b-4ee8e29b87b5
idx_date = findfirst(isequal(Date(m.captures[1][1:8],"YYYYmmdd")),map(i->first(df_by_date[i,:date]),1:size(df_by_date,1)))

# ╔═╡ e4dc87e4-f7d4-462a-a24f-d22ee036ba9b
files

# ╔═╡ 24ff107b-f337-4667-a09d-4e0342b085de
pwd()

# ╔═╡ 34f82d2f-4740-4422-9a17-e35702b51556
manifest[!,:solar_hour_angle] = SolarRotation.get_solar_hour_angle(manifest[!,:bjd], obs=:WIYN)
#manifest[!,:solar_hour_angle] = map(i->SolarRotation.get_solar_hour_angle(manifest[i,:bjd], obs=:WIYN),1:size(manifest,1))

# ╔═╡ 451d154c-29ad-4862-894a-d96dc97e8f12
names(manifest)

# ╔═╡ 7fde57bb-8bb9-499b-89f5-3f31fc5050ce
df_out = select(manifest, :Filename => (x->basename.(x))=>:filename, :bjd, :pyrheliometer_flux_over_model, :pyrheliometer_rms_flux, :alt_sun, :solar_hour_angle, :Δv_diff_ext, :Δfwhm²)

# ╔═╡ a013b123-05c8-4054-bc69-eedab7da28cc
CSV.write(datestr * "_solar_info.csv", df_out)

# ╔═╡ cb47c02c-cd8f-4a05-a3f4-d7baaf6afb8d
size(manifest)

# ╔═╡ b5afbe2c-b3c6-4825-9804-9df9d22c9a56
df_by_date[idx_date,:time]|>collect

# ╔═╡ a83c43ba-fc15-4756-bebb-53d29635db5c
 SolarRotation.get_solar_info(manifest.bjd[1], obs=:WIYN)

# ╔═╡ f5bfccb9-109b-48ae-ab4c-39b72ff070c8
SolarRotation.get_solar_hour_angle(manifest.bjd[1], obs=:WIYN)

# ╔═╡ fe8a34e8-afc6-11eb-1ff5-dfac1585fd04
#=
begin
	local plt3 = plot(legend=:none)
	local pal = palette(:berlin, length(df_plt))
	for day_idx in idx_clear_days[rms_list.<=0.015]
		if size(df_plt[day_idx],1) > 1
		local t = Time.(df_plt[day_idx][:,:times]).-Hour(7)
		local model = flux_model.(df_plt[day_idx][!,:airmass],df_plt[day_idx][!,:sol_dist])
		local mask = (0.95 .<= df_plt[day_idx][!,:flux]./model .<= 1.05)
		if sum(mask) >= 4
			local coeff_daily = fit(Dates.value.(t[mask].-t[1]),df_plt[day_idx][mask,:flux]./model[mask],2)

			model .*= coeff_daily.(Dates.value.(t.-t[1]))
			local norm = quantile(df_plt[day_idx][mask,:flux]./model[mask],0.75)
			model .*= norm
			local mask2 = (0.995 .<= df_plt[day_idx][:,:flux]./model .<= 1.005)
			#plot!(plt3, df_plt[day_idx][!,:airmass], model, color=pal[day_idx])
			scatter!(plt3,df_plt[day_idx][!,:airmass],df_plt[day_idx][mask,:flux]./(model[mask]) ,ms=1.0, label=:none)
			scatter!(plt3,df_plt[day_idx][!,:airmass],df_plt[day_idx][mask2,:flux]./(model[mask2]) ,ms=2.5, label=:none)
		end
		end
		ylims!(plt3,0.98,1.015)
	end
	plt3
end
=#

# ╔═╡ d56c42dc-66b9-4730-99a6-fdb18debf798
function get_idx_in_pyro_data(d_start::DateTime, tvec::AVD; exp_time::Real = 55) where {AVD <: AbstractVector{DateTime} }
	idx_lo = searchsortedfirst(tvec,d_start)
	d_stop = julian2datetime(datetime2julian(d_start)+exp_time/(24*60*60))
	idx_hi = searchsortedlast(tvec,d_stop)
	if idx_lo >1 idx_lo -= 1 end
	if idx_hi <length(tvec) idx_hi += 1 end
	return (idx = idx_lo:idx_hi, time_start = tvec[idx_lo], time_stop=tvec[idx_hi])
end

# ╔═╡ fff4ec82-2dd3-4d56-9faf-4b309bfee108
get_idx_in_pyro_data(DateTime(m.captures[1],"YYYYmmddTHHMMSS"), df_by_date[idx_date,:time] )

# ╔═╡ d9f7d500-7f43-4217-9946-7c20647d3940
idx_sec = get_idx_in_pyro_data(julian2datetime(manifest[3,:bjd]), (df_by_date[idx_date,:time]) )

# ╔═╡ 3c9725b1-ea31-40ab-b61c-90e1f07c23ea
function get_alt_airmass(d::DateTime)
	t = datetime2julian(d)
	return get_alt_airmass(t)
end

# ╔═╡ 72c08181-2a6f-4475-abfa-e89ca5a4965e
function get_alt_airmass(t::Real)
	result = SolarRotation.get_solar_info(t, obs=:WIYN)
	return (alt = result[:alt], airmass = result[:airmass], sol_dist = result[:sol_dist_au])
end

# ╔═╡ 51c3761e-0a85-4b37-90b0-499538120f4c
df_obs = DataFrame(get_alt_airmass.(df_by_date[idx_date,:time][idx_sec]))

# ╔═╡ 86c3d69a-57d2-4828-bb75-f9b88f6036a5
function flux_model(airmass::Real, dist::Real = 1)
	poly = Polynomial([7.061980163705281, -0.07628424166629734])
	exp(poly(airmass))/ dist^2
end

# ╔═╡ cf5b8699-6957-4067-a4a1-f32780ef2644
function flux_model_obs(dt::DateTime; exp_time::Real = 55)
	d = Date(dt)
	idx_date = findfirst(isequal(d),map(i->first(df_by_date[i,:date]),1:size(df_by_date,1)))
	(idx_sec, tlo, thi) = get_idx_in_pyro_data(dt, dropna(df_by_date[idx_date,:time]), exp_time=exp_time )
	#return (1.0, length(idx_sec))
	#(alts, airmasses, soldists) = get_alt_airmass.(df_by_date[idx_date,:time][idx_sec])
	#fm = flux_model.(df_obs.airmass,df_obs.sol_dist)
	df_obs_start = get_alt_airmass(df_by_date[idx_date,:time][first(idx_sec)])
	df_obs_stop  = get_alt_airmass(df_by_date[idx_date,:time][last(idx_sec)])
	flux_start = flux_model(df_obs_start[:airmass],df_obs_start[:sol_dist])
	flux_stop  = flux_model(df_obs_stop[:airmass],df_obs_stop[:sol_dist])
	flux_mean = 0.5*(flux_start+flux_stop)
	#flux_slope = (flux_stop-flux_start)/Dates.value(Dates.Second(thi-tlo))
	#return df_by_date[idx_date,:flux][idx_sec] |> collect
	#return dropna(df_by_date[idx_date,:flux])[idx_sec], std(df_by_date[idx_date,:flux][idx_sec] |> collect)
	flux_ratio = dropna(df_by_date[idx_date,:flux])
	flux_ratio = flux_ratio[idx_sec]./flux_mean
	mean_flux_ratio = mean(flux_ratio)
	#flux_ratio = dropna(df_by_date[idx_date,:flux])[idx_sec] ./flux_mean
	#flux_ratio = mean(dropna(df_by_date[idx_date,:flux][idx_sec])) ./flux_mean
	#return flux_ratio #, std(flux_ratio) )
	#NaNMath.std(flux_ratio)
	return (ave = mean_flux_ratio, rms = std(flux_ratio) )
	#slope = flux_slope, idx_date=idx_date, idx_sec=idx_sec, flux_ratio = flux_ratio)
end

# ╔═╡ 5af436fe-bb4a-4069-bb82-42259d7d57c0
mean(flux_model.(df_obs.airmass,df_obs.sol_dist))

# ╔═╡ 23bd0a0b-2ef3-4acd-8ab6-d61f8cdd9c2b
df_by_date[188

# ╔═╡ 63de83e7-91c8-41aa-91de-51f280534c67
function flux_model_obs(fn::String)
	m = match(r"neidL1_(\d+T\d+)\.fits",basename(fn))
	dt = DateTime(m.captures[1],"yyyymmddTHHMMSS")
	flux_model_obs(dt)
end

# ╔═╡ 8f484f4a-c194-474e-825c-4301aae07e92
function normalized_flux_obs(df::DataFrame)
	flux_ratios = zeros(size(df,1))
	std_flux_ratios = zeros(size(df,1))
	#idx_rngs = fill(0:1,size(df,1))
	for (i,row) in enumerate(eachrow(df))
		#flux_ratio = flux_model_obs.(df[!,:Filename])
		(mean_flux_ratio, std_flux_ratio) = flux_model_obs(row.Filename)
		flux_ratios[i] = mean_flux_ratio
		std_flux_ratios[i] = std_flux_ratio
	end

		local mask = (0.95 .<= flux_ratios .<= 1.05)
		if sum(mask) >= 4
			coeff_daily = fit((df[mask,:bjd].-manifest[1,:bjd]),dropna(flux_ratios[mask]),2)
			flux_ratios ./= coeff_daily.(df[:,:bjd].-df[1,:bjd])
			flux_ratios ./= quantile(dropna(flux_ratios[mask]),0.75)
		end
	return (flux_ratios, std_flux_ratios)
end

# ╔═╡ b671c156-4ea7-4d4b-93fa-69f632e3da18
for (i,f) in enumerate(files)
	m = match(r"(\d{8})$", f)
	if m == nothing continue end
	datestr = m.captures[1]
	manifest_fn = joinpath(path,datestr,"manifest.csv")
	println("# Processing ",manifest_fn)
	if !isfile(manifest_fn)
		println("# Can't find ", manifest_fn)
		continue
	end
	output_fn = joinpath(path,datestr * "_solar_info.csv")
	if filesize(output_fn) > 0 continue end
	manifest = CSV.read(manifest_fn, DataFrame)
	println("# Read ", size(manifest,1), " rows.")
	if size(manifest,1) == 0 continue end
	(fm, fs) = normalized_flux_obs(manifest)
	manifest[!,:pyrheliometer_flux_over_model] = fm
	manifest[!,:pyrheliometer_rms_flux] = fs
	manifest[!,:solar_hour_angle] = SolarRotation.get_solar_hour_angle(manifest[!,:bjd], obs=:WIYN)
	df_out = select(manifest, :Filename => (x->basename.(x))=>:filename, :bjd, :pyrheliometer_flux_over_model, :pyrheliometer_rms_flux, :alt_sun, :solar_hour_angle, :Δv_diff_ext, :Δfwhm²)

	CSV.write(output_fn, df_out)
	#if i>3 break end
end

# ╔═╡ 89beb337-444f-4953-a6ed-a51c75d22cd9
(fm, fs) = normalized_flux_obs(manifest)

# ╔═╡ 0a6793ad-99e4-4d2a-9ec6-a82f2e1278da
scatter(manifest[!,:bjd],fm)

# ╔═╡ e34fbb33-05c0-48bf-b4b8-3b5d74d77d57
scatter(fm,fs,xlims=(0.9,1.02),ylims=(0.,0.1))

# ╔═╡ 4471d7dd-f26d-4a0d-9f37-de86d736579d
manifest[!,:pyrheliometer_flux_over_model] = fm

# ╔═╡ b96e3c2e-87e4-4311-9887-f5692d5e2259
manifest[!,:pyrheliometer_rms_flux] = fs

# ╔═╡ b4ba194d-8123-4ddf-b80b-f79470926d6a
size(fs)

# ╔═╡ 657c2fbc-5d68-4549-93ac-aabee96a4304
#DateTime(Date(now()),Time(Dates.hour(now()),Dates.minute(now()),Dates.second(now())+55))
#(Date(now()) +
	Dates.Second(Time(0,30,46)-Time(0,23,46))
#Dates.minute(now())

# ╔═╡ c9e6b8fe-df4d-4c56-85f2-1e8d9b38a70b
df_by_date[idx_date,:flux][idx_sec]./1000

# ╔═╡ Cell order:
# ╠═723b4419-4413-431a-9a2b-5a501cbd2607
# ╠═67b3208b-428b-4f48-94e5-3cb11beee67d
# ╠═b2c1942e-7d8a-4b5e-afba-d42aff7f6471
# ╠═b677c0ee-0c58-4a48-8d19-d1c54f7771bb
# ╠═90172c7e-3339-49da-bd3c-27a7fb0dd3f0
# ╠═003ff73f-aad3-40bf-a673-7c28091e51e4
# ╠═8b662f93-b270-4ce4-ace0-ac3f3bb1dc20
# ╠═1cffcbe4-fa66-4cbb-835b-4ee8e29b87b5
# ╠═1276658b-97fc-4859-84e2-8df7cd43d1f6
# ╠═fff4ec82-2dd3-4d56-9faf-4b309bfee108
# ╠═d9f7d500-7f43-4217-9946-7c20647d3940
# ╠═0a6793ad-99e4-4d2a-9ec6-a82f2e1278da
# ╠═e34fbb33-05c0-48bf-b4b8-3b5d74d77d57
# ╠═b671c156-4ea7-4d4b-93fa-69f632e3da18
# ╠═e4dc87e4-f7d4-462a-a24f-d22ee036ba9b
# ╠═24ff107b-f337-4667-a09d-4e0342b085de
# ╠═89beb337-444f-4953-a6ed-a51c75d22cd9
# ╠═4471d7dd-f26d-4a0d-9f37-de86d736579d
# ╠═b96e3c2e-87e4-4311-9887-f5692d5e2259
# ╠═34f82d2f-4740-4422-9a17-e35702b51556
# ╠═451d154c-29ad-4862-894a-d96dc97e8f12
# ╠═7fde57bb-8bb9-499b-89f5-3f31fc5050ce
# ╠═a013b123-05c8-4054-bc69-eedab7da28cc
# ╠═b4ba194d-8123-4ddf-b80b-f79470926d6a
# ╠═cb47c02c-cd8f-4a05-a3f4-d7baaf6afb8d
# ╠═8f484f4a-c194-474e-825c-4301aae07e92
# ╠═cf5b8699-6957-4067-a4a1-f32780ef2644
# ╠═51c3761e-0a85-4b37-90b0-499538120f4c
# ╠═5af436fe-bb4a-4069-bb82-42259d7d57c0
# ╠═b5afbe2c-b3c6-4825-9804-9df9d22c9a56
# ╠═a83c43ba-fc15-4756-bebb-53d29635db5c
# ╠═f5bfccb9-109b-48ae-ab4c-39b72ff070c8
# ╠═fe8a34e8-afc6-11eb-1ff5-dfac1585fd04
# ╠═d56c42dc-66b9-4730-99a6-fdb18debf798
# ╠═3c9725b1-ea31-40ab-b61c-90e1f07c23ea
# ╠═72c08181-2a6f-4475-abfa-e89ca5a4965e
# ╠═86c3d69a-57d2-4828-bb75-f9b88f6036a5
# ╠═23bd0a0b-2ef3-4acd-8ab6-d61f8cdd9c2b
# ╠═63de83e7-91c8-41aa-91de-51f280534c67
# ╠═657c2fbc-5d68-4549-93ac-aabee96a4304
# ╠═c9e6b8fe-df4d-4c56-85f2-1e8d9b38a70b


text:
showprogram(core::AsapCore) = show(STDOUT, core.program, core.pc)

# Method for quickly summarizing a core.
function summarize(core::AsapCore)
    # Get the stall signals since these are evaluated on each clock cycle and
    # not stored in core state.

    # Display various processor state.
    show_hwstate(core)
    println()
    show_condexec(core)
    println()
    show_io(core)
    println()
    show_pipeline(core)
end

################################################################################
#                           SHOW ALU FLAGS
################################################################################
function Base.show(io::IO, f::ALUFlags)
    print(io, "ALU Flags: ")
    print(io, "C = $(f.carry), ")
    print(io, "N = $(f.negative), ")
    print(io, "O = $(f.overflow), ")
    print(io, "Z = $(f.zero) ")
end

################################################################################
#                           SHOW PIPELINE ENTRY
################################################################################
function Base.show(io::IO, p::PipelineEntry)
    # Show the instruction.
    show(io, p.instruction)

    # Show the values
    print(io, "src1_v=$(p.src1_value) ")
    print(io, "src2_v=$(p.src2_value) ")
    print(io, "result=$(p.result) ")
end

################################################################################
#                          SHOW ADDRESS GENERATOR
################################################################################
function Base.show(io::IO, ag::AddressGenerator)
    # start:step:stop -- current
    print(io, "$(ag.start):$(ag.stride):$(ag.stop) -- $(ag.current)")
end

################################################################################
#                        SHOW CONDITIONAL EXECUTION
################################################################################
function Base.show(io::IO, c::CondExec)
    print(io, "Flag = $(c.flag), ")
    # Mask has 14 meaningful bits. Show the bits directly.
    print(io, "Mask = $(bits(c.mask)[end - 14:end]), ")
    # Print Reduction
    print(io, "Reduction = $(c.unary_op), ")

    # Interpret the bits for faster decoding.
    mask = c.mask

    # Only print out the summary if there's something to show.
    print(io, "Mask Summary: ")
    isbitset(mask, 13) && print(io, "EK ")
    # Bits 11 and 12 encode the unary op, which is already diaplayed
    isbitset(mask, 10) && print(io, "MO ")
    isbitset(mask,  9) && print(io, "PO ")
    isbitset(mask,  8) && print(io, "OBUF ")
    isbitset(mask,  7) && print(io, "MI ")
    isbitset(mask,  6) && print(io, "PI ")
    isbitset(mask,  5) && print(io, "IBUF1 ")
    isbitset(mask,  4) && print(io, "IBUF0 ")
    isbitset(mask,  3) && print(io, "O ")
    isbitset(mask,  2) && print(io, "C ")
    isbitset(mask,  1) && print(io, "Z ")
    isbitset(mask,  0) && print(io, "N ")

end

################################################################################
#                          SHOW PROCESSOR STATE
################################################################################

function show_hwstate(core::AsapCore)
    printstyled("HW State\n"; bold = true, color = :white)
    println("PC: $(core.pc)")

    # Create an auto-align object for printing information horizontally.
    aa = AutoAligns.AutoAlign(align = Dict(:default => AutoAligns.left))

    # Line 1
    print(aa, "Pending Nops: ", core.pending_nops, "   ")
    print(aa, "HW Pointer 0: ", core.pointers[1], "   ")
    print(aa, "AG 0: ", core.address_generators[1], "   ")
    println(aa)

    # Line 2
    print(aa, "Return Address: ", core.return_address, "   ")
    print(aa, "HW Pointer 1: ",   core.pointers[2], "   ")
    print(aa, "AG 1: ", core.address_generators[2], "   ")
    println(aa)

    # Line 3
    print(aa, "Repeat Start: ", core.repeat_start, "   ")
    print(aa, "HW Pointer 2: ", core.pointers[3], "   ")
    print(aa, "AG 2: ", core.address_generators[3], "   ")
    println(aa)

    # Line 4
    print(aa, "Repeat End: ", core.repeat_end, "   ")
    print(aa, "HW Pointer 3: ", core.pointers[4], "   ")
    println(aa)

    # Line 5
    print(aa, "Repeat Count: ", core.repeat_count, "   ")
    println(aa)

    # Print the AutoAlign object
    print(STDOUT, aa)

    # --- ALU Flags and Accumulator --- #
    println()
    println(core.aluflags)
    println()
    println("Accumulator: $(core.accumulator)")
    println()

    # Print out branch mispredict
    if checkbranch(core)
        printstyled("Branch Mispredicted\n", color = :red)
    else
        printstyled("No Branch Mispredict\n", color = :green)
    end

    # Print out stall reason if the core is stalled.
    stall_reason = stall_fifo_check(core)
    if stall_reason != NoStall
        println()
        println("Stall reason: $stall_reason")
    end
end

function show_condexec(core::AsapCore)
    printstyled("Conditional Execution\n", bold = true, color = :white)
    for (index, condexec) in enumerate(core.condexec)
        println("CX$(index-1): $condexec")
    end
end

# Show the status of the input and output fifos.
function show_io(core::AsapCore)
    printstyled("Input Fifo Occupancy\n"; bold = true, color = :white)
    for (index, fifo) in enumerate(core.fifos)
        color = isreadready(fifo) ? :green : :red
        printstyled("Fifo $(index -1): $(read_occupancy(fifo))\n", color = color)
    end

    println()
    printstyled("Output Fifo Occupancy\n"; bold = true, color = :white)
    for (k,v) in core.outputs
        println("Output $k: $(write_occupancy(v))")
    end
    # Declare boldly if no output fifos have been connected.
    if length(core.outputs) == 0
        printstyled("No Connected Outputs\n", color = :red)
    end

    # --- Show the output mask --- #
    println()
    printstyled("OBUF Mask\n", bold = true, color = :white)

    # Convert the obuf mask to a string
    obuf_mask_string = join([x ? "1" : "0" for x in core.obuf_mask])
    println(obuf_mask_string)

    # Print out selected indices for easier viewing.
    for (index, val) in enumerate(core.obuf_mask)
        val && print(io, index-1, " ")
    end
end

# Show the pipeline to STDOUT. Take stall signals as an argument to color
# stages according to if they are stalled or not.
#
# Color has the following meaning:
#
# WHITE: Stage is executing normally
# RED:   Stage is stalled
# CYAN:  Stage is inserting NOPS.
function show_pipeline(core::AsapCore)
    stall_signals = stall_check(core)
    mispredict = checkbranch(core)

    printstyled("Pipeline\n"; bold = true, color = :white)

    # Show stages 0 and 1. For Stage 1, create a pipe stage for the current
    # value of the PC so we can se what's next.
    color = stall_signals.stall_01 ? (:red) : (:white)

    # Print stage 0
    printstyled("Stage 0: No Instruction\n", color = color)

    # Run Stage 0 to get the next program counter.
    s0_nextstate = pipeline_stage0(core, stall_signals.stall_01, mispredict)

    # Create a instruction for stage 1. Make NOP if PC is out of bounds.
    printstyled("Stage 1: $(core.pipeline.stage1)\n", color = color)

    # Select a color for Stage 2
    if stall_signals.stall_234
        color = :red
    elseif stall_signals.nop_2
        color = :cyan
    else
        color = :white
    end
    printstyled("Stage 2: $(core.pipeline.stage2)\n")

    # Select a color for stages 3 and 4.
    color = stall_signals.stall_234 ? (:red) : (:white)
    printstyled("Stage 3: $(core.pipeline.stage3)\n", color = color)
    printstyled("Stage 4: $(core.pipeline.stage4)\n", color = color)

    # Print stage 5 with its myriad of colors.
    if stall_signals.nop_5
        color = :cyan
    elseif stall_signals.stall_5678
        color = :red
    else
        color = :white
    end
    printstyled(color, "Stage 5: $(core.pipeline.stage5)\n", color = color)

    # Default to Stalled and Unstalled.
    color = stall_signals.stall_5678  ? (:red) : (:white)
    printstyled("Stage 6: $(core.pipeline.stage6)\n", color = color)
    printstyled("Stage 7: $(core.pipeline.stage7)\n", color = color)
end

ext: jl


text:
#imports

using JuMP, Gurobi, DelimitedFiles, Random

function atualize_lowerBound(A,c,x_res,indCost)

    (m,n) = size(A)
    I = (1:m)
    J = (1:n)
    x_resLower = zeros(n,1)
    global x_res_temp = zeros(n,1)

    for indice_j in indCost
        indice_j = floor(Int,indice_j)
        if x_res[indice_j] == 1
            x_res_temp[indice_j] = 1
            value_mult = A*x_res_temp

            flag = 0
            for i in I
                if value_mult[i] > 1

                    flag = 1
                    #break
                end
            end

            if flag == 0
                x_resLower[indice_j] = 1
            end

            global x_res_temp = copy(x_resLower)

        end
    end

    return c'*x_resLower

end

function RelaxLagrang(A, c, maxiter = 50, tol = 1e-1)
    spp   = Model(Gurobi.Optimizer);
    (m,n) = size(A)

    I = (1:m)
    J = (1:n)

    @variable(spp, x[1:n], Bin);
    u     = 1*ones(m,1);
    zs = []
    irep    = 0;
    iter    = 0;
    lower_b = 0;
    upper_b = Inf;

    step_size = 0
    grad      = zeros(m,1);
    indCost   = zeros(1,n)

    for j1 in J
        maior_valor = -Inf
        for j2 in J
            if (c[j2] >=  maior_valor) && ((j2 in indCost) == false)
                global maior_indice = j2
                 maior_valor = c[j2]
            end
        end
        indCost[j1] = maior_indice
    end

    while iter <= maxiter && abs(upper_b - lower_b) >= tol

        iter += 1
        u = max.(0, u - step_size*grad)

        @objective(spp, Max, sum(c.*x) - sum(u[i]* (sum(A[i,j]*x[j] for j in (1:n)) - 1) for i in (1:m)))
        optimize!(spp)

        x_res    = value.(x)
        soma_res = A*x_res
        z        = objective_value(spp);
        append!(zs, z)

        best_ub_old = upper_b;
        upper_b = min(upper_b, z);

        if best_ub_old == upper_b
            irep += 1;
        else
            irep = 0;
            lower_b = atualize_lowerBound(A,c,x_res,indCost)[1];

        end

        grad = ones(m,1) - soma_res + ones(m,1)*rand(1)/10
        sum_g2 = sum(grad[i]^2 for i in I)
        #step_size = (z - lower_b)/sum_g2
        step_size = 1/iter
        #step_size  = 0.008

    end

    return zs, iter
end

function loadSPP(fname)
    f=open(fname)

    m, n = parse.(Int, split(readline(f)))
    C = parse.(Int, split(readline(f)))
    A=zeros(Int, m, n)
    for i=1:m
        readline(f)
        for valor in split(readline(f))
          j = parse(Int, valor)
          A[i,j]=1
        end
    end
    close(f)
    return C, A
end

function setSPP(solverSelected, C, A)
  m, n = size(A)
  ip = Model(solverSelected)
  @variable(ip, x[1:n], Bin)
  @objective(ip, Max, sum(C.*x))
  @constraint(ip , cte[i=1:m], sum(A[i,j] * x[j] for j=1:n) <= 1)
  return ip, x
end

function setSPP_relaxlin(solverSelected, C, A)
  m, n = size(A)
  ip = Model(solverSelected)
  @variable(ip, 0 <= x[1:n] <= 1)
  @objective(ip, Max, sum(C.*x))
  @constraint(ip , cte[i=1:m], sum(A[i,j] * x[j] for j=1:n) <= 1)
  return ip, x
end

function integer_test(x)
    if floor(x)^2 - x^2 == 0
        return 1
    end
    return 0
end

function branch_ceil(solverSelected, C, A, idxs = [])
  m, n = size(A)
  ip = Model(solverSelected)
  @variable(ip, 0 <= x[1:n] <= 1)
  for idx in idxs
    @constraint(ip, x[idx] == 1)
  end
  @objective(ip, Max, sum(C.*x))
  @constraint(ip , cte[i=1:m], sum(A[i,j] * x[j] for j=1:n) <= 1)
  return ip, x
end

function branch_floor(solverSelected, C, A, idxs = [])
  m, n = size(A)
  ip = Model(solverSelected)
  @variable(ip, 0 <= x[1:n] <= 1)
  for idx in idxs
    @constraint(ip, x[idx] == 0)
  end
  @objective(ip, Max, sum(C.*x))
  @constraint(ip , cte[i=1:m], sum(A[i,j] * x[j] for j=1:n) <= 1)
  return ip, x
end

function BranchNBound(A, c, maxiter = 50, tol = 1e-1)

    rp, rp_x = setSPP_relaxlin(Gurobi.Optimizer, c, A)
    nodes    = [[rp,[]]]
    s,       = size(nodes)
    values   = []
    idxs     = []

    while s >= 0

        if s == 0
            break
        end

        p,idxs   = nodes[1]

        if s == 1
            nodes = []
            s     = 0
        else
            nodes = nodes[2:s]
            s,     = size(nodes)
        end

        try
           optimize!(p)
        catch
           println(" Modelo Infeasible. Cortem as árvores! ")
           ss,       = size(idxs)
           idxs      = idxs[1:ss-1]

        end
        append!(values, objective_value(p))

        for (i,xs) in enumerate(value.((all_variables(p))))
            if integer_test(xs) == 0

                push!(idxs, i)
                print("----------> ", idxs, " <---------")
                print(" BRANCHING ")

                floor_prob, f_x = branch_floor(Gurobi.Optimizer, c, A,idxs)
                push!(nodes, floor_prob)

                ceil_prob, c_x  =  branch_ceil(Gurobi.Optimizer, c, A,idxs)
                push!(nodes, ceil_prob)

                s += 2
                print("--> ", values, " <---")
                break
            end
        end
    end

    return objective_value(rp), value.(rp_x), values
end

function Branch_Bound(A,c)

    rp, rp_x = setSPP(Gurobi.Optimizer, c, A)

    println("Solving...");
    set_optimizer_attribute(rp, "Presolve", 0)
    set_optimizer_attribute(rp, "Cuts", 0)
    set_optimizer_attribute(rp, "Heuristics", 0)
    optimize!(rp)
end

function Branch_Cut(A,c)

    p, p_x = setSPP(Gurobi.Optimizer, c, A)

    println("Solving...");
    set_optimizer_attribute(p, "Presolve", 0)
    set_optimizer_attribute(p, "Cuts", 0)
    set_optimizer_attribute(p, "Heuristics", 0)
    set_optimizer_attribute(p, "CliqueCuts", 1)
    optimize!(p)

end


text:
using JetPack, Jets, Test

n1,n2,n3 = 11,22,33

@testset "JopGradient2D, linearity test, T=$T" for T in (Float32,Float64)
    A = JopGradient(JetSpace(T,n1,n2), (rand(T),rand(T)))
    lhs,rhs = linearity_test(A)
    @test lhs ≈ rhs
end

@testset "Gradient3D, linearity test, T=$(T)" for T in (Float32,Float64)
    A = JopGradient(JetSpace(T,n1,n2,n3), (rand(T),rand(T),rand(T)))
    lhs,rhs = linearity_test(A)
    @test lhs ≈ rhs
end

@testset "Gradient2D, dot product test, T=$(T)" for T in (Float32,Float64)
    A = JopGradient(JetSpace(T,n1,n2), (rand(T),rand(T)))
    lhs,rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs rtol=1e-4
end

@testset "Gradient3D dot product test, T=$(T)" for T in (Float64,Float32)
    T=Float64
    A = JopGradient(JetSpace(T,n1,n2,n3), (rand(T),rand(T),rand(T)))
    lhs,rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs rtol=1e-4
end

@testset "Gradient2D correctness test, T=$(T)" for T in (Float64,Float32)
    dom = JetSpace(T,n1,n2)
    delta = (rand(T), rand(T))
    op = JopGradient(dom, delta)
    x = zeros(domain(op))
    α = rand(T)
    β = rand(T)
    for k1 = 1:size(x,1)
        for k2 = 1:size(x,2)
            x[k1,k2] = α * k1 + β * k2
        end
    end
    y = op * x
    diff1 = 0.0
    diff2 = 0.0
    for k1 = 1:size(x,1)
        for k2 = 1:size(x,2)
            if k1 < n1
                diff1 += α - y[k1,k2,1]
            end
            if k2 < n2
                diff2 += β - y[k1,k2,2]
            end
        end
    end
    diff1 /= length(x)
    diff2 /= length(x)
    @test diff1 < 1000 * eps(T)
    @test diff2 < 1000 * eps(T)
end

@testset "Gradient3D correctness test, T=$(T)" for T in (Float64,Float32)
    dom = JetSpace(T,n1,n2,n3)
    delta = (rand(T), rand(T), rand(T))
    op = JopGradient(dom, delta)
    x = zeros(domain(op))
    α = rand(T)
    β = rand(T)
    γ = rand(T)
    for k1 = 1:size(x,1)
        for k2 = 1:size(x,2)
            for k3 = 1:size(x,3)
                x[k1,k2,k3] = α * k1 + β * k2 + γ * k3
            end
        end
    end
    y = op * x
    diff1 = 0.0
    diff2 = 0.0
    diff3 = 0.0
    for k1 = 1:size(x,1)
        for k2 = 1:size(x,2)
            for k3 = 1:size(x,3)
                if k1 < n1
                    diff1 += α - y[k1,k2,k3,1]
                end
                if k2 < n2
                    diff2 += β - y[k1,k2,k3,2]
                end
                if k3 < n3
                    diff3 += γ - y[k1,k2,k3,3]
                end
            end
        end
    end
    diff1 /= length(x)
    diff2 /= length(x)
    diff3 /= length(x)
    @test diff1 < 1000 * eps(T)
    @test diff2 < 1000 * eps(T)
end


text:
@testset "format with positional arguments" begin

    @test format("{1}", 10) == "10"
    @test format("αβγ {1}", 10) == "αβγ 10"
    @test format("{1} efg", 10) == "10 efg"
    @test format("{1} ϵζη", 10) == "10 ϵζη"
    @test format("abc {1} efg", 10) == "abc 10 efg"
    @test format("αβγ{1}ϵζη", 10) == "αβγ10ϵζη"
    @test format("{1} + {2}", 10, "xyz") == "10 + xyz"
    @test format("{1} + {2}", 10, "χψω") == "10 + χψω"
    @test format("abc {1} + {2}", 10, "xyz") == "abc 10 + xyz"
    @test format("αβγ {1} + {2}", 10, "χψω") == "αβγ 10 + χψω"
    @test format("{1} + {2} efg", 10, "xyz") == "10 + xyz efg"
    @test format("{1} + {2} ϵζη", 10, "χψω") == "10 + χψω ϵζη"
    @test format("abc {1} + {2} efg", 10, "xyz") == "abc 10 + xyz efg"
    @test format("αβγ {1} + {2} ϵζη", 10, "χψω") == "αβγ 10 + χψω ϵζη"
    @test format("abc {1}{2} efg", 10, "xyz") == "abc 10xyz efg"
    @test format("αβγ {1}{2} ϵζη", 10, "χψω") == "αβγ 10χψω ϵζη"

    @test format("{1:d} + {2:s}", 10, "xyz") == "10 + xyz"
    @test format("{1:d} + {2:s}", 10, "χψω") == "10 + χψω"
    @test format("{1:04d} + {2:*>5}", 10, "xyz") == "0010 + **xyz"
    @test format("{1:04d} + {2:⋆>5}", 10, "χψω") == "0010 + ⋆⋆χψω"
    @test format("let {2:<5} := {1:.4f};", 12.3, "var") == "let var   := 12.3000;"
    @test format("let {2:<5} := {1:.4f};", 12.3, "χψω") == "let χψω   := 12.3000;"

    @test format("{}", 10) == "10"
    @test format("{} + {}", 10, 20) == "10 + 20"
    @test format("{} + {:04d}", 10, 20) == "10 + 0020"
    @test format("{:03d} + {}", 10, 20) == "010 + 20"
    @test format("{:03d} + {:04d}", 10, 20) == "010 + 0020"
    @test_throws(ErrorException, format("{1} + {}", 10, 20) )
    @test_throws(ErrorException, format("{} + {1}", 10, 20) )
    @test_throws(ErrorException, format("{1", 10) )

    f = FormatExpr("{1:s}")
    printfmt(f,"")
    printfmtln(f,"")
    io = IOBuffer()
    printfmt(io,f,"foobar")
    @test io.size == 6
    printfmtln(io,f,"ugly")
    @test io.size == 11
end

@testset "format escape {{ and }}" begin
    @test format("{{}}") == "{}"
    @test format("{{{1}}}", 10) == "{10}"
    @test format("v: {{{2}}} = {1:.4f}", 1.2, "ab") == "v: {ab} = 1.2000"
    @test format("χ: {{{2}}} = {1:.4f}", 1.2, "αβ") == "χ: {αβ} = 1.2000"
end

@testset "format with filter" begin
    @test format("{1|>abs2} + {2|>abs2:.2f}", 2, 3) == "4 + 9.00"
end


text:
using DataValues
using LinearAlgebra
using Test

@testset "DataValueArray: DataMatrix" begin

#----- test Base.diag -----#
A = reshape([1:25...], 5, 5)
m = fill(false, 5, 5)
m[1] = true
m[25] = true
X = DataValueArray(A)
Y = DataValueArray(A, m)

@test isequal(diag(X), DataValueArray([1, 7, 13, 19, 25]))
@test isequal(diag(Y), DataValueArray([1, 7, 13, 19, 25],
                                        [true, false, false, false, true]))

end


text:
# Spawn a thread (using pthreads on Unix/OSX and Windows threads on Windows)
# to implement the "heartbeat" message channel, which is just a ZMQ
# socket that echoes all messages.  This is implemented with the zmq_proxy
# call in libzmq, which simply blocks forever, so the usual lack of
# thread safety in Julia should not be an issue here.

import Libdl

const threadid = zeros(Int, 128) # sizeof(uv_thread_t) <= 8 on Linux, OSX, Win
const zmq_proxy = Ref(C_NULL)

# entry point for new thread
function heartbeat_thread(sock::Ptr{Cvoid})
    ccall(zmq_proxy[], Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
          sock, sock, C_NULL)
    nothing
end

function start_heartbeat(sock)
    zmq_proxy[] = Libdl.dlsym(Libdl.dlopen(ZMQ.libzmq), :zmq_proxy)
    heartbeat_c = @cfunction(heartbeat_thread, Cvoid, (Ptr{Cvoid},))
    ccall(:uv_thread_create, Cint, (Ptr{Int}, Ptr{Cvoid}, Ptr{Cvoid}),
          threadid, heartbeat_c, sock)
end


text:
function invoke(f::Function, pos::Int, ::Type{Float64}, args...)
    out = Ref{Float64}(0.0) # should we use Cfloat here instead?

    args = collect(args)
    insert!(args, pos-1, out)

    r = f(args...)
    r != 0 ? throw(XpressError(r, "Unable to invoke $f")) : out[]
end

function invoke(f::Function, pos::Int, ::Type{Int}, args...)
    out = Ref{Cint}(0)

    args = collect(args)
    insert!(args, pos-1, out)

    r = f(args...)
    r != 0 ? throw(XpressError(r, "Unable to invoke $f")) : out[]
end

function invoke(f::Function, pos::Int, ::Type{String}, args...)
    out = Cstring(pointer(Array{Cchar}(undef, 1024)))

    args = collect(Any, args)
    insert!(args, pos-1, out)

    r = f(args...)
    r != 0 ? throw(XpressError(r, "Unable to invoke $f")) : unsafe_string(out)
end

"""
    @invoke expr

Lets you invoke a lower level `Lib` function.
Xpress' library API expects the caller to pre-allocate memory.
Use this macro to minimize repetition and increase readability.

This macro expects syntax that mimics the call to the `Lib.function`.
The `_` argument must be used for all arguments where the `Lib.function` expects the caller to manage memory.
Additionally, the return type declaration must be used.

Examples:

    @invoke Lib.XPRSgetversion(_)::String
    @invoke Lib.XPRSgetbanner(_)::String
    @invoke Lib.XPRSgetprobname(prob, _)::String

As an example of what @invoke expands to:

```
julia> @macroexpand @invoke Lib.XPRSgetversion(_)::String
:(invoke(Lib.XPRSgetversion, Xpress.String))
```

"""
macro invoke(expr)
    @assert expr.head == :(::) "macro argument must have return type declaration"

    # macro return type must be a valid type that exists in Xpress or Julia
    return_type = expr.args[2]

    f = expr.args[1]

    @assert f.head == :call "macro argument must have a function call"
    @assert :_ in f.args[2:end] "macro argument must have an underscore argument"

    # Remove all `_` arguments
    # TODO: note the positions and pass positions to the invoke function
    indices = findall(x->x==:_, f.args)
    filter!(x->x≠:_, f.args)

    # Call invoke function at macro call site instead
    pushfirst!(f.args, :(invoke))

    f.args = esc.(f.args)

    # invoke function takes the position of return type as the first argument
    # invoke function takes the return type as the first argument
    # invoke function uses this argument to dispatch to the correct method
    if length(indices) == 1
        insert!(f.args, 3, return_type)
        insert!(f.args, 3, indices[1])
    else
        # TODO: Implement multiple return types
        error("Not implemented @invoke macro for multiple `_`")
    end

    return f
end

macro checked(expr)
    @assert expr.head == :call "Can only use @checked on function calls"
    @assert ( expr.args[1].head == :(.) ) && ( expr.args[1].args[1] == :Lib) "Can only use @checked on Lib.\$function"
    f = replace(String(expr.args[1].args[2].value), "XPRS"=>"")
    return esc(quote
        r = $(expr)
        r != 0 ? throw(XpressError(r, "Unable to call $($f)")) : nothing
    end)
end


text:
const e0_coefficent = (

# 1-10
( 0,  0,  0,  0,  1,  0,  0,  0),
( 0,  0,  0,  0,  2,  0,  0,  0),
( 0,  0,  2, -2,  3,  0,  0,  0),
( 0,  0,  2, -2,  1,  0,  0,  0),
( 0,  0,  2, -2,  2,  0,  0,  0),
( 0,  0,  2,  0,  3,  0,  0,  0),
( 0,  0,  2,  0,  1,  0,  0,  0),
( 0,  0,  0,  0,  3,  0,  0,  0),
( 0,  1,  0,  0,  1,  0,  0,  0),
( 0,  1,  0,  0, -1,  0,  0,  0),

# 11-20
( 1,  0,  0,  0, -1,  0,  0,  0),
( 1,  0,  0,  0,  1,  0,  0,  0),
( 0,  1,  2, -2,  3,  0,  0,  0),
( 0,  1,  2, -2,  1,  0,  0,  0),
( 0,  0,  4, -4,  4,  0,  0,  0),
( 0,  0,  1, -1,  1, -8, 12,  0),
( 0,  0,  2,  0,  0,  0,  0,  0),
( 0,  0,  2,  0,  2,  0,  0,  0),
( 1,  0,  2,  0,  3,  0,  0,  0),
( 1,  0,  2,  0,  1,  0,  0,  0),

# 21-30
( 0,  0,  2, -2,  0,  0,  0,  0),
( 0,  1, -2,  2, -3,  0,  0,  0),
( 0,  1, -2,  2, -1,  0,  0,  0),
( 0,  0,  0,  0,  0,  8,-13, -1),
( 0,  0,  0,  2,  0,  0,  0,  0),
( 2,  0, -2,  0, -1,  0,  0,  0),
( 1,  0,  0, -2,  1,  0,  0,  0),
( 0,  1,  2, -2,  2,  0,  0,  0),
( 1,  0,  0, -2, -1,  0,  0,  0),
( 0,  0,  4, -2,  4,  0,  0,  0),

# 31-33
( 0,  0,  2, -2,  4,  0,  0,  0),
( 1,  0, -2,  0, -3,  0,  0,  0),
( 1,  0, -2,  0, -1,  0,  0,  0)
)

const e0_arg = (

    (2640.96e-6, -0.39e-6 ),
    (  63.52e-6, -0.02e-6 ),
    (  11.75e-6,  0.01e-6 ),
    (  11.21e-6,  0.01e-6 ),
    (  -4.55e-6,  0.00e-6 ),
    (   2.02e-6,  0.00e-6 ),
    (   1.98e-6,  0.00e-6 ),
    (  -1.72e-6,  0.00e-6 ),
    (  -1.41e-6, -0.01e-6 ),
    (  -1.26e-6, -0.01e-6 ),
    (  -0.63e-6,  0.00e-6 ),
    (  -0.63e-6,  0.00e-6 ),
    (   0.46e-6,  0.00e-6 ),
    (   0.45e-6,  0.00e-6 ),
    (   0.36e-6,  0.00e-6 ),
    (  -0.24e-6, -0.12e-6 ),
    (   0.32e-6,  0.00e-6 ),
    (   0.28e-6,  0.00e-6 ),
    (   0.27e-6,  0.00e-6 ),
    (   0.26e-6,  0.00e-6 ),
    (  -0.21e-6,  0.00e-6 ),
    (   0.19e-6,  0.00e-6 ),
    (   0.18e-6,  0.00e-6 ),
    (  -0.10e-6,  0.05e-6 ),
    (   0.15e-6,  0.00e-6 ),
    (  -0.14e-6,  0.00e-6 ),
    (   0.14e-6,  0.00e-6 ),
    (  -0.14e-6,  0.00e-6 ),
    (   0.14e-6,  0.00e-6 ),
    (   0.13e-6,  0.00e-6 ),
    (  -0.11e-6,  0.00e-6 ),
    (   0.11e-6,  0.00e-6 ),
    (   0.11e-6,  0.00e-6 )

)
const e1_coefficent = (
( 0,  0,  0,  0,  1,  0,  0,  0),
)
const e1_arg = (
(-0.87e-6,  0.00e-6 ),)


text:
# e_ident[EIBCLASS]
const ELFCLASSNONE = 0 #Invalid class
const ELFCLASS32 = 1 #32-bit objects
const ELFCLASS64 = 2 #64-bit objects

# e_ident[EI_DATA]
const ELFDATANONE = 0
const ELFDATA2LSB = 1
const ELFDATA2MSB = 2

# e_ident[EI_VERSION]
const EV_NONE = 0       #Invalid Version
const EV_CURRENT = 1    #Current Version

# e_ident[EI_OSABI]
const ELFOSABI_NONE         = 0   # No extensions or unspecified
const ELFOSABI_HPUX         = 1   # Hewlett-Packard HP-UX
const ELFOSABI_NETBSD       = 2   # NetBSD
const ELFOSABI_GNU          = 3   # GNU
const ELFOSABI_LINUX        = 3   # Linux historical - alias for ELFOSABI_GNU
const ELFOSABI_SOLARIS      = 6   # Sun Solaris
const ELFOSABI_AIX          = 7   # AIX
const ELFOSABI_IRIX         = 8   # IRIX
const ELFOSABI_FREEBSD      = 9   # FreeBSD
const ELFOSABI_TRU64        = 10  # Compaq TRU64 UNIX
const ELFOSABI_MODESTO      = 11  # Novell Modesto
const ELFOSABI_OPENBSD      = 12  # Open BSD
const ELFOSABI_OPENVMS      = 13  # Open VMS
const ELFOSABI_NSK          = 14  # Hewlett-Packard Non-Stop Kernel
const ELFOSABI_AROS         = 15  # Amiga Research OS
const ELFOSABI_FENIXOS      = 16  # FenixOS
const ELFOSABI_CLOUDABI     = 17  # Nuxi CloudABI
const ELFOSABI_OPENVOS      = 18  # Stratus Technologies OpenVOS
const ELFOSABI_C6000_ELFABI = 64  # Bare-metal TMS320C6000
const ELFOSABI_C6000_LINUX  = 65  # Linux TMS320C6000
const ELFOSABI_ARM	        = 97  # ARM
const ELFOSABI_STANDALONE   = 255 # Standalone (embedded) application

# e_type constants
@constants ET_TYPES "ET_" begin
    const ET_NONE = 0   #No file type
    const ET_REL  = 1   #Relocatable file
    const ET_EXEC = 2   #Executable file
    const ET_DYN  = 3   #Shared object file
    const ET_CORE = 4   #Core file
end
@constants ET_RANGES "ET_" begin
    const ET_LOOS       = 0xfe00  #Operating system-specific
    const ET_HIOS       = 0xfeff  #Operating system-specific
    const ET_LOPROC     = 0xff00  #Processor-specific
    const ET_HIPROC     = 0xffff  #Processor-specific
end

@constants EM_MACHINES "EM_" begin
    # e_machine constants
    const EM_NONE = 0   #No machine
    const EM_M32  = 1   #AT&T WE 32100
    const EM_SPARC    = 2   #SPARC
    const EM_386  = 3   #Intel 80386
    const EM_68K  = 4   #Motorola 68000
    const EM_88K  = 5   #Motorola 88000
    #6                  Reserved for future use (was EM_486)
    const EM_860  = 7   #Intel 80860
    const EM_MIPS = 8   #MIPS I Architecture
    const EM_S370 = 9   #IBM System/370 Processor
    const EM_MIPS_RS3_LE  = 10  #MIPS RS3000 Little-endian
    #11-14              Reserved for future use
    const EM_PARISC   = 15  #Hewlett-Packard PA-RISC
    #16  Reserved for future use
    const EM_VPP500   = 17  #Fujitsu VPP500
    const EM_SPARC32PLUS  = 18  #Enhanced instruction set SPARC
    const EM_960  = 19  #Intel 80960
    const EM_PPC  = 20  #PowerPC
    const EM_PPC64    = 21  #64-bit PowerPC
    const EM_S390 = 22  #IBM System/390 Processor
    const EM_SPU  = 23  #IBM SPU/SPC
    #24-35   Reserved for future use
    const EM_V800 = 36  #NEC V800
    const EM_FR20 = 37  #Fujitsu FR20
    const EM_RH32 = 38  #TRW RH-32
    const EM_RCE  = 39  #Motorola RCE
    const EM_ARM  = 40  #ARM 32-bit architecture (AARCH32)
    const EM_ALPHA    = 41  #Digital Alpha
    const EM_SH   = 42  #Hitachi SH
    const EM_SPARCV9  = 43  #SPARC Version 9
    const EM_TRICORE  = 44  #Siemens TriCore embedded processor
    const EM_ARC  = 45  #Argonaut RISC Core, Argonaut Technologies Inc.
    const EM_H8_300   = 46  #Hitachi H8/300
    const EM_H8_300H  = 47  #Hitachi H8/300H
    const EM_H8S  = 48  #Hitachi H8S
    const EM_H8_500   = 49  #Hitachi H8/500
    const EM_IA_64    = 50  #Intel IA-64 processor architecture
    const EM_MIPS_X   = 51  #Stanford MIPS-X
    const EM_COLDFIRE = 52  #Motorola ColdFire
    const EM_68HC12   = 53  #Motorola M68HC12
    const EM_MMA  = 54  #Fujitsu MMA Multimedia Accelerator
    const EM_PCP  = 55  #Siemens PCP
    const EM_NCPU = 56  #Sony nCPU embedded RISC processor
    const EM_NDR1 = 57  #Denso NDR1 microprocessor
    const EM_STARCORE = 58  #Motorola Star*Core processor
    const EM_ME16 = 59  #Toyota ME16 processor
    const EM_ST100    = 60  #STMicroelectronics ST100 processor
    const EM_TINYJ    = 61  #Advanced Logic Corp. TinyJ embedded processor family
    const EM_X86_64   = 62  #AMD x86-64 architecture
    const EM_PDSP = 63  #Sony DSP Processor
    const EM_PDP10    = 64  #Digital Equipment Corp. PDP-10
    const EM_PDP11    = 65  #Digital Equipment Corp. PDP-11
    const EM_FX66 = 66  #Siemens FX66 microcontroller
    const EM_ST9PLUS  = 67  #STMicroelectronics ST9+ 8/16 bit microcontroller
    const EM_ST7  = 68  #STMicroelectronics ST7 8-bit microcontroller
    const EM_68HC16   = 69  #Motorola MC68HC16 Microcontroller
    const EM_68HC11   = 70  #Motorola MC68HC11 Microcontroller
    const EM_68HC08   = 71  #Motorola MC68HC08 Microcontroller
    const EM_68HC05   = 72  #Motorola MC68HC05 Microcontroller
    const EM_SVX  = 73  #Silicon Graphics SVx
    const EM_ST19 = 74  #STMicroelectronics ST19 8-bit microcontroller
    const EM_VAX  = 75  #Digital VAX
    const EM_CRIS = 76  #Axis Communications 32-bit embedded processor
    const EM_JAVELIN  = 77  #Infineon Technologies 32-bit embedded processor
    const EM_FIREPATH = 78  #Element 14 64-bit DSP Processor
    const EM_ZSP  = 79  #LSI Logic 16-bit DSP Processor
    const EM_MMIX = 80  #Donald Knuth's educational 64-bit processor
    const EM_HUANY    = 81  #Harvard University machine-independent object files
    const EM_PRISM    = 82  #SiTera Prism
    const EM_AVR  = 83  #Atmel AVR 8-bit microcontroller
    const EM_FR30 = 84  #Fujitsu FR30
    const EM_D10V = 85  #Mitsubishi D10V
    const EM_D30V = 86  #Mitsubishi D30V
    const EM_V850 = 87  #NEC v850
    const EM_M32R = 88  #Mitsubishi M32R
    const EM_MN10300  = 89  #Matsushita MN10300
    const EM_MN10200  = 90  #Matsushita MN10200
    const EM_PJ   = 91  #picoJava
    const EM_OPENRISC = 92  #OpenRISC 32-bit embedded processor
    const EM_ARC_COMPACT  = 93  #ARC International ARCompact processor (old spelling/synonym: EM_ARC_A5)
    const EM_XTENSA   = 94  #Tensilica Xtensa Architecture
    const EM_VIDEOCORE    = 95  #Alphamosaic VideoCore processor
    const EM_TMM_GPP  = 96  #Thompson Multimedia General Purpose Processor
    const EM_NS32K    = 97  #National Semiconductor 32000 series
    const EM_TPC  = 98  #Tenor Network TPC processor
    const EM_SNP1K    = 99  #Trebia SNP 1000 processor
    const EM_ST200    = 100 #STMicroelectronics (www.st.com) ST200 microcontroller
    const EM_IP2K = 101 #Ubicom IP2xxx microcontroller family
    const EM_MAX  = 102 #MAX Processor
    const EM_CR   = 103 #National Semiconductor CompactRISC microprocessor
    const EM_F2MC16   = 104 #Fujitsu F2MC16
    const EM_MSP430   = 105 #Texas Instruments embedded microcontroller msp430
    const EM_BLACKFIN = 106 #Analog Devices Blackfin (DSP) processor
    const EM_SE_C33   = 107 #S1C33 Family of Seiko Epson processors
    const EM_SEP  = 108 #Sharp embedded microprocessor
    const EM_ARCA = 109 #Arca RISC Microprocessor
    const EM_UNICORE  = 110 #Microprocessor series from PKU-Unity Ltd. and MPRC of Peking University
    const EM_EXCESS   = 111 #eXcess: 16/32/64-bit configurable embedded CPU
    const EM_DXP  = 112 #Icera Semiconductor Inc. Deep Execution Processor
    const EM_ALTERA_NIOS2 = 113 #Altera Nios II soft-core processor
    const EM_CRX  = 114 #National Semiconductor CompactRISC CRX microprocessor
    const EM_XGATE    = 115 #Motorola XGATE embedded processor
    const EM_C166 = 116 #Infineon C16x/XC16x processor
    const EM_M16C = 117 #Renesas M16C series microprocessors
    const EM_DSPIC30F = 118 #Microchip Technology dsPIC30F Digital Signal Controller
    const EM_CE   = 119 #Freescale Communication Engine RISC core
    const EM_M32C = 120 #Renesas M32C series microprocessors
    #121-130 Reserved for future use
    const EM_TSK3000  = 131 #Altium TSK3000 core
    const EM_RS08 = 132 #Freescale RS08 embedded processor
    const EM_SHARC    = 133 #Analog Devices SHARC family of 32-bit DSP processors
    const EM_ECOG2    = 134 #Cyan Technology eCOG2 microprocessor
    const EM_SCORE7   = 135 #Sunplus S+core7 RISC processor
    const EM_DSP24    = 136 #New Japan Radio (NJR) 24-bit DSP Processor
    const EM_VIDEOCORE3   = 137 #Broadcom VideoCore III processor
    const EM_LATTICEMICO32    = 138 #RISC processor for Lattice FPGA architecture
    const EM_SE_C17   = 139 #Seiko Epson C17 family
    const EM_TI_C6000 = 140 #The Texas Instruments TMS320C6000 DSP family
    const EM_TI_C2000 = 141 #The Texas Instruments TMS320C2000 DSP family
    const EM_TI_C5500 = 142 #The Texas Instruments TMS320C55x DSP family
    #143-159 Reserved for future use
    const EM_MMDSP_PLUS   = 160 #STMicroelectronics 64bit VLIW Data Signal Processor
    const EM_CYPRESS_M8C  = 161 #Cypress M8C microprocessor
    const EM_R32C = 162 #Renesas R32C series microprocessors
    const EM_TRIMEDIA = 163 #NXP Semiconductors TriMedia architecture family
    const EM_QDSP6    = 164 #QUALCOMM DSP6 Processor
    const EM_8051 = 165 #Intel 8051 and variants
    const EM_STXP7X   = 166 #STMicroelectronics STxP7x family of configurable and extensible RISC processors
    const EM_NDS32    = 167 #Andes Technology compact code size embedded RISC processor family
    const EM_ECOG1    = 168 #Cyan Technology eCOG1X family
    const EM_ECOG1X   = 168 #Cyan Technology eCOG1X family
    const EM_MAXQ30   = 169 #Dallas Semiconductor MAXQ30 Core Micro-controllers
    const EM_XIMO16   = 170 #New Japan Radio (NJR) 16-bit DSP Processor
    const EM_MANIK    = 171 #M2000 Reconfigurable RISC Microprocessor
    const EM_CRAYNV2  = 172 #Cray Inc. NV2 vector architecture
    const EM_RX   = 173 #Renesas RX family
    const EM_METAG    = 174 #Imagination Technologies META processor architecture
    const EM_MCST_ELBRUS  = 175 #MCST Elbrus general purpose hardware architecture
    const EM_ECOG16   = 176 #Cyan Technology eCOG16 family
    const EM_CR16 = 177 #National Semiconductor CompactRISC CR16 16-bit microprocessor
    const EM_ETPU = 178 #Freescale Extended Time Processing Unit
    const EM_SLE9X    = 179 #Infineon Technologies SLE9X core
    const EM_L10M = 180 #Intel L10M
    const EM_K10M = 181 #Intel K10M
    #182 Reserved for future Intel use
    const EM_AARCH64  = 183 #ARM 64-bit architecture (AARCH64)
    #184 Reserved for future ARM use
    const EM_AVR32    = 185 #Atmel Corporation 32-bit microprocessor family
    const EM_STM8 = 186 #STMicroeletronics STM8 8-bit microcontroller
    const EM_TILE64   = 187 #Tilera TILE64 multicore architecture family
    const EM_TILEPRO  = 188 #Tilera TILEPro multicore architecture family
    const EM_MICROBLAZE   = 189 #Xilinx MicroBlaze 32-bit RISC soft processor core
    const EM_CUDA = 190 #NVIDIA CUDA architecture
    const EM_TILEGX   = 191 #Tilera TILE-Gx multicore architecture family
    const EM_CLOUDSHIELD  = 192 #CloudShield architecture family
    const EM_COREA_1ST    = 193 #KIPO-KAIST Core-A 1st generation processor family
    const EM_COREA_2ND    = 194 #KIPO-KAIST Core-A 2nd generation processor family
    const EM_ARC_COMPACT2 = 195 #Synopsys ARCompact V2
    const EM_OPEN8    = 196 #Open8 8-bit RISC soft processor core
    const EM_RL78 = 197 #Renesas RL78 family
    const EM_VIDEOCORE5   = 198 #Broadcom VideoCore V processor
    const EM_78KOR    = 199 #Renesas 78KOR family
    const EM_56800EX  = 200 #Freescale 56800EX Digital Signal Controller (DSC)
    const EM_BA1  = 201 #Beyond BA1 CPU architecture
    const EM_BA2  = 202 #Beyond BA2 CPU architecture
    const EM_XCORE    = 203 #XMOS xCORE processor family
    const EM_MCHP_PIC = 204 #Microchip 8-bit PIC(r) family
end

# Special Section Indices
const SHN_UNDEF   = 0
const SHN_LORESERVE   = 0xff00
const SHN_LOPROC  = 0xff00
const SHN_HIPROC  = 0xff1f
const SHN_LOOS    = 0xff20
const SHN_HIOS    = 0xff3f
const SHN_ABS     = 0xfff1
const SHN_COMMON  = 0xfff2
const SHN_XINDEX  = 0xffff
const SHN_HIRESERVE   = 0xffff

# sh_type constants
@constants SHT_TYPES "SHT_" begin
    const SHT_NULL          =  0
    const SHT_PROGBITS      =  1
    const SHT_SYMTAB        =  2
    const SHT_STRTAB        =  3
    const SHT_RELA          =  4
    const SHT_HASH          =  5
    const SHT_DYNAMIC       =  6
    const SHT_NOTE          =  7
    const SHT_NOBITS        =  8
    const SHT_REL           =  9
    const SHT_SHLIB         = 10
    const SHT_DYNSYM        = 11
    const SHT_INIT_ARRAY    = 14
    const SHT_FINI_ARRAY    = 15
    const SHT_PREINIT_ARRAY = 16
    const SHT_GROUP         = 17
    const SHT_SYMTAB_SHNDX  = 18
end
@constants SHT_MASKS "SHT_" begin
    const SHT_LOOS          = 0x60000000
    const SHT_HIOS          = 0x6fffffff
    const SHT_LOPROC        = 0x70000000
    const SHT_HIPROC        = 0x7fffffff
    const SHT_LOUSER        = 0x80000000
    const SHT_HIUSER        = 0xffffffff
end

# Section Attribute Flags
@constants SHF_FLAGS "SHF_" begin
    const SHF_WRITE                 =   0x1
    const SHF_ALLOC                 =   0x2
    const SHF_EXECINSTR             =   0x4
    const SHF_MERGE                 =  0x10
    const SHF_STRINGS               =  0x20
    const SHF_INFO_LINK             =  0x40
    const SHF_LINK_ORDER            =  0x80
    const SHF_OS_NONCONFORMING      = 0x100
    const SHF_GROUP                 = 0x200
    const SHF_TLS                   = 0x400
end
@constants SHF_MASKS "SHF_" begin
    const SHF_MASKOS                = 0x0f000000
    const SHF_MASKPROC              = 0xf0000000
end

# Section Group Flags
const GRP_COMDAT  = 0x1
const GRP_MASKOS  = 0x0ff00000
const GRP_MASKPROC    = 0xf0000000

# Symbol Binding
const STB_LOCAL   = 0
const STB_GLOBAL  = 1
const STB_WEAK    = 2
const STB_LOOS    = 10
const STB_HIOS    = 12
const STB_LOPROC  = 13
const STB_HIPROC  = 15

# Symbol Types
const STT_NOTYPE  = 0
const STT_OBJECT  = 1
const STT_FUNC    = 2
const STT_SECTION = 3
const STT_FILE    = 4
const STT_COMMON  = 5
const STT_TLS = 6
const STT_LOOS    = 10
const STT_HIOS    = 12
const STT_LOPROC  = 13
const STT_HIPROC  = 15

# Symbol Visibility
const STV_DEFAULT = 0
const STV_INTERNAL    = 1
const STV_HIDDEN  = 2
const STV_PROTECTED   = 3

# ptype constants
@constants P_TYPE "PT_" begin
    const PT_NULL           = 0
    const PT_LOAD           = 1
    const PT_DYNAMIC        = 2
    const PT_INTERP         = 3
    const PT_NOTE           = 4
    const PT_SHLIB          = 5
    const PT_PHDR           = 6
    const PT_TLS            = 7
    const PT_LOOS           = 0x60000000
    const PT_GNU_EH_FRAME   = 0x6474e550
    const PT_GNU_STACK      = 0x6474e551
    const PT_GNU_RELRO      = 0x6474e552
    const PT_PAX_FLAGS      = 0x65041580
    const PT_HIOS           = 0x6fffffff
    const PT_LOPROC         = 0x70000000
    const PT_HIPROC         = 0x7fffffff
end

# p_flags constants
const PF_X    = 0x1 #Execute
const PF_W    = 0x2 #Write
const PF_R    = 0x4 #Read
const PF_MASKOS   = 0x0ff00000  #Unspecified
const PF_MASKPROC = 0xf0000000  #Unspecified

# X86_64 relocations
@constants R_X86_64 "" begin
    const R_X86_64_NONE              =   0
    const R_X86_64_64                =   1
    const R_X86_64_PC32              =   2
    const R_X86_64_GOT32             =   3
    const R_X86_64_PLT32             =   4
    const R_X86_64_COPY              =   5
    const R_X86_64_GLOB_DAT          =   6
    const R_X86_64_JUMP_SLOT         =   7
    const R_X86_64_RELATIVE          =   8
    const R_X86_64_GOTPCREL          =   9
    const R_X86_64_32                =  10
    const R_X86_64_32S               =  11
    const R_X86_64_16                =  12
    const R_X86_64_PC16              =  13
    const R_X86_64_8                 =  14
    const R_X86_64_PC8               =  15
    const R_X86_64_DTPMOD64          =  16
    const R_X86_64_DTPOFF64          =  17
    const R_X86_64_TPOFF64           =  18
    const R_X86_64_TLSGD             =  19
    const R_X86_64_TLSLD             =  20
    const R_X86_64_DTPOFF32          =  21
    const R_X86_64_GOTTPOFF          =  22
    const R_X86_64_TPOFF32           =  23
    const R_X86_64_PC64              =  24
    const R_X86_64_GOTOFF64          =  25
    const R_X86_64_GOTPC32           =  26
    const R_X86_64_GOT64             =  27
    const R_X86_64_GOTPCREL64        =  28
    const R_X86_64_GOTPC64           =  29
    const R_X86_64_GOTPLT64          =  30
    const R_X86_64_PLTOFF64          =  31
    const R_X86_64_SIZE32            =  32
    const R_X86_64_SIZE64            =  33
    const R_X86_64_GOTPC32_TLSDESC   =  34
    const R_X86_64_TLSDESC_CALL      =  35
    const R_X86_64_TLSDESC           =  36
    const R_X86_64_IRELATIVE         =  37
end

@constants AUXV_TYPE "AT_" begin
    const AT_NULL         = 0
    const AT_IGNORE       = 1
    const AT_EXECFD       = 2
    const AT_PHDR         = 3
    const AT_PHENT        = 4
    const AT_PHNUM        = 5
    const AT_PAGESZ       = 6
    const AT_BASE         = 7
    const AT_FLAGS        = 8
    const AT_ENTRY        = 9
    const AT_NOTELF       = 10
    const AT_UID          = 11
    const AT_EUID         = 12
    const AT_GID          = 13
    const AT_EGID         = 14
    const AT_PLATFORM     = 15
    const AT_HWCAP        = 16
    const AT_CLKTCK       = 17
    const AT_SYSINFO      = 32
    const AT_SYSINFO_EHDR = 33
end

@constants DYNAMIC_TYPE "DT_" begin
    const DT_NULL         = 0
    const DT_NEEDED       = 1
    const DT_PLTRELSZ     = 2
    const DT_PLTGOT       = 3
    const DT_HASH         = 4
    const DT_STRTAB       = 5
    const DT_SYMTAB       = 6
    const DT_RELA         = 7
    const DT_RELASZ       = 8
    const DT_RELAENT      = 9
    const DT_STRSZ        = 10
    const DT_SYMENT       = 11
    const DT_INIT         = 12
    const DT_FINI         = 13
    const DT_SONAME       = 14
    const DT_RPATH        = 15
    const DT_SYMBOLIC     = 16
    const DT_REL          = 17
    const DT_RELSZ        = 18
    const DT_RELENT       = 19
    const DT_PLTREL       = 20
    const DT_DEBUG        = 21
    const DT_TEXTREL      = 22
    const DT_JMPREL       = 23
    const DT_BIND_NOW     = 24
    const DT_INIT_ARRAY	  = 25
    const DT_FINI_ARRAY	  = 26
    const DT_INIT_ARRAYSZ =	27
    const DT_FINI_ARRAYSZ =	28
    const DT_RUNPATH	  = 29
    const DT_FLAGS	      = 30
    const DT_ENCODING     = 32
    const OLD_DT_LOOS     = 0x60000000
    const DT_LOOS         = 0x6000000d
    const DT_HIOS         = 0x6ffff000
    const DT_VALRNGLO     = 0x6ffffd00
    const DT_VALRNGHI     = 0x6ffffdff
    const DT_ADDRRNGLO    = 0x6ffffe00
    const DT_CONFIG       = 0x6ffffefa
    const DT_DEPAUDIT     = 0x6ffffefb
    const DT_AUDIT        = 0x6ffffefc
    const DT_ADDRRNGHI    = 0x6ffffeff
    const DT_VERSYM       = 0x6ffffff0
    const DT_RELACOUNT    = 0x6ffffff9
    const DT_RELCOUNT     = 0x6ffffffa
    const DT_FLAGS_1      = 0x6ffffffb
    const DT_VERDEF       = 0x6ffffffc
    const DT_VERDEFNUM    = 0x6ffffffd
    const DT_VERNEED      = 0x6ffffffe
    const DT_VERNEEDNUM   = 0x6fffffff
    const OLD_DT_HIOS     = 0x6fffffff
    const DT_LOPROC       = 0x70000000
    const DT_AUXILIARY   = 0x7ffffffd
    const DT_FILTER       = 0x7fffffff
end


text:
module FieldsTests

using Test
using Gridap

using ..FieldsMocks
using ..FieldsMocks: GradMockField

f = MockField(2,Int)

g = ∇(f)

@test isa(g,GradMockField)

end # module


text:
using Test
using StringArrayEditor

@testset "StringArrayEditor" begin
    f = load("./data/testfile01.txt")
    @testset "Loading" begin
        @test length(f.data) == 9
    end
    l1 = Line(f,3)
    l2 = Line(f,5)
    @testset "Reference Lines" begin
        @test l1.ln == 3
        @test l1.file == f
        @test_throws AssertionError Line(f,11)
        @test_throws AssertionError Line(f,0)
        @test length(f.references) == 2
        StringArrayEditor.value(l1) == "Line 3"
        StringArrayEditor.value(l2) == "Line 5"
    end
    same = copy(l1)
    @testset "Copy lines" begin
        @test l1.ln == same.ln
        @test l1.file === same.file
        @test !(l1 === same)
    end
    before = Line(f,1)
    @testset "Delete Lines" begin
        delete!(l1)
        @test l1.file == nothing
        @test l1.ln == nothing
        @test StringArrayEditor.value(l2) == "Line 5"
        @test l2.ln == 4
        @test before.ln == 1
        @test StringArrayEditor.value(before) == "Line 1"
        @test same.file == nothing
        @test same.ln == nothing
        @test length(f.references) == 2
        @test length(f.data) == 8
    end
    after = Line(f,7)
    @testset "Insert Lines" begin
        li = insert!(l2,"Hello World")
        @test isa(li,Line)
        @test value(li) == "Hello World"
        @test l2.ln == 5
        @test f.data[l2.ln-1] == "Hello World"
        @test after.ln == 8
        @test StringArrayEditor.value(after) == "Line 7"
        insert!(l2,after)
        @test f.data[l2.ln-1] == "Line 7"
        ri = insert!(l2,["Vec 1","Vec 2", "Vec 3"])
        @test isa(ri,Range)
        @test value(ri) == ["Vec 1","Vec 2","Vec 3"]
        @test l2.ln == 9
        @test f.data[l2.ln-4:l2.ln-1] == ["Line 7", "Vec 1","Vec 2", "Vec 3"]
        @test after.ln == 12
    end
    f = load("./data/testfile01.txt")
    l1 = Line(f,2)
    l2 = Line(f,3)
    l3 = Line(f,7)
    lend = Line(f,length(f.data))
    @testset "Append Lines" begin
        a1 = append!(l2,"Hello World")
        @test value(a1) == "Hello World"
        @test a1.ln == 4
        @test l2.ln == 3
        @test StringArrayEditor.value(l2) == "Line 3"
        @test l1.ln == 2
        @test StringArrayEditor.value(l1) == "Line 2"
        @test l3.ln == 8
        @test StringArrayEditor.value(l3) == "Line 3"
        append!(lend,"asdf")
        @test f.data[end-1:end] == ["Line 8", "asdf"]
        nlend = Line(f,length(f.data))
        a2 = append!(nlend,["ffff","gggg"])
        @test value(a2) == ["ffff","gggg"]
        @test a2.from == nlend.ln+1
        @test f.data[end-3:end] == ["Line 8", "asdf", "ffff", "gggg"]
        append!(l2,["ffff","gggg"])
        @test l3.ln == 10
        @test StringArrayEditor.value(l3) == "Line 3"
        @test l2.ln == 3
        @test StringArrayEditor.value(l2) == "Line 3"
        @test f.data[l2.ln:l2.ln+2] == ["Line 3", "ffff", "gggg"]
    end
    @testset "Range Creation" begin
        f = load("./data/testfile01.txt")
        r1 = Range(f,2,5)
        @test r1.from == 2
        @test r1.to == 5
        @test_throws AssertionError Range(f,3,2)
        @test_throws AssertionError Range(f,2,100)
        @test length(r1) == 4
        @test collect(r1) == [2,3,4,5]
        @test length(f.references) == 1
        @test value(r1) == ["Line 2","Line 3","Line 4", "Line 5"]
        @test !(r1 === copy(r1))
    end
    @testset "Delete Range" begin
        f = load("./data/testfile01.txt")
        lbefore = Line(f,2)
        lafter = Line(f,8)
        rbefore = Range(f,1,2)
        rafter = Range(f,8,9)
        rdel = Range(f,3,6)
        lin1 = Line(f,4)
        lin2 = Line(f,6)
        rin = Range(f,4,5)
        roverlap = Range(f,5,7)
        rsurround = Range(f,2,9)
        @test length(f.references) == 10
        delete!(rdel)
        @test length(f.references) == 5
        @test lbefore.ln == 2
        @test length(rbefore) == 2
        @test rbefore.from == 1
        @test rbefore.to == 2
        @test lafter.ln == 4
        @test length(rafter) == 2
        @test rafter.from == 4
        @test rafter.to == 5
        @test length(rsurround) == 4
        @test rsurround.from == 2
        @test rsurround.to == 5
        @test StringArrayEditor.isdestroyed(roverlap)
        @test StringArrayEditor.isdestroyed(lin1)
        @test StringArrayEditor.isdestroyed(lin2)
        @test StringArrayEditor.isdestroyed(rdel)
        @test StringArrayEditor.isdestroyed(rin)
        @test f.data == [ "Line 1",
                          "Line 2",
                          "Line 3",
                          "Line 7",
                          "Line 8"]
        delete!(lafter)
        @test f.data == [ "Line 1",
                          "Line 2",
                          "Line 3",
                          "Line 8"]
        @test lbefore.ln == 2
        @test rbefore.from == 1
        @test rbefore.to == 2
        @test StringArrayEditor.isdestroyed(lafter)
        @test length(rsurround) == 3
        @test length(f.references) == 4
    end
    @testset "Insert Range" begin
        f = load("./data/testfile01.txt")
        lbefore = Line(f,2)
        lafter = Line(f,8)
        rbefore = Range(f,1,2)
        rafter = Range(f,8,9)
        rins = Range(f,3,6)
        lin1 = Line(f,4)
        lin2 = Line(f,6)
        rin = Range(f,4,5)
        roverlap = Range(f,5,7)
        rsurround = Range(f,2,9)
        i2 = insert!(rbefore,"asdf")
        @test value(i2) == "asdf"
        @test i2.ln == rbefore.from-1
        @test lafter.ln == 9
        @test rafter.from == 9
        i3 = insert!(rin,["asdf","ffff"])
        @test value(i3) == ["asdf","ffff"]
        @test rin.from == 7
        @test rin.to == 8
        @test rsurround.from == 3
        @test rsurround.to == 12
        @test lin2.ln == 9
        insert!(lin1,"fff")
        @test rafter.from == 12
        insert!(lin1,rbefore)
        @test f.data[8] == "Line 1"
        @test f.data[9] == "Line 2"
        insert!(lin1,lin1)
        @test f.data[10] == "Line 4"
    end
@testset "Append Range" begin
    f = load("./data/testfile01.txt")
    lbefore = Line(f,2)
    lafter = Line(f,8)
    rbefore = Range(f,1,2)
    rafter = Range(f,8,9)
    rins = Range(f,3,6)
    lin1 = Line(f,4)
    lin2 = Line(f,6)
    rin = Range(f,4,5)
    roverlap = Range(f,5,7)
    rsurround = Range(f,2,9)
    ra = append!(rbefore,"asdf")
    @test value(ra) == "asdf"
    @test rbefore.from == 1
    @test rbefore.to == 2
    @test lafter.ln == 9
    @test rafter.from == 9
    ra2 = append!(rin,["asdf","ffff"])
    @test value(ra2) == ["asdf","ffff"]
    @test rin.from == 5
    @test rin.to == 6
    @test rsurround.from == 3
    @test rsurround.to == 12
    @test lin2.ln == 9
    append!(lin1,"fff")
    @test rafter.from == 12
    lin1
    append!(lin1,rbefore)
    @test f.data[6] == "Line 1"
    @test f.data[7] == "Line 2"
    append!(lin1,lin1)
    @test f.data[6] == "Line 4"
end
@testset "Replace" begin

    f = load("./data/testfile01.txt")
    lbefore = Line(f,2)
    lafter = Line(f,8)
    rbefore = Range(f,1,2)
    rafter = Range(f,8,9)
    rins = Range(f,3,6)
    lin1 = Line(f,4)
    lin2 = Line(f,6)
    rin = Range(f,4,5)
    roverlap = Range(f,5,7)
    rsurround = Range(f,2,9)
    replace!(lbefore,"ASDF")
    @test value(lbefore) == "ASDF"

    rnew = replace!(lbefore,["GGG","AAA"])
    @test f.data[2:3] == ["GGG","AAA"]
    @test lafter.ln == 9
    @test rbefore.from == 1
    @test rbefore.to == 3
    @test rafter.from == 9
    @test rafter.to == 10
    @test rins.from == 4
    @test rins.to == 7
    @test lin1.ln == 5
    @test lin2.ln == 7
    @test f.data[rsurround.from] == "GGG"

    f = load("./data/testfile01.txt")
    lbefore = Line(f,2)
    lafter = Line(f,8)
    rbefore = Range(f,1,2)
    rafter = Range(f,8,9)
    rins = Range(f,3,6)
    lin1 = Line(f,3)
    lin2 = Line(f,6)
    rin = Range(f,4,5)
    roverlap = Range(f,5,7)
    rsurround = Range(f,2,9)
    f.data
    replace!(lbefore,"ASDF")
    length(f.references)
    lnew = replace!(rins,"GGG")
    @test lnew.ln == 3
    @test StringArrayEditor.isdestroyed(lin1)
    @test StringArrayEditor.isdestroyed(lin2)
    @test StringArrayEditor.isdestroyed(roverlap)
    @test StringArrayEditor.isdestroyed(rin)
    @test lafter.ln == 8-3
    @test rsurround.to == 9-3
    @test rsurround.from == 2
    @test length(f.references) == 6
    @test f.data == [
        "Line 1",
        "ASDF"  ,
        "GGG"   ,
        "Line 3",
        "Line 7",
        "Line 8"
    ]
    l = replace!(lbefore,lafter)
    @test value(l) == "Line 7"
    l = replace!(lbefore,rafter)
    @test StringArrayEditor.isdestroyed(lbefore)
    @test value(l) == ["Line 7","Line 8"]
end

@testset "Search Functions" begin

    f = load("./data/testfile01.txt")
    l = Line(f,r"Line 2")
    @test l.ln == 2
    l2 = Line(f,r"Line 2",after=r"Line 5")
    @test l2.ln == 6
    l3 = Line(f,r"Line 2",after=l)
    @test l3.ln == 6
    @test_throws ErrorException Line(f,r"Line 7",before=l2)
    r = Range(f,from=r"Line 2",to=r"Line 2")
    @test r.from == 2
    @test r.to == 6
    r = Range(f,after=l,to=r"Line 3")
    @test r.from == 1
    @test r.to == 7
    r = Range(f,from=l,after=l,to=r"Line 3")
    @test r.from == 2
    @test r.to == 7
    replace!(l3,"ASDF")
    r = Range(f,from=l,until=r"Line")
    @test r.from == 2
    @test r.to == 5
end

@testset "Save" begin
    f = load("./data/testfile01.txt")
    save(f,"test.txt")
    @test readlines("test.txt") == f.data
    rm("test.txt")
end
@testset "Line and Range manipulations" begin
    f = load("./data/testfile01.txt")
    l1 = Line(f,1)
    l2 = l1+3
    @test l2.ln == 4
    l3 = l2-1
    @test l3.ln == 3
    r1 = Range(f,2,5)
    l2 = r1-1
    @test l2.ln == 1
    l3 = r1+2
    @test l3.ln == 7
    r2 = r1[2:end-1]
    @test r2.from == 3
    @test r2.to == 4
    @test_throws BoundsError r2[-1:end]
    @test_throws BoundsError r2[0:end]
    @test_throws BoundsError r2[1:20]
end
@testset "Move" begin
    f = load("./data/testfile01.txt")
    l = Line(f,1)
    r = replace!(l,["foo","bar","baz"])
    l = replace!(Line(f,8),"ASDF")
    move!(r,l)
    @test l.ln == 8
    f.data[l.ln-3:l.ln-1] == ["foo","bar","baz"]
    @test r.from == l.ln-3
    @test r.to == l.ln-1
    @test value(r) == ["foo","bar","baz"]
    move!(l,r)
    @test l.ln == 5
    @test value(l) == "ASDF"

end
@testset "Range by lines" begin
    f = load("./data/testfile01.txt")
    r = Range(f,from=Line(f,1),to=Line(f,4))
    @test value(r) == f.data[1:4]
end
end


text:
using Unitful, Dates

"""
    AbstractMeshArray{T, N}

Subtype of AbstractArray{T, N}
"""
abstract type AbstractMeshArray{T, N} <: AbstractArray{T, N} end

"""
    gcmgrid

gcmgrid data structure. Available constructors:

```
gcmgrid(path::String, class::String, nFaces::Int,
        fSize::Array{NTuple{2, Int},1},
        ioSize::Union{NTuple{2, Int},Array{Int64,2}},
        ioPrec::Type, read::Function, write::Function)
```

The `class` can be set to "LatLonCap", "CubeSphere", "PeriodicChannel", "PeriodicDomain".

For example, A periodic channel (periodic in the x direction) of size 360 by 160, can be defined as follows.

```
pth=MeshArrays.GRID_LL360
class="PeriodicChannel"
ioSize=(360, 160)
ioPrec=Float32

γ=gcmgrid(pth, class, 1, [ioSize], ioSize, ioPrec, read, write)

Γ=GridLoad(γ)
```

Please refer to `GridSpec` and `UnitGrid` for more info related to `gcmgrid` options.
"""
struct gcmgrid
  path::String
  class::String
  nFaces::Int
  fSize::Array{NTuple{2, Int},1}
  ioSize::Union{NTuple{2, Int},Array{Int64,2}}
  ioPrec::Type
  read::Function
  write::Function
end

"""
    varmeta

varmeta data structure. By default, `unit` is `missing` (non-dimensional), `position`
is `fill(0.5,3)` (cell center), time is `missing`, and `name` / `long_name` is `unknown`.

Available constructors:

```
varmeta(unit::Union{Unitful.Units,Number,Missing},position::Array{Float64,1},
        time::Union{DateTime,Missing,Array{DateTime,1}},name::String,long_name::String)
```

And:

```defaultmeta = varmeta(missing,fill(0.5,3),missing,"unknown","unknown")```

"""
struct varmeta
  unit::Union{Unitful.Units,Number,Missing}
  position::Array{Float64,1}
  time::Union{DateTime,Missing,Array{DateTime,1}}
  name::String
  long_name::String
end

defaultmeta = varmeta(missing,fill(0.5,3),missing,"unknown","unknown")

## concrete types and MeshArray alias:

include("Type_gcmfaces.jl");
OuterArray{T,N}=Array{T,N} where {T,N}
InnerArray{T,N}=Array{T,N} where {T,N}
include("Type_gcmarray.jl");
include("Type_gcmvector.jl");

"""
    MeshArray

Alias to `gcmarray` or `gcmfaces` concrete type
"""
MeshArray=gcmarray

## Methods that apply to all AbstractMeshArray types

import Base: maximum, minimum, sum, fill, fill!

function maximum(a::AbstractMeshArray)
  c=-Inf;
  for I in eachindex(a)
    c = max(c,maximum(a[I]))
  end
  return c
end

function minimum(a::AbstractMeshArray)
  c=Inf;
  for I in eachindex(a)
    c = min(c,minimum(a[I]))
  end
  return c
end

function sum(a::AbstractMeshArray)
  c=0.0
  for I in eachindex(a)
    c = c + sum(a[I])
  end
  return c
end

function fill(val::Any,a::AbstractMeshArray)
  c=similar(a)
  for I in eachindex(a.f)
    c.f[I] = fill(val,size(a.f[I]))
  end
  return c
end

function fill!(a::AbstractMeshArray,val::Any)
  for I in eachindex(a.f)
    fill!(a.f[I],val)
  end
  return a
end

import Base: +, -, *, /

function +(a::AbstractMeshArray,b::AbstractMeshArray)
  c=similar(a)
  for I in eachindex(a.f)
    c.f[I] = a.f[I] + b.f[I]
  end
  return c
end

function -(a::AbstractMeshArray,b::AbstractMeshArray)
  c=similar(a)
  for I in eachindex(a.f)
    c.f[I] = a.f[I] - b.f[I]
  end
  return c
end

function /(a::AbstractMeshArray,b::AbstractMeshArray)
  c=similar(a)
  for I in eachindex(a.f)
    c.f[I] = a.f[I] ./ b.f[I]
  end
  return c
end

function *(a::AbstractMeshArray,b::AbstractMeshArray)
  c=similar(a)
  for I in eachindex(a.f)
    c.f[I] = a.f[I] .* b.f[I]
  end
  return c
end


text:
s = open("input.txt") do file
    readlines(file);
end

nums = map(x-> parse(Int, x), s);

preambleLength = 25;

function ValidNum(nums, i)
    lookingFor = Set();
    for j = i-preambleLength:i-1
        if in(nums[j], lookingFor)
            return true;
        end

        if nums[i] - nums[j] > 0
            push!(lookingFor, nums[i] - nums[j])
        end
    end

    return false;
end

for i = preambleLength+1:length(nums)
    if !ValidNum(nums, i)
        println(nums[i])
    end
end


text:
module HoboSort
include("Extract.jl")
include("TMixtures2.jl")
include("HoboPlot.jl")
include("cholutils.jl")
using Base.Dates, Base.LinAlg, Distributions, PDMats, .Extract, .Mixtures, .HoboPlot

immutable SortParams
    init_period::Second
    inspect_period::Second
    history_period::Second
    min_separation_periods::Int
    split_score_improvement::Float64
    merge_score_improvement::Float64
    min_cluster_size::Int
    initial_clusters::Int
    rerun_em::Bool
    verbose::Bool
    plot::Bool
    movie_file::UTF8String
    graph_file::UTF8String
end
SortParams(; init_period::TimePeriod=Minute(10), inspect_period::TimePeriod=Minute(5),
           history_period::TimePeriod=Minute(10), min_separation_periods::Int=6,
           merge_score_improvement::Real=0, split_score_improvement::Real=0,
           min_cluster_size::Integer=120, initial_clusters::Integer=10, rerun_em::Bool=true,
           verbose::Bool=false, plot::Bool=false, movie_file::AbstractString="",
           graph_file::AbstractString="") =
	SortParams(init_period, inspect_period, history_period, min_separation_periods,
               split_score_improvement, merge_score_improvement,
               min_cluster_size, initial_clusters, rerun_em, verbose, plot, movie_file, graph_file)

function fit_tdist{T}(X::AbstractMatrix{T}, ν;
                      maxiter::Int=100, tol=1e-3, μ_init::Vector{T}=vec(mean(X, 2)),
                      Σ_init::Union{PDMat{T},Void}=nothing)
    μ = copy(μ_init)
    Xmμ = X .- μ
    Σ::PDMat{T,Matrix{T}} = Σ_init === nothing ? PDMat(LinAlg.copytri!(Base.LinAlg.BLAS.syrk('U', 'N', 1/size(X, 2), Xmμ), 'U')) : deepcopy(Σ_init)
    mahal = η = zeros(eltype(X), size(X, 2))
    dist = MvTDist(ν, μ, Σ)
    ll = -Inf
    for t = 1:maxiter
        oldll = ll
        ll = 0.0

        # E step
        broadcast!(-, Xmμ, X, μ)
        mahal = invquad!(mahal, Σ, Xmμ)
        shdfhdim, v = Distributions.mvtdist_consts(dist)
        for i = 1:size(X, 2)
            ll += v - shdfhdim * log1p(mahal[i] / ν)
            η[i] = (ν + size(X, 1))/(ν + mahal[i])
        end

        # M step
        Mixtures.weighted_mean_cov!(μ, Σ.mat, Xmμ, X, η)
        copy!(Σ.chol.factors, Σ.mat)
        LinAlg.chol!(Σ.chol.factors, Val{:U})

        abs(ll - oldll) < tol && return (dist, η, ll)
    end
    warning("T distribution fit did not converge")
    (dist, η, ll)
end

type Cluster{T,S<:AbstractMatrix}
    active::Bool
    first_appeared::Int
    split_from::Int
    separation::Vector{T}
    wf::S
    η::Vector{T}
    sumη::T
    sum::Vector{T}
    comoments::Matrix{T}
    comoment_chol::Cholesky{T,Matrix{T}}
    indexes::Vector{Int}
    dist::MvTDist
    scratch::Vector{T}
end

function Cluster{T}(ν::T, wf::AbstractMatrix{T}, indexes::Vector{Int}, first_appeared::Int, split_from::Int=-1)
	@assert issorted(indexes)
    n = length(indexes)
    cluster_wf = wf[:, indexes]

    dist, η = fit_tdist(cluster_wf, ν)
    sumη = sum(η)

    Cluster(true, first_appeared, split_from, T[], wf, η, sumη, dist.μ*sumη, dist.Σ.mat*sumη,
            Cholesky{T,Matrix{T}}(dist.Σ.chol.factors*sqrt(sumη), 'U'),
            indexes, dist, zeros(T, size(wf, 1)))
end

"""
    refit!(cluster::Cluster)

Refit distribution based on current parameters, returning log likelihood
"""
function refit!{T}(cluster::Cluster{T})
    cluster_wf = cluster.wf[:, cluster.indexes]
    dist, η, ll = fit_tdist(cluster_wf, cluster.dist.df; μ_init=cluster.dist.μ, Σ_init=cluster.dist.Σ)
    sumη = sum(η)
    cluster.η = η
    cluster.sumη = sumη
    cluster.sum = dist.μ*sumη
    cluster.comoments = dist.Σ.mat*sumη
    cluster.comoment_chol = Cholesky{T,Matrix{T}}(dist.Σ.chol.factors*sqrt(sumη), 'U')
    cluster.dist = dist
    ll
end

"""
    updatedist!(cluster::Cluster)

Update distribution parameters based on our forms, which are more convenient for online updates
"""
function updatedist!(cluster)
    n = cluster.sumη
    broadcast!(*, cluster.dist.μ, cluster.sum, 1/n)
    broadcast!(*, cluster.dist.Σ.mat, cluster.comoments, 1/n)
    broadcast!(*, cluster.dist.Σ.chol.factors, cluster.comoment_chol.factors, 1/sqrt(n))
    # mu = sum(cluster.wf[:, cluster.indexes].*cluster.η', 2)/cluster.sumη
    # demeaned = (cluster.wf[:, cluster.indexes] .- mu)
    # @assert isapprox(cluster.dist.μ, mu)
    # @assert isapprox(cluster.dist.Σ.mat, (demeaned.*cluster.η')*demeaned'/cluster.sumη)
    # @assert isapprox(cluster.dist.Σ.chol[:U], cholfact(cluster.dist.Σ.mat)[:U])

    return cluster
end

"""
    addspike!(cluster::Cluster, spike::AbstractVector)

Add spike to cluster
"""
function addspike!{T}(cluster::Cluster{T}, spike::AbstractVector, index::Int)
    scratch = cluster.scratch
    broadcast!(-, scratch, spike, cluster.dist.μ)
    mahal = invquad(cluster.dist.Σ, scratch)
    shdfhdim, v = Distributions.mvtdist_consts(cluster.dist)
    η = (cluster.dist.df + size(cluster.wf, 1))/(cluster.dist.df + mahal)

    oldsumη = cluster.sumη
    sumη = cluster.sumη += η
    sum_ = cluster.sum
    push!(cluster.indexes, index)
    push!(cluster.η, η)

    # Update sum and mean
    factor = sqrt(η*oldsumη/sumη)
    for i = 1:size(spike, 1)
        scratch[i] *= factor
        sum_[i] += η*spike[i]
    end

    # Update comoments
    LinAlg.copytri!(BLAS.syr!('U', one(T), scratch, cluster.comoments), 'U')
    lowrankupdate!(cluster.comoment_chol, scratch)

    # Update everything else
    updatedist!(cluster)
end

"""
    prune!(cluster::Cluster, wf::AbstractMatrix{T}, before::Real)

Remove old spikes from cluster
"""
function prune!{T}(cluster::Cluster{T}, before::Int)
    wf = cluster.wf
    sum_ = cluster.sum
    scratch = cluster.scratch
    indexes = cluster.indexes
    η = cluster.η
    for iprune = 1:searchsortedlast(indexes, before)
        n = length(indexes)

        # Make sure cluster is still full rank
        n > size(wf, 1)+1 || (cluster.active = false; return cluster)

        spike = sub(wf, :, shift!(indexes))
        oldsumη = cluster.sumη
        spikeη = shift!(η)
        sumη = cluster.sumη -= spikeη
        factor = sqrt(spikeη*sumη/oldsumη)

        # Update sum and mean
        for i = 1:size(wf, 1)
            sum_[i] -= spikeη*spike[i]
            scratch[i] = (spike[i] - sum_[i]/sumη)*factor
        end

        # Update comoments
        LinAlg.copytri!(BLAS.syr!('U', -one(T), scratch, cluster.comoments), 'U')
        lowrankdowndate!(cluster.comoment_chol, scratch)
    end
    @assert minimum(indexes) > before

    # Update everything else
    updatedist!(cluster)
end

function active_cluster_indexes(clusters)
    active_clusters = Int[]
    for icluster = 1:length(clusters)
        if clusters[icluster].active
            push!(active_clusters, icluster)
        end
    end
    active_clusters
end

function mixture_model_from_clusters(clusters, X, λ)
    dists = MvTDist[cluster.dist for cluster in clusters]
    logπ = Float64[length(cluster.indexes) for cluster in clusters]
    nspikes = sum(logπ)
    for i = 1:length(clusters)
        logπ[i] = log(logπ[i]/nspikes)
    end
    TMixture(dists, X; logπ=logπ, λ=λ, fix_ν=true, use_uniform=false)
end

function cluster_separation(clusters, X, λ)
    m = mixture_model_from_clusters(clusters, X, λ)
    γ = Mixtures.computeγ(m)
    s = zeros(length(clusters), length(clusters))
    n = zeros(Int, length(clusters), length(clusters))

    # This is only approximately right, since it's not necessarily
    # _exactly_ the cluster the spike is in, but should be close enough
    ass = zeros(Int, size(X, 2))
    maxes = fill(-Inf, size(X, 2))
    for k = 1:length(clusters)
        for i = 1:length(ass)
            if γ[i, k] > maxes[i]
                ass[i] = k
                maxes[i] = γ[i, k]
            end
        end
    end

    for j = 1:length(clusters), i = 1:length(ass)
        icluster = ass[i]
        s[icluster, j] += γ[i, j]
        n[icluster, j] += 1
    end
    s./n
end

function check_when_split_from(clusters, cluster, split_from)
    while cluster.split_from != -1
        @show cluster.split_from
        if cluster.split_from == split_from
            return cluster.first_appeared
        end
        cluster = clusters[cluster.split_from]
    end
    return -1
end

immutable SplitAction
    oldcluster::Int
    newcluster::Int
    iinspect::Int
end
immutable MergeAction
    deadcluster::Int
    targetcluster::Int
    iinspect::Int
end

using PyPlot
function spikesort{T}(wf::AbstractMatrix{T}, times::Vector, params::SortParams)
    λ = 1e-6
    size(wf, 2) == length(times) || throw(ArgumentError("second dimension of wf must match times"))
    verbose = params.verbose
    local assignments_::Vector{Int}

    pca_in = similar(wf)
    rg = 1:searchsortedlast(times, first(times)+Int(params.history_period))
    s = sum(sub(wf, :, rg), 2)
    pca_in[:, rg] = sub(wf, :, rg) .- s/length(rg)
    v = sumabs2(sub(pca_in, :, rg))
    scale!(sub(pca_in, :, rg), length(rg)/v)
    cur_first = 1
    n = length(rg)
    for j = last(rg)+1:size(wf, 2)
        last_index = searchsortedlast(times, times[j]-Int(params.history_period))
        while cur_first < last_index
            for k = 1:size(wf, 1)
                c = wf[k, cur_first]
                olds = s[k]
                news = s[k] -= c
                v -= (c - olds/n)*(c - news/(n-1))
            end
            n -= 1
            cur_first += 1
        end
        for k = 1:size(wf, 1)
            c = wf[k, j]
            olds = s[k]
            news = s[k] += c
            v += (c - olds/n)*(c - news/(n+1))
        end
        for k = 1:size(wf, 1)
            pca_in[k, j] = ((n+1)*wf[k, j] - s[k])/v
        end
        # @assert isapprox(s, sum(sub(wf, :, last_index:j), 2))
        # @assert isapprox(v, sumabs2(sub(wf, :, last_index:j) .- mean(sub(wf, :, last_index:j), 2)))
        n += 1
    end
    C = Base.covzm(pca_in, vardim=2)
    vals, vecs = eig(Symmetric(C))
    reverse!(vals)
    vecs = flipdim(vecs, 2)
    nkeep = searchsortedfirst(cumsum(vals./sum(vals)), 0.9)
    verbose && println("Using $nkeep PCs")
    keepvecs = vecs[:, 1:nkeep]
    cluster_space::Matrix{T} = keepvecs'*pca_in
    # Drop spikes that are more than 10 SDs away from the mean along any dimension
    # Chebyshev's inequality guarantees this drops no more than 1% of spikes, but
    # in practice it will be far less.
    sd = std(cluster_space, 2)
    mu = mean(cluster_space, 2)
    dropped = Set()
    for j = 1:size(cluster_space, 2), i = 1:size(cluster_space, 1)
        if abs(cluster_space[i, j] - mu[i]) > 10*sd[i]
            cluster_space[i, j] = mu[i]
            push!(dropped, j)
        end
    end
    verbose && println("Dropping $(length(dropped)) spikes >10 SD away")
    broadcast!(-, cluster_space, cluster_space, minimum(cluster_space, 2))
    broadcast!(/, cluster_space, cluster_space, maximum(cluster_space, 2))

    # Initialize using free split-merge GMM
    verbose && println("Initializing Clusters...")
    init_rg = 1:searchsortedlast(times, first(times)+Int(params.init_period))
    init_wft = cluster_space[:, init_rg]
    m = TMixture(params.initial_clusters, init_wft; λ=λ, use_uniform=true)::TMixture{T,Matrix{T}}
    m = fsmem!(m, maxiter=500, verbose=:op)::TMixture{T,Matrix{T}}

    # Drop clusters that are too small
    assignments_ = assignments(m)
    init_nclusters = length(m.components)
    ν = m.components[1].df
    verbose && println("Using $init_nclusters clusters")
    counts = hist(assignments_, 0:init_nclusters)[2]
    verbose && println("Cluster counts: $counts")
    valid = counts .>= params.min_cluster_size
    if !all(valid)
        remap = zeros(Int, init_nclusters+1)
        remap[find(valid)+1] = 1:sum(valid)
        assignments_ = [remap[x+1] for x in assignments_]
        init_nclusters = sum(valid)
        verbose && println("Dropping $(sum(!valid)) clusters with fewer than $(params.min_cluster_size) spikes")
    end
    counts = hist(assignments_, 0:init_nclusters)[2]
    verbose && println("Cluster counts: $counts")

    # Compute sum and comoment for each cluster
    clusters = Array(Cluster{T,typeof(cluster_space)}, init_nclusters)
    for icluster = 1:length(clusters)
        indexes = find(assignments_ .== icluster)
        clusters[icluster] = Cluster(ν, cluster_space, indexes, 1)
    end

    noise_indexes = find(assignments_ .== 0)+1
    resize!(assignments_, size(wf, 2))

    if params.plot
        hp = HoboPlotter(params.movie_file)
        plotclusters(hp, wf[:, init_rg], sub(assignments_, init_rg), Int(params.history_period))
    end

    connections = Tuple{Int,Int}[]
    connections_time = Tuple{ASCIIString,ASCIIString}[]

    iinspect = 0
    next_inspect = zero(eltype(times))
    prune_index = 0
    for ispike = last(init_rg)+1:size(wf, 2)
        if ispike in dropped
            assignments_[ispike] = 0
            continue
        end

        # TODO fix Distributions so this can use SubArrays
        cur_wf = cluster_space[:, ispike]
        cur_time = times[ispike]
        last_prune_index = prune_index
        prune_index = searchsortedlast(times, cur_time - Int(params.history_period))

        # See if we need to split/merge clusters
        if cur_time >= next_inspect
            println("Inspecting...")
            iinspect += 1
            rg = last_prune_index+1:ispike-1
            sep_rg = rg[assignments_[rg] .!= 0]
            sep_wft = cluster_space[:, sep_rg]

            if params.rerun_em
                active_clusters = active_cluster_indexes(clusters)
                m = mixture_model_from_clusters(clusters[active_clusters], sep_wft, λ)
                em!(m)
                for icluster = 1:length(active_clusters)
                    cluster = clusters[active_clusters[icluster]]
                    indexes = sep_rg[assignments(m) .== icluster]
                    assignments_[indexes] = active_clusters[icluster]
                    cluster.indexes = indexes
                    if length(indexes) < 2*size(cluster_space, 1)
                        println("Cluster $(active_clusters[icluster]) deactivated by EM")
                        cluster.active = false
                    else
                        refit!(cluster)
                    end
                end
            else
                for icluster = 1:length(clusters)
                    !clusters[icluster].active && continue
                    refit!(clusters[icluster])
                end
            end

            changed = false
            more_iter = true
            while more_iter
                more_iter = false

                # Try merging clusters
                for iclust1 = 1:length(clusters)
                    clust1 = clusters[iclust1]
                    !clust1.active && continue
                    for iclust2 = iclust1+1:length(clusters)
                        clust2 = clusters[iclust2]
                        !clust2.active && continue
                        indexes = [clust1.indexes; clust2.indexes]

                        # mix = Mixture([MvTDist(10., clust1.μ, deepcopy(clust1.Σ)),
                        #                MvTDist(10., clust2.μ, deepcopy(clust2.Σ))])
                        # γ_1 = zeros(length(rgindexes))
                        # γ_2 = zeros(length(rgindexes))
                        # γ_1[1:length(clust1.indexes)] = 1
                        # γ_2[length(clust1.indexes)+1:end] = 1
                        # mtest = em!(SoftMixtureFit(mix, wft; λ=0.01, γ=Vector{Float64}[γ_1, γ_2]))
                        # mtest = em!(HardMixtureFit(mix, wft; λ=λ,
                        #                       assignments=[fill(UInt8(1), length(clust1.indexes)); fill(UInt8(2), length(clust2.indexes))]))
                        mtest = em!(TMixture(typeof(clust1.dist)[deepcopy(clust1.dist), deepcopy(clust2.dist)], cluster_space[:, indexes];
                                    λ=λ, logπ=[log(length(clust1.indexes)/length(indexes)), log(length(clust2.indexes)/length(indexes))],
                                    fix_ν=true, use_uniform=false))
                        mmerged = merge(mtest)

                        if mmerged !== mtest && (score(mmerged) - score(mtest)) > params.merge_score_improvement
                            changed = more_iter = true
                            newclust = Cluster(ν, cluster_space, sort!([clust1.indexes; clust2.indexes]), iinspect)
                            sep1 = maximum(keepvecs'*(abs(clust1.dist.μ) - sqrt(diag(clust1.dist.Σ))))
                            sep2 = maximum(keepvecs'*(abs(clust2.dist.μ) - sqrt(diag(clust2.dist.Σ))))
                            if sep2 < sep1
                                clusters[iclust1] = newclust
                                clust2.active = false
                                inewclust = iclust1
                                ideadclust = iclust2
                            else
                                clusters[iclust2] = newclust
                                clust1.active = false
                                inewclust = iclust2
                                ideadclust = iclust1
                            end
                            assignments_[clusters[ideadclust].indexes] = inewclust
                            if verbose
                                println("Merging cluster $ideadclust into $inewclust")
                                println("    cluster $iclust1 (split from $(clust1.split_from)) had $(length(clust1.indexes)) waveforms with mean separation $(sep1)")
                                println("    cluster $iclust2 (split from $(clust2.split_from)) had $(length(clust2.indexes)) waveforms with mean separation $(sep2)")
                                println("    score $(score(mtest)) -> $(score(mmerged))")
                            end
                            push!(connections, (ideadclust,inewclust))
                            clust1_split_from_clust2 = check_when_split_from(clusters, clust1, iclust2)
                            clust2_split_from_clust1 = check_when_split_from(clusters, clust2, iclust1)
                            if clust1_split_from_clust2 != -1 || clust2_split_from_clust1 != -1
                                if clust1_split_from_clust2 != -1
                                    split_periods = iinspect - clust1_split_from_clust2
                                    ioriginal = iclust2
                                    isplit = iclust1
                                else
                                    split_periods = iinspect - clust2_split_from_clust1
                                    ioriginal = iclust1
                                    isplit = iclust2
                                end
                                println("    split for $split_periods periods")
                                if split_periods < params.min_separation_periods
                                    println("    min number of separation periods is $(params.min_separation_periods), so re-merging")
                                    assignments_[assignments_ .== isplit] = ioriginal
                                end
                            end
                            break;
                        end
                    end
                end

                # Try splitting clusters
                for icluster = 1:length(clusters)
                    oldclust = clusters[icluster]
                    indexes = oldclust.indexes
                    !oldclust.active && continue
                    length(indexes) <= params.min_cluster_size && continue

                    mtest = em!(TMixture(typeof(oldclust.dist)[deepcopy(oldclust.dist)], cluster_space[:, indexes]; λ=λ, fix_ν=true, use_uniform=false))
                    msplit = split(mtest)

                    if msplit !== mtest && (score(msplit) - score(mtest)) > params.split_score_improvement
                        assign = assignments(msplit)
                        assign1 = assign .== 1
                        assign2 = assign .== 2
                        cluster1_large_enough = sum(assign1) >= params.min_cluster_size
                        cluster2_large_enough = sum(assign2) >= params.min_cluster_size
                        changed = more_iter = true
                        if cluster1_large_enough && cluster2_large_enough

                            # Compute separation for new clusters
                            active_clusters = active_cluster_indexes(clusters)
                            deleteat!(active_clusters, findfirst(active_clusters, icluster))
                            splitclust = Cluster(ν, cluster_space, indexes[assign1], iinspect)
                            newclust = Cluster(ν, cluster_space, indexes[assign2], iinspect)
                            seps = cluster_separation([splitclust; newclust; clusters[active_clusters]], sep_wft, λ)

                            splitsep = -sum(seps[1, 3:end])
                            newsep = -sum(seps[2, 3:end])
                            if splitsep > newsep
                                # The new cluster should be the better separated of the two
                                splitclust, newclust = newclust, splitclust
                                splitsep, newsep = newsep, splitsep
                            end
                            splitclust.split_from = oldclust.split_from
                            newclust.split_from = icluster

                            for i = 1:length(indexes)
                            end

                            clusters[icluster] = splitclust
                            push!(clusters, newclust)
                            assignments_[newclust.indexes] = length(clusters)
                            push!(connections, (icluster,length(clusters)))
                            if verbose
                                println("Splitting cluster $icluster")
                                println("    new cluster $(length(clusters)) has $(length(newclust.indexes)) waveforms with mean separation $(newsep)")
                                println("    old cluster $(icluster) has $(length(splitclust.indexes)) waveforms with mean separation $(splitsep)")
                                println("    score $(score(mtest)) -> $(score(msplit))")
                            end
                        else
                            if !cluster1_large_enough && !cluster2_large_enough
                                # Deactivate original cluster
                                noise_idx = oldclust.indexes
                                oldclust.active = false
                                if verbose
                                    println("Cluster $icluster is small but beneficial to split further; assigning to noise")
                                end
                            else
                                if cluster1_large_enough
                                    oldclust.indexes = indexes[assign1]
                                    noise_idx = indexes[assign2]
                                    append!(noise_indexes, noise_idx)
                                elseif cluster2_large_enough
                                    oldclust.indexes = indexes[assign2]
                                    noise_idx = indexes[assign1]
                                    append!(noise_indexes, noise_idx)
                                end
                                if verbose
                                    println("Assigning $(length(noise_idx)) spikes from $icluster to noise")
                                end
                            end
                            assignments_[noise_idx] = 0
                            append!(noise_indexes, noise_idx)
                            sort!(noise_indexes)
                        end
                    end
                end
            end

            # Add separation to the cluster info
            active_clusters = active_cluster_indexes(clusters)
            seps = cluster_separation(clusters[active_clusters], sep_wft, λ)
            for i = 1:length(active_clusters)
                push!(clusters[active_clusters[i]].separation, seps[i, i])
            end

            if params.plot
                plotclusters(hp, wf[:, rg], sub(assignments_, rg), Int(params.history_period))
            end

            next_inspect = cur_time+Int(params.inspect_period)
        end

        # Find best cluster
        deleteat!(noise_indexes, 1:searchsortedlast(noise_indexes, prune_index))
        best_cluster = 0
        nspikes = ispike - last_prune_index - 1
        best_logpdf = log((length(noise_indexes)+1)/(nspikes+1))
        for icluster = 1:length(clusters)
            cluster = clusters[icluster]
            !cluster.active && continue
            logpdf_ = logpdf(cluster.dist, cur_wf) + log(length(cluster.indexes)/(nspikes+1))
            prune!(cluster, prune_index)
            if logpdf_ > best_logpdf
                best_cluster = icluster
                best_logpdf = logpdf_
            end
        end

        # Add to cluster
        assignments_[ispike] = best_cluster
        if best_cluster == 0
            push!(noise_indexes, ispike)
        else
            addspike!(clusters[best_cluster], cur_wf, ispike)
        end
    end

    @show connections
    params.plot && finish!(hp)
    if params.graph_file != ""
        writegraph(params.graph_file, connections, init_nclusters)
    end
    assignments_
end

function plot_assignments(wf, times, assignments_, time_step::TimePeriod=Minute(10); movie_file="")
    hp = HoboPlotter(movie_file)
    time_step_seconds = Int(Second(time_step))
    for t = 0:time_step_seconds:last(times)
        rg = searchsortedfirst(times, t):searchsortedlast(times, t+time_step_seconds)
        plotclusters(hp, wf[:, rg], sub(assignments_, rg), times[last(rg)]-times[first(rg)])
    end
end

end # module


text:
using Base64
using Sockets

function main(tries)
  str_size = 131072
  str = repeat("a", str_size)

  str2 = base64encode(str)
  print("encode $(str[1:4])... to $(str2[1:4]): ")

  t = time()
  s = 0
  for i in 1:tries
    str2 = base64encode(str)
    s += length(str2)
  end
  print("$s, $(time() - t)\n")

  str3 = base64decode(str2)
  print("decode $(str2[1:4])... to $(String(str3[1:4])): ")

  t = time()
  s = 0
  for i in 1:tries
    str3 = base64decode(str2)
    s += length(str3)
  end
  print("$s, $(time() - t)\n")
end

println("JIT warming up")
main(5)

println("bench")
try
  socket = connect("localhost", 9001)
  write(socket, "Julia")
  close(socket)
catch
  # standalone usage
end

x = @timed main(8192)
println("Elapsed: $(x[2]), Allocated: $(x[3]), GC Time: $(x[4])")

max_stars_repo_name: wmeddie/benchmarks

max_issues_repo_name: wmeddie/benchmarks

max_forks_repo_name: wmeddie/benchmarks


text:
"""
# script bench_eq

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-08-20

Script to microbenchmark model equations.
"""

using DrWatson
@quickactivate "MEngProject"
using MEngProject,
    CUDA,
    DifferentialEquations,
    PyPlot,
    NNlib,
    ImageFiltering,
    Images,
    MEngProject.LaminartKernels,
    MEngProject.LaminartInitFunc,
    MEngProject.Utils,
    BenchmarkTools,
    Test

using OrdinaryDiffEq,
    ParameterizedFunctions, LSODA, Sundials, DiffEqDevTools, Noise

batch = 1

files = readdir(datadir("img"))

tspan = (0.0f0, 800f0)

batch_ = string(batch, "_", rand(1000:9999))
mkdir(plotsdir(string("bench_eq", batch_)))
file = "kan_sq_cont_l.png"

p = LaminartInitFunc.parameterInit_conv_gpu(
    datadir("img", file),
    Parameters.parameters_f32,
);

u0 = cu(reshape(
    zeros(Float32, p.dim_i, p.dim_j * (5 * p.K + 2)),
    p.dim_i,
    p.dim_j,
    5 * p.K + 2,
    1,
))

arr1 = similar(u0[:, :, 1:2, :])
arr2 = similar(u0[:, :, 1:1, :])

f = LaminartGPU.LamFunction_equ(
    similar(arr1), #x
    similar(arr1), #m
    similar(arr1), #s
    similar(arr2), #x_lgn,
    similar(arr1), #C,
    similar(arr1), #H_z,
    similar(arr1), # dy_temp,
    similar(arr1), # dm_temp,
    similar(arr1), # dz_temp,
    similar(arr1), # ds_temp,
    similar(arr2), # dv_temp,
    similar(arr1), # H_z_temp,
    similar(arr2), #  V_temp_1,
    similar(arr2), #  V_temp_2,
    similar(arr1), #  Q_temp,
    similar(arr1), #   P_temp
);
prob = ODEProblem(f, u0, tspan, p)
bm = @benchmark solve(prob)
push!(benches, bm)