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

using Markdown
using InteractiveUtils

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

# ╔═╡ 42b4e6fe-a897-11eb-2676-f7fd96f35a22
begin
	using NeidSolarScripts
	using RvSpectMLBase
	using EchelleInstruments
	using EchelleCCFs
	using RvSpectML
	using DataFrames, Query
	using CSV, JLD2, FileIO
	using Dates
	using PlutoUI
	using Plots, ColorSchemes
	using Statistics, NaNMath
	using Polynomials
end

# ╔═╡ 87597e02-19c7-418a-8eb3-ad25695bb5bf
md"""
# Process Daily NEID Solar CCFs
"""

# ╔═╡ ba3a32ee-9664-4279-b27b-0fed992e4d1c
md"## Code "

# ╔═╡ f64f9468-25ee-4799-b776-b1cfbe5b7f57
md"""
### Parameters for processing
"""

# ╔═╡ c259acbf-da3d-4b73-af2b-8fc22f2548d7
begin
	ccf_path = "/mnt/data_simons/NEID/v0.7develop20210501/outputs3"
	daily_filename = "daily_ccfs.jld2"
	max_days_to_use = 365
	order_idx_to_use_for_snr = 60
	order_snr_threshold_max = 36779.37219973966
	min_obs_in_day = 5
end;

# ╔═╡ 15741d2d-5fc2-4dc0-baf1-3ba3a72b5399
begin  # Parameters for picking which days to calculate RVs for
	order_snr_threshold = order_snr_threshold_max*0.94
	solar_hour_angle_threshold = 2.0
end;

# ╔═╡ d6219a57-7c1d-4b23-9b80-ae96f005b762
begin   # Parameters for building CCF template
	solar_hour_angle_threshold_template = 1.0
	order_snr_threshold_template = order_snr_threshold_max*0.96
end;

# ╔═╡ 60559545-8f85-4027-8471-4ce3d7cc952a
md"""
### Script
"""

# ╔═╡ 21d22e3c-a56a-4c03-a9aa-226acea40c72
md"""
# Results
## One Day
"""

# ╔═╡ fe0b7b47-39bb-45f8-873e-9a4a07d7d1b3
md"""
Order RV slope limits (lo, hi) $(@bind order_slope_lo Slider(-100.0:0.0; default=-20))  $(@bind order_slope_hi Slider(0.0:100.0; default=20.0))

"""

# ╔═╡ b6da427f-0c72-4de7-bdca-b85688ac9bed
md"""
## Combined RVs
#### Order indices to include in RV calculation
1 $(@bind use_ord_1 CheckBox(default=true))
2 $(@bind use_ord_2 CheckBox(default=false))
3 $(@bind use_ord_3 CheckBox(default=false))
4 $(@bind use_ord_4 CheckBox(default=true))
5 $(@bind use_ord_5 CheckBox(default=true))
6 $(@bind use_ord_6 CheckBox(default=true))
7 $(@bind use_ord_7 CheckBox(default=false))
8 $(@bind use_ord_8 CheckBox(default=true))
9 $(@bind use_ord_9 CheckBox(default=true))
10 $(@bind use_ord_10 CheckBox(default=true))
11 $(@bind use_ord_11 CheckBox(default=false))
12 $(@bind use_ord_12 CheckBox(default=true))
13 $(@bind use_ord_13 CheckBox(default=false))
14 $(@bind use_ord_14 CheckBox(default=false))
15 $(@bind use_ord_15 CheckBox(default=false))
16 $(@bind use_ord_16 CheckBox(default=false))
17 $(@bind use_ord_17 CheckBox(default=false))
18 $(@bind use_ord_18 CheckBox(default=true))
19 $(@bind use_ord_19 CheckBox(default=true))
20 $(@bind use_ord_20 CheckBox(default=true))

21 $(@bind use_ord_21 CheckBox(default=true))
22 $(@bind use_ord_22 CheckBox(default=true))
23 $(@bind use_ord_23 CheckBox(default=false))
24 $(@bind use_ord_24 CheckBox(default=false))
25 $(@bind use_ord_25 CheckBox(default=false))
26 $(@bind use_ord_26 CheckBox(default=false))
27 $(@bind use_ord_27 CheckBox(default=false))
28 $(@bind use_ord_28 CheckBox(default=false))
29 $(@bind use_ord_29 CheckBox(default=false))
30 $(@bind use_ord_30 CheckBox(default=false))
31 $(@bind use_ord_31 CheckBox(default=false))
32 $(@bind use_ord_32 CheckBox(default=false))
33 $(@bind use_ord_33 CheckBox(default=false))
34 $(@bind use_ord_34 CheckBox(default=false))
35 $(@bind use_ord_35 CheckBox(default=false))
36 $(@bind use_ord_36 CheckBox(default=false))
37 $(@bind use_ord_37 CheckBox(default=false))
38 $(@bind use_ord_38 CheckBox(default=false))
39 $(@bind use_ord_39 CheckBox(default=false))
40 $(@bind use_ord_40 CheckBox(default=false))

41 $(@bind use_ord_41 CheckBox(default=false))
42 $(@bind use_ord_42 CheckBox(default=false))
43 $(@bind use_ord_43 CheckBox(default=false))
44 $(@bind use_ord_44 CheckBox(default=false))
45 $(@bind use_ord_45 CheckBox(default=false))
46 $(@bind use_ord_46 CheckBox(default=false))
47 $(@bind use_ord_47 CheckBox(default=false))
48 $(@bind use_ord_48 CheckBox(default=false))
49 $(@bind use_ord_49 CheckBox(default=false))
50 $(@bind use_ord_50 CheckBox(default=false))
51 $(@bind use_ord_51 CheckBox(default=false))
52 $(@bind use_ord_52 CheckBox(default=false))
53 $(@bind use_ord_53 CheckBox(default=false))
Manual override $(@bind use_ord_manual CheckBox(default=false))
Highlight $(@bind order_to_highlight NumberField(0:53; default=8))
"""

# ╔═╡ b94f222a-0c35-49d7-94f2-a2347efb6b2d
order_idx_for_rvs_manual = vcat(1,4:6,10:11,19:23);

# ╔═╡ e7541e93-8b7a-4c26-901e-d887931b8103
begin
	if use_ord_manual
		order_idx_for_rvs = order_idx_for_rvs_manual
	else
		order_idx_for_rvs = Int64[] #order_idx_for_rvs
	if use_ord_1 push!(order_idx_for_rvs,1) end
	if use_ord_2 push!(order_idx_for_rvs,2) end
	if use_ord_3 push!(order_idx_for_rvs,3) end
	if use_ord_4 push!(order_idx_for_rvs,4) end
	if use_ord_5 push!(order_idx_for_rvs,5) end
	if use_ord_6 push!(order_idx_for_rvs,6) end
	if use_ord_7 push!(order_idx_for_rvs,7) end
	if use_ord_8 push!(order_idx_for_rvs,8) end
	if use_ord_9 push!(order_idx_for_rvs,9) end
	if use_ord_10 push!(order_idx_for_rvs,10) end
	if use_ord_11 push!(order_idx_for_rvs,11) end
	if use_ord_12 push!(order_idx_for_rvs,12) end
	if use_ord_13 push!(order_idx_for_rvs,13) end
	if use_ord_14 push!(order_idx_for_rvs,14) end
	if use_ord_15 push!(order_idx_for_rvs,15) end
	if use_ord_16 push!(order_idx_for_rvs,16) end
	if use_ord_17 push!(order_idx_for_rvs,17) end
	if use_ord_18 push!(order_idx_for_rvs,18) end
	if use_ord_19 push!(order_idx_for_rvs,19) end
	if use_ord_20 push!(order_idx_for_rvs,20) end
	if use_ord_21 push!(order_idx_for_rvs,21) end
	if use_ord_22 push!(order_idx_for_rvs,22) end
	if use_ord_23 push!(order_idx_for_rvs,23) end
	if use_ord_24 push!(order_idx_for_rvs,24) end
	if use_ord_25 push!(order_idx_for_rvs,25) end
	if use_ord_26 push!(order_idx_for_rvs,26) end
	if use_ord_27 push!(order_idx_for_rvs,27) end
	if use_ord_28 push!(order_idx_for_rvs,28) end
	if use_ord_29 push!(order_idx_for_rvs,29) end
	if use_ord_30 push!(order_idx_for_rvs,30) end
	if use_ord_31 push!(order_idx_for_rvs,31) end
	if use_ord_32 push!(order_idx_for_rvs,32) end
	if use_ord_33 push!(order_idx_for_rvs,33) end
	if use_ord_34 push!(order_idx_for_rvs,34) end
	if use_ord_35 push!(order_idx_for_rvs,35) end
	if use_ord_36 push!(order_idx_for_rvs,36) end
	if use_ord_37 push!(order_idx_for_rvs,37) end
	if use_ord_38 push!(order_idx_for_rvs,38) end
	if use_ord_39 push!(order_idx_for_rvs,39) end
	if use_ord_40 push!(order_idx_for_rvs,40) end
	if use_ord_41 push!(order_idx_for_rvs,41) end
	if use_ord_42 push!(order_idx_for_rvs,42) end
	if use_ord_43 push!(order_idx_for_rvs,43) end
	if use_ord_44 push!(order_idx_for_rvs,44) end
	if use_ord_45 push!(order_idx_for_rvs,45) end
	if use_ord_46 push!(order_idx_for_rvs,46) end
	if use_ord_47 push!(order_idx_for_rvs,47) end
	if use_ord_48 push!(order_idx_for_rvs,48) end
	if use_ord_49 push!(order_idx_for_rvs,49) end
	if use_ord_50 push!(order_idx_for_rvs,50) end
	if use_ord_51 push!(order_idx_for_rvs,51) end
	if use_ord_52 push!(order_idx_for_rvs,52) end
	if use_ord_53 push!(order_idx_for_rvs,53) end
	end
	order_idx_for_rvs
end;

# ╔═╡ 473302a3-2541-49a1-8f21-906bc2c9d279
md"""
Order RV RMS limits (lo, hi)
$(@bind order_rv_rms_lo Slider(0.0:0.1:2.0; default=0))
$(@bind order_rv_rms_hi Slider(0.5:0.2:200.0; default=10))
"""

# ╔═╡ 9833f763-ece1-4739-a6d1-0f3182d32d5d
md"### Checking for Chromatic Effects"

# ╔═╡ 2c88f409-6763-41d8-b1cc-168df855728b
md"""
# Implementation Details
"""

# ╔═╡ 38bb2810-3a30-4a39-99dd-83518c093be2
md"""
## Functions used above
"""

# ╔═╡ c23c4063-947e-436a-910b-c43afdb9c244
function make_dataframe_of_ccf_files(path::String,daily_filename::String)
	datafiles = DataFrame(:datestr=>String[],:filename=>String[])
	for (root, dirs, files) in walkdir(path)
	    for file in files
			fn = joinpath(root, file)
			if !(isfile(fn) && filesize(fn)>0) continue end
			m = match(r"(\d+)/" * daily_filename,fn)
			if m == nothing continue end
			datestr = m.captures[1]
			d = Dict(:datestr=>datestr, :filename=>fn)
			push!(datafiles,d)
	    end
	end
	return datafiles
end

# ╔═╡ 3c4b1028-bcf7-45b0-bba2-dab5a63724d0
datafiles = make_dataframe_of_ccf_files(ccf_path,daily_filename);

# ╔═╡ 26caf850-c46e-453c-aad9-c7bdb71f858b
function select_days_to_analyze_add_manifest_new(df::DataFrame; solar_hour_angle_threshold::Real = 3.0, order_snr_threshold::Real = order_snr_threshold_max*0.9, order_idx_to_use_for_snr::Integer = 60, min_obs_in_day::Integer = 5, max_days_to_use::Integer = 10)
	df[!,:nobs_in_file] .= 0
	#df[!,:nobs_use] .= 0
	df[!,:analyze_day] .= false
	df[!,:manifest] = fill(DataFrame(),size(df,1))

	num_days_to_use_found = 0
	for (i,day) in enumerate(eachrow(df))
		solar_info_fn = joinpath(dirname(day.filename),"..",day.datestr * "_solar_info.csv")
		println("Looking for ", solar_info_fn)
		@assert isfile(solar_info_fn) && filesize(solar_info_fn)>0
		df_solarinfo = CSV.read(solar_info_fn, DataFrame)
		manifest_all = load(day.filename,"manifest")
		df_combo = innerjoin(manifest_all, df_solarinfo, on = :bjd, makeunique = true)
		println("# Sizes: manifest= ", size(manifest_all,1), " solarinfo= ", size(df_solarinfo), " df_combo= ", size(df_combo,1))
		#orders_to_use = load(day.filename,"orders_to_use")
		#=
		manifest_use = manifest_all |>
			@filter(abs(_.solar_hour_angle)<= solar_hour_angle_threshold) |>
			@filter(_.order_snrs[order_idx_to_use_for_snr] >= order_snr_threshold ) |>
				DataFrame
		=#
		manifest_use = df_combo |>
			@filter(abs(_.solar_hour_angle)<= solar_hour_angle_threshold) |>
			#@filter(_.order_snrs[order_idx_to_use_for_snr] >= order_snr_threshold ) |>
			@filter( 0.995 <= _.pyrheliometer_flux_over_model <= 1.015 ) |>
			@filter( _.pyrheliometer_rms_flux <= 0.005 ) |>
				DataFrame

		day.nobs_in_file = size(manifest_all,1)
		#day.
		nobs_use = size(manifest_use,1)
		#day.manifest = manifest_all
		day.manifest = df_combo
		#day.manifest[!,:pyrheliometer_flux_over_model] = manifest_use[!,pyrheliometer_flux_over_model]
		#day.manifest[!,:pyrheliometer_rms_flux] = manifest_use[!,:pyrheliometer_rms_flux]
		if nobs_use <= min_obs_in_day
			continue
		end
		day.analyze_day = true
		num_days_to_use_found += 1
		if num_days_to_use_found >= max_days_to_use
			break
		end

	end
	df_out = df |> @filter(_.analyze_day) |> @orderby_descending(_.datestr)|>  @take(max_days_to_use) |> @orderby(_.datestr) |> DataFrame
 	return df_out
end

# ╔═╡ 785a2042-bb12-4a2a-8117-c2cfb396d29c
#datafiles_use = select_days_to_analyze_add_manifest(datafiles,
datafiles_use = select_days_to_analyze_add_manifest_new(datafiles,
					max_days_to_use=max_days_to_use,
					solar_hour_angle_threshold=solar_hour_angle_threshold,
					order_snr_threshold=order_snr_threshold);

# ╔═╡ 433c64f5-3f0c-4442-a5ea-5d57e14a2038
size(datafiles_use)

# ╔═╡ 4ae8bce1-5890-4647-94e5-e0dfd14e9eb5
names(datafiles_use)

# ╔═╡ 188cf9ab-f839-490b-99fa-04611c355bfe
datafiles_use

# ╔═╡ 0e1f6708-53d8-4cca-acf3-e7b2cac1dd7a
begin   # Get size of order ccfs
	v_grid =  load(first(datafiles_use.filename),"v_grid")
	num_v = length(v_grid)
	ex_order_ccfs = load(first(datafiles_use.filename),"order_ccfs")
	num_order_idx = size(ex_order_ccfs,2)
end;

# ╔═╡ ebf32f98-3ba0-4650-83ca-4d7a5b310c05
md"""
## Orders RVs
Order index (min, max, step)
$(@bind order_idx_lo NumberField(1:num_order_idx; default=4))
$(@bind order_idx_hi NumberField(1:num_order_idx; default=14))
$(@bind order_idx_step NumberField(1:10; default=1))
Use order range? $(@bind use_order_range CheckBox(default=true))
Plot bad days? $(@bind plt_all_days CheckBox(default=false))
"""

# ╔═╡ c487c08a-93a1-42dc-b465-3268be5954df
datafiles_use[79,:manifest]

# ╔═╡ bcbd90c4-60b1-4c90-aa4b-69b0634e2fcd
datafiles_use.manifest[13]

# ╔═╡ 0b6ec47b-3610-453f-b0c3-a99d5fb83ab4
histogram(datafiles_use.manifest[1][!,:pyrheliometer_rms_flux],bins=range(0,stop=0.006,length=40))

# ╔═╡ c1bf8bf8-fb05-4731-8140-c0f6332496bd
histogram(datafiles_use.manifest[1][!,:pyrheliometer_flux_over_model],bins=range(0.995,stop=1.015,length=40))

# ╔═╡ b957ac39-458f-48b6-8683-229d09ecaea4
function select_days_to_analyze_add_manifest(df::DataFrame; solar_hour_angle_threshold::Real = 3.0, order_snr_threshold::Real = order_snr_threshold_max*0.9, order_idx_to_use_for_snr::Integer = 60, min_obs_in_day::Integer = 5, max_days_to_use::Integer = 10)
	df[!,:nobs_in_file] .= 0
	#df[!,:nobs_use] .= 0
	df[!,:analyze_day] .= false
	df[!,:manifest] = fill(DataFrame(),size(df,1))

	num_days_to_use_found = 0
	for (i,day) in enumerate(eachrow(df))
		manifest_all = load(day.filename,"manifest")
		#orders_to_use = load(day.filename,"orders_to_use")
		manifest_use = manifest_all |>
			@filter(abs(_.solar_hour_angle)<= solar_hour_angle_threshold) |>
			@filter(_.order_snrs[order_idx_to_use_for_snr] >= order_snr_threshold ) |>
				DataFrame
		day.nobs_in_file = size(manifest_all,1)
		#day.
		nobs_use = size(manifest_use,1)
		day.manifest = manifest_all
		if nobs_use <= min_obs_in_day
			continue
		end
		day.analyze_day = true
		num_days_to_use_found += 1
		if num_days_to_use_found >= max_days_to_use
			break
		end

	end
	df_out = df |> @filter(_.analyze_day) |> @orderby_descending(_.datestr)|>  @take(max_days_to_use) |> @orderby(_.datestr) |> DataFrame
 	return df_out
end

# ╔═╡ d9877f46-c461-4615-bb0c-0c7541238899
function compute_mean_order_ccf(df::DataFrame; solar_hour_angle_threshold::Real = 1.0, order_snr_threshold::Real = order_snr_threshold_max*0.95, order_idx_to_use_for_snr::Integer = 60 )
	ex_order_ccfs = load(first(df.filename),"order_ccfs")
	mean_order_ccf = zeros(size(ex_order_ccfs,1),size(ex_order_ccfs,2))
	weights_order_ccf = zeros(size(ex_order_ccfs,1),size(ex_order_ccfs,2))
	for row in eachrow(df)
		sha_mask = abs.(row.manifest.solar_hour_angle).<= solar_hour_angle_threshold
	    snr_mask = map(i->row.manifest.order_snrs[i][order_idx_to_use_for_snr],1:length(row.manifest.order_snrs)) .>= order_snr_threshold
		mask = sha_mask .& snr_mask
		nobs_use = sum(mask)
		v_grid_this =  load(row.filename,"v_grid")
		@assert all(v_grid .== v_grid_this)
		order_ccfs = load(row.filename,"order_ccfs")
		order_ccf_vars = load(row.filename,"order_ccf_vars")
		@assert size(order_ccfs,1) == size(mean_order_ccf,1)
		@assert size(order_ccfs,2) == size(mean_order_ccf,2)
		num_obs = size(order_ccfs,3)
		weights = zeros(num_obs)
		for i in 1:num_v
			for j in 1:num_order_idx
				weights .= 1.0 ./ order_ccf_vars[i,j,:].^2
				mean_order_ccf[i,j] += sum(order_ccfs[i,j,mask] .* weights[mask])
				weights_order_ccf[i,j] += sum(weights[mask])
			end
		end
	end
	mean_order_ccf ./= weights_order_ccf
	return mean_order_ccf
end

# ╔═╡ de7e3ebb-f9d2-4830-9bec-f034d965ab76
mean_order_ccf = compute_mean_order_ccf(datafiles_use, solar_hour_angle_threshold=solar_hour_angle_threshold_template, order_snr_threshold=order_snr_threshold_template, order_idx_to_use_for_snr=order_idx_to_use_for_snr);

# ╔═╡ b12592bf-ec54-435f-b795-183cf84bcb4c
begin
	mrv = []
	for ord_idx in 1:num_order_idx
		templ = vec(mean_order_ccf[:,ord_idx])
		if all(isnan.(templ))
			push!(mrv,missing)
			continue
		end
		mrv_this_order = RVFromCCF.MeasureRvFromCCFTemplate(v_grid=v_grid,template=templ, frac_of_width_to_fit = 1.0, measure_width_at_frac_depth=0.5)
		push!(mrv,mrv_this_order)
	end
	mrv
end

# ╔═╡ 93465fa0-83bc-413b-b46e-7ed0d3d0d1a7
function calc_all_rvs(fn::String) #, mask::BitVector)
	v_grid_this =  load(fn,"v_grid")
	@assert all(v_grid .== v_grid_this)
	order_ccfs = load(fn,"order_ccfs")
	order_ccf_vars = load(fn,"order_ccf_vars")

	num_obs = size(order_ccfs,3) #sum(mask)
	@assert num_obs <= size(order_ccfs,3)
	rv_matrix = zeros(num_order_idx, num_obs)
	σ_rv_matrix = zeros(num_order_idx, num_obs)
	for (i,obs_idx) in enumerate((1:size(order_ccfs,3))) #[mask])
		for ord_idx in 1:num_order_idx
			if typeof(mrv[ord_idx]) == Missing  continue end
			(rv, σ_rv) = mrv[ord_idx](v_grid,order_ccfs[:,ord_idx,obs_idx],order_ccf_vars[:,ord_idx,obs_idx])
			rv_matrix[ord_idx,i] = rv
			σ_rv_matrix[ord_idx,i] = σ_rv
		end
	end
	return (rvs=rv_matrix, σ_rvs = σ_rv_matrix, filename=fn)
end

# ╔═╡ 91d634c4-e4ec-49b8-a035-5df0740f9fb0
all_rvs = DataFrame(map(i->calc_all_rvs(datafiles_use.filename[i]),1:length(datafiles_use.filename)));

# ╔═╡ e00fc7cb-05a7-4215-b7f6-5bad647ab090
df_rvs = innerjoin(datafiles_use, all_rvs, on = :filename);

# ╔═╡ 022f7788-a2d0-4489-af0b-71692f800fdf
function find_good_day_idx(df::DataFrame; solar_hour_angle_threshold::Real = 1.0, order_snr_threshold::Real = order_snr_threshold_max*0.96, order_idx_to_use_for_snr::Integer = 60)
	argmax(
		df |> @map( sum((abs.(_.manifest.solar_hour_angle).<solar_hour_angle_threshold) .& map(i->_.manifest.order_snrs[i][order_idx_to_use_for_snr] > order_snr_threshold,1:size(_.manifest,1))) ) |> collect )
end

# ╔═╡ 8cb9eca6-2ebe-4faa-b6c6-5ea267e64dc9
obs_idx_good_day = find_good_day_idx(datafiles_use)

# ╔═╡ 34bda3ea-4238-4670-9f6b-9207f63047cb
function calc_order_weights(df::DataFrame, obs_idx::Integer)
	order_weights = 1.0 ./ df.σ_rvs[obs_idx_good_day]
	order_weights[isinf.(order_weights)] .= 0
	order_weights[isnan.(order_weights)] .= 0
	order_weights = vec(sum(order_weights,dims=2))
	order_weights ./= sum(order_weights)
	return order_weights
end

# ╔═╡ 1174095b-a339-47b2-9f21-3f5b2c52b4eb
order_weights = calc_order_weights(all_rvs,obs_idx_good_day)

# ╔═╡ 72522f83-2b3c-4e40-a1de-c3508d50fdab
function summarize_day_orders(day_to_plt::Integer; df::DataFrame, orders, weights)
	local δrvs = df.manifest[day_to_plt][!,:Δv_diff_ext]
	local daily_means = map(ord->NaNMath.mean(df.rvs[day_to_plt][ord,:].+δrvs),orders)
	local daily_rmss = map(ord->NaNMath.std(df.rvs[day_to_plt][ord,:].+δrvs),orders)
	local times = df.manifest[day_to_plt].bjd

	local daily_rmss_binned =  map(ord->
		NaNMath.std(bin_times_and_rvs_max_Δt(times=times,
			rvs=df.rvs[day_to_plt][ord,:].+δrvs,
				Δt_threshold= 6.0/(24*60))[2] ),orders)

	local shas = df.manifest[day_to_plt].solar_hour_angle ./24
	local tsha0 = roots(Polynomials.fit(times.-2.459e6,shas,1))[1].+2.459e6
	local daily_fits = map(ord->Polynomials.fit(times.-tsha0,df.rvs[day_to_plt][ord,:].+δrvs,1),orders)
	local daily_rv_sha0s = map(i->daily_fits[i][0],1:length(orders))
	local daily_slopes = map(i->daily_fits[i][1],1:length(orders))

	return Dict(:mean_order_rv=>daily_means,
		:rms_order_rv=>daily_rmss, :t_sha0=>tsha0, :Δt=>times.-tsha0,
		:order_rv_slope=>daily_slopes, :order_rv_sha0=>daily_rv_sha0s,
		:rms_order_rv_binned=>daily_rmss_binned )
end

# ╔═╡ 4c62a111-112c-4b8d-979d-dea84e35232a
day_orders = DataFrame(summarize_day_orders.(1:size(datafiles_use,1),df=df_rvs,orders=1:num_order_idx, weights= order_weights))

# ╔═╡ 2d56d9c0-b23a-4399-a86e-77edbf25d8b5
median_order_rv = map(ord->NaNMath.median(map(i->day_orders.mean_order_rv[i][ord],(1:size(day_orders,1)))),1:num_order_idx)

# ╔═╡ 701ebc68-2809-4544-be85-c4fadae6c391
begin
	daily_rms_across_orders = map(obsid->NaNMath.std((day_orders.order_rv_sha0[obsid].-median_order_rv)[vcat(1,4:6,8:10,12)]),1:size(day_orders,1))
end

# ╔═╡ 6fb4b458-6192-4000-8efd-f6fb02488861
function summarize_day_total(day_to_plt::Integer; df::DataFrame, orders, weights)
	rv_tot = vec(sum(df.rvs[day_to_plt][orders,:].*weights[orders],dims=1) ./ 					   sum(weights[orders],dims=1))
	rv_tot .+= df.manifest[day_to_plt][!,:Δv_diff_ext]
	times = df.manifest[day_to_plt].bjd
	shas = df.manifest[day_to_plt].solar_hour_angle ./24
	tsha0 = roots(Polynomials.fit(times.-2.459e6,shas,1))[1].+2.459e6

	sha_mask = abs.(df.manifest[day_to_plt].solar_hour_angle).<= solar_hour_angle_threshold
	snr_mask = map(i->df.manifest[day_to_plt].order_snrs[i][order_idx_to_use_for_snr],1:length(df.manifest[day_to_plt].order_snrs)) .>= order_snr_threshold
	pyrh_mask = (0.995 .<= df.manifest[day_to_plt].pyrheliometer_flux_over_model .<= 1.015 ) .& ( df.manifest[day_to_plt].pyrheliometer_rms_flux .<= 0.006 )
	#mask = sha_mask .& snr_mask
	mask = sha_mask .& pyrh_mask
	if sum(mask) >= 1
		rms_rv = std(rv_tot[mask])
		(t_binned, rvs_binned, num_obs_in_bin) =
			bin_times_and_rvs_max_Δt(times=times[mask],rvs=rv_tot[mask],Δt_threshold= 6.01/(24*60))
		mask_bin_full = num_obs_in_bin .== 4
		rms_rv_binned = sum(mask_bin_full)>3 ? std(rvs_binned[mask_bin_full]) : NaN
		daily_fit = Polynomials.fit(times[mask].-tsha0,rv_tot[mask],1)
		rms_rv_binned_fit =  sum(mask_bin_full)>3 ? std(rvs_binned[mask_bin_full] .-daily_fit.(t_binned[mask_bin_full].-tsha0)) : NaN
		fit_constant = daily_fit[0]
		fit_slope = daily_fit[1]
		num_obs_in_bin = zeros(Int64,0)
	else
		rms_rv = NaN
		t_binned = zeros(0)
		rvs_binned = zeros(0)
		num_obs_in_bin = zeros(Int64,0)
		rms_rv_binned = NaN
		rms_rv_binned_fit = NaN
		fit_constant = NaN
		fit_slope = NaN
	end

	Dict( :rv=>rv_tot, :rms_rv=>rms_rv,
		  :t_sha0=>tsha0, :Δt => times.-tsha0,
		  :rv_slope => fit_slope, :rv_sha0 =>fit_constant,
		  :Δt_binned => t_binned.-tsha0, :rv_binned=>rvs_binned, :num_obs_in_bin=>num_obs_in_bin,
		  :rms_rv_binned=>rms_rv_binned, :rms_rv_binned_fit =>rms_rv_binned_fit,
		  :mask=>mask
		)
end

# ╔═╡ 343f3226-a6e1-4bf5-a816-f92cb74cf434
begin
	day_comb = DataFrame(summarize_day_total.(1:size(datafiles_use,1),df=df_rvs,orders=order_idx_for_rvs, weights= order_weights));
end;

# ╔═╡ 592b8b00-91a0-4ab4-b26b-9fc358978e9a
md"""
Day to plot
$(@bind obs_idx_plt NumberField(1:size(day_comb,1); default=obs_idx_good_day))
"""

# ╔═╡ c863a4fc-270b-4473-beaa-d34ba2d719de
begin
	mask_decent_days = (length.(day_comb.rv_binned).>=20) #.&
		#(daily_rms_across_orders .< 6) #.&
	 	#(map(obsidx->median(map(ordidx->day_orders[obsidx,:rms_order_rv][ordidx],order_idx_for_rvs)),1:size(day_orders,1)) .< 6)
	all_days_to_plt = plt_all_days ? (1:size(day_orders,1)) : (1:size(day_orders,1))[mask_decent_days]
	num_all_days_to_plt = length(all_days_to_plt)
end

# ╔═╡ fb17d03b-880c-4be6-bd50-f31bdd3e2fd4
md"""
Date range (min, max, step)
$(@bind days_to_plt_idx_lo Slider(1:num_all_days_to_plt; default=1))
$(@bind days_to_plt_idx_hi Slider(1:num_all_days_to_plt; default=num_all_days_to_plt))

"""

# ╔═╡ bbf59177-36a8-480c-af39-1998b37032d2
begin
		default_order_rv_lo = plt_all_days ? -500 : -6
		default_order_rv_hi = plt_all_days ? 500 : 6
		default_order_rv_step = plt_all_days ? 1 : 0.1
		days_to_plt = all_days_to_plt[days_to_plt_idx_lo:days_to_plt_idx_hi]
end;

# ╔═╡ a79fb3f3-a07c-43b8-9be6-72e7347ff40e
size(days_to_plt)

# ╔═╡ 92175df7-b060-4f6f-8acb-265bdbe6713a
begin
	max_order_daily_rms_hi = NaNMath.maximum(daily_rms_across_orders[days_to_plt])
	default_order_daily_rms_hi = 10 #NaNMath.maximum(daily_rms_across_orders[mask_decent_days])
end;

# ╔═╡ 8af9d40f-8d50-4376-82e4-01981a54893f
md"""
Order RV RMS limits (lo, hi)
$(@bind order_daily_rms_lo Slider(0.0:0.2:5; default=0.1))
$(@bind order_daily_rms_hi Slider(0.0:0.2:30; default=default_order_daily_rms_hi))

"""

# ╔═╡ f7db71f7-9bac-4aad-a636-5c60ae4c259b
begin
	local plt = plot()
	local orders_to_plt = use_order_range ? (order_idx_lo:order_idx_step:order_idx_hi) : order_idx_for_rvs
	local pal = palette(:berlin, length(orders_to_plt))
	for (i,ordidx) in enumerate(orders_to_plt)
		scatter!(plt,Date.(Dates.julian2datetime.(day_orders.t_sha0[days_to_plt])),	map(obsidx->day_orders[obsidx,:rms_order_rv][ordidx],1:size(day_orders,1)),ms=1.0,label=:none,c=pal[i])
		scatter!(plt,Date.(Dates.julian2datetime.(day_orders.t_sha0[days_to_plt])),	map(obsidx->day_orders[obsidx,:rms_order_rv][ordidx],days_to_plt),label=:none,c=pal[i])
	end
	xlabel!(plt,"Time (d)")
	ylabel!(plt,"RMS of Orders RVs within each day")
	ylims!(order_daily_rms_lo,order_daily_rms_hi)
	title!(plt,"Daily RMS of Order RVs")
end

# ╔═╡ 9dc77012-f845-4a5b-a059-fd13c1308b3d
begin
	histogram(daily_rms_across_orders,
		bins=range(order_daily_rms_lo,stop=order_daily_rms_hi,length=20),alpha=0.25, label="All usable days")
	histogram!(daily_rms_across_orders[days_to_plt],
		bins=range(order_daily_rms_lo,stop=order_daily_rms_hi,length=20),alpha=0.5, label="Good days")
	xlabel!("RMS (m/s)")
	ylabel!("Count")
	title!("RMS in daily RV over orders used for RVs")
	xlims!(order_daily_rms_lo,order_daily_rms_hi)
end

# ╔═╡ f8e4b341-1613-4993-a193-437f2adf4a21
md"""
Order RV limits (lo, hi) $(@bind order_rv_lo Slider(-500.0:0.0; default=default_order_rv_lo))  $(@bind order_rv_hi Slider(0.0:500.0; default=default_order_rv_hi))
"""

# ╔═╡ 0d9dd1ea-9ef8-460c-8d95-20a46eb1ebb6
begin
	local plt1 = plot(legend=:none)
	local t_day =  map(i->day_orders[i,:t_sha0],1:size(day_orders,1))
	local orders_to_plt =use_order_range ?  (order_idx_lo:order_idx_step:order_idx_hi) :  order_idx_for_rvs
	local pal = palette(:berlin, length(orders_to_plt))
	for (i,ord) in enumerate(orders_to_plt)
		ord_rv_vs_day = map(i->day_orders[i,:order_rv_sha0][ord]-median_order_rv[ord],days_to_plt)
		if all(ord_rv_vs_day .== 0) continue end
		#scatter!(plt1,Date.(Dates.julian2datetime.(t_day)),map(i->day_orders[i,:order_rv_sha0][ord]-median_order_rv[ord],1:size(day_orders,1)),c=pal[i],ms=1.0)
		scatter!(plt1,Date.(Dates.julian2datetime.(t_day[days_to_plt])),ord_rv_vs_day,c=pal[i])
	end
	title!(plt1,"Order RV at midday")
	xlabel!(plt1,"BJD")
	ylabel!(plt1,"ΔRV (m/s)")
	ylims!(plt1,order_rv_lo,order_rv_hi)
	plt1
end

# ╔═╡ cfce959c-c171-4923-b5f9-29666d9c20e4
begin
	local plt2 = plot(legend=:none)
	local t_day =  map(i->day_orders[i,:t_sha0],1:size(day_orders,1))
	local orders_to_plt = use_order_range ? (order_idx_lo:order_idx_step:order_idx_hi) : order_idx_for_rvs
	local pal = palette(:berlin, length(orders_to_plt))
	for (j,ord) in enumerate(orders_to_plt)
		ord_rv_vs_day = map(i->day_orders[i,:order_rv_slope][ord],days_to_plt)
		if all(ord_rv_vs_day .== 0) continue end
		scatter!(plt2,Date.(Dates.julian2datetime.(t_day[days_to_plt])),ord_rv_vs_day,c=pal[j])
	end
	title!(plt2,"Order RV Daily Slope")
	xlabel!(plt2,"BJD")
	ylabel!(plt2,"RV (m/s/day)")
	ylims!(plt2,order_slope_lo,order_slope_hi)
	plt2
end

# ╔═╡ 0b07ac3c-d61e-43ff-acff-f907019835d6
begin
	local plt = plot()
	local pal = palette(:berlin, length(order_idx_for_rvs))
	for (i,ord) in enumerate(order_idx_for_rvs)
		#histogram!(plt,day_orders[ord,:rms_order_rv_binned],bin=range(0.1,stop=order_rv_rms_hi,length=40), color=pal[i],label=string(ord))
		histogram!(plt,map(obsid->day_orders[obsid,:rms_order_rv_binned][ord],days_to_plt),
				bin=range(0.1,stop=order_rv_rms_hi,length=40), color=pal[i],alpha=0.4,label=string(ord))
	end
	if 1<= order_to_highlight <= size(day_orders,1)
		histogram!(plt,map(obsid->day_orders[obsid,:rms_order_rv_binned][order_to_highlight],days_to_plt),
				bin=range(0.1,stop=order_rv_rms_hi,length=40), color=:yellow,alpha=0.5,label=string(order_to_highlight))
	end
	xlabel!("RMS (m/s)")
	ylabel!("Count")
	xlims!(order_rv_rms_lo,order_rv_rms_hi)
	plt
end

# ╔═╡ 5bc2c3e1-62dc-4f64-9339-87c6cbca8a8b
size(day_comb[44,:rv_binned])

# ╔═╡ e858ca22-5ddd-45aa-8815-c3acc3f9741e
date_mask = datetime2julian(Date(2021,01,01)) .<= day_comb.t_sha0[days_to_plt] .<= datetime2julian(Date(2021,05,15));

# ╔═╡ c8621274-00c9-4f67-b285-c624d1346527
begin
	default_daily_rv_rms_lo = NaNMath.minimum(day_comb.rms_rv_binned_fit)
	default_daily_rv_rms_hi = NaNMath.maximum(day_comb.rms_rv)
end;

# ╔═╡ 9b81f986-ea3c-49f5-be04-b12b392aacd6
md"""
Daily RMS RV limits (lo, hi)
$(@bind daily_rv_rms_lo Slider(0.0:0.05:1.0; default=default_daily_rv_rms_lo))
$(@bind daily_rv_rms_hi Slider(0.0:0.2:5.0; default=1.2))
"""

# ╔═╡ 34397405-e939-464c-a384-5d6c8fb9ca83
begin
	histogram(day_comb.rms_rv[days_to_plt], bins=daily_rv_rms_lo:(round((daily_rv_rms_hi-daily_rv_rms_lo)/40,digits=3)):daily_rv_rms_hi, alpha=0.2,label="Unbinned")
	histogram!(day_comb.rms_rv_binned[days_to_plt], bins=daily_rv_rms_lo:(round((daily_rv_rms_hi-daily_rv_rms_lo)/40,digits=3)):daily_rv_rms_hi,alpha=0.6, label="Binned")
	histogram!(day_comb.rms_rv_binned_fit[days_to_plt], bins=daily_rv_rms_lo:(round((daily_rv_rms_hi-daily_rv_rms_lo)/40,digits=3)):daily_rv_rms_hi, alpha=0.4, label="Binned - linear fit")

	xlims!(daily_rv_rms_lo,daily_rv_rms_hi)
	xlabel!("RMS (m/s)")
	ylabel!("Count")
end

# ╔═╡ 7b2e6266-6dfe-47b5-8f28-3b9d76f57a42
begin
	plot(legend=:none)
	scatter!(Date.(Dates.julian2datetime.(day_comb.t_sha0)),day_comb.rms_rv,ms=1.5,color=:blue)
	scatter!(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt])),day_comb.rms_rv[days_to_plt], color=:blue)
	scatter!(Date.(Dates.julian2datetime.(day_comb.t_sha0)),day_comb.rms_rv_binned, ms=2,color=:red)
	scatter!(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt])),day_comb.rms_rv_binned[days_to_plt], color=:red)
	title!("Daily RV RMS")
	xlabel!("Day")
	ylabel!("RMS RV within day (m/s)")
	xlims!(extrema(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt]))))
	ylims!(daily_rv_rms_lo,daily_rv_rms_hi)
end

# ╔═╡ 7d63e3eb-6132-439a-bc68-05e73cd5b068
Δrv = day_comb.rv_sha0.-NaNMath.median(day_comb.rv_sha0[days_to_plt][date_mask])

# ╔═╡ d4e7987c-94bf-4194-8e16-88c2581832d1
begin
	default_daily_rv_lo = NaNMath.minimum(Δrv[days_to_plt])
	default_daily_rv_hi = NaNMath.maximum(Δrv[days_to_plt])
end;

# ╔═╡ 72e90532-02da-4d31-88c3-7766393d9f92
md"""
Daily RV limits (lo, hi)
$(@bind daily_rv_lo Slider(-20.0:0.2:0.0; default=default_daily_rv_lo))
$(@bind daily_rv_hi Slider(0.0:0.2:20.0; default=default_daily_rv_hi))

"""

# ╔═╡ fca2fa7d-99fe-44bd-a34a-1b1552a3ae48
begin
	local plt = plot(legend=:none)
	local fit_mask = date_mask .&  (-10.0 .<Δrv[days_to_plt].< 10.0 )
	long_term_daily_rv_fit = Polynomials.fit(day_comb.t_sha0[days_to_plt][fit_mask],Δrv[days_to_plt][fit_mask],1)
	scatter!(plt,Date.(Dates.julian2datetime.(day_comb.t_sha0)),Δrv,ms=2,color=:black)
	scatter!(plt,Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt])), Δrv[days_to_plt])
	local pred = long_term_daily_rv_fit.(day_comb.t_sha0)
	plot!(plt,Date.(Dates.julian2datetime.(day_comb.t_sha0)),pred)
	rms_resid = std(Δrv[days_to_plt][fit_mask].-pred[days_to_plt][fit_mask])
	rms_resid_0 = std(Δrv[days_to_plt][fit_mask])
	title!(plt,"RV fit at midday.  Linear slope: " * string(round(long_term_daily_rv_fit[1],digits=4)) * " m/s/d.\nRMS (fit): " * string(round(rms_resid,digits=3)) * "m/s"  * "  RMS (const): " * string(round(rms_resid_0,digits=3)) * "m/s")
	xlabel!(plt,"Day")
	ylabel!(plt,"RV at solar noon (m/s)")
	xlims!((minimum(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt].-1))),maximum(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt].+1)))))
	ylims!(plt,daily_rv_lo,daily_rv_hi)
	plt
end

# ╔═╡ 1c42a7ce-3e97-48b3-8ddd-89cab6177bfe
"Daily RVs:  Linear slope: " * string(round(long_term_daily_rv_fit[1],digits=4)) * " m/s/d    " * string("RMS about fit: " * string(round(rms_resid,digits=3)) * " m/s") * " Median RMS Binned RVs " * string(round(NaNMath.median(day_comb.rms_rv_binned[days_to_plt]),digits=3)) * " m/s"

# ╔═╡ 0e9b78b0-208f-4b29-96c3-4b4129416bdd
md"""
Daily RV slope limits (lo, hi)
$(@bind daily_slope_lo Slider(-20.0:0.0; default=minimum(day_comb.rv_slope[days_to_plt])))
$(@bind daily_slope_hi Slider(0.0:20.0; default=maximum(day_comb.rv_slope[days_to_plt])))

"""

# ╔═╡ 4c71c3d6-1f9e-4812-9f36-2e65a43c8207
begin
	plot(legend=:none)
	scatter!(Date.(Dates.julian2datetime.(day_comb.t_sha0)),day_comb.rv_slope,ms=2,color=:black)
	scatter!(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt])),day_comb.rv_slope[days_to_plt])
	title!("Daily RV slope")
	xlabel!("Day")
	ylabel!("RV Slope (m/s/day)")
	xlims!(extrema(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt]))))
	ylims!(daily_slope_lo,daily_slope_hi)
end

# ╔═╡ b2b81585-2144-4b37-b0a7-286c862bb812
begin
	plot(legend=:none)
	histogram(day_comb.rv_slope[days_to_plt],bins=range(daily_slope_lo,stop=daily_slope_hi,length=30),alpha=0.5,label=:none)
	#local idx_good = length.(day_comb.rv_binned[days_to_plt]).>=25
	#histogram!(day_comb.rv_slope[days_to_plt][idx_good],bins=range(daily_slope_lo,stop=daily_slope_hi,length=30),alpha=0.5,label="Many points")
	title!("Distribution of Daily RV slopes")
	xlabel!("RV Slope (m/s/day)")
	xlims!(daily_slope_lo,daily_slope_hi)
end

# ╔═╡ e9f30692-375c-49fc-bf0d-f045951ae11f
begin
	day_farred = DataFrame(summarize_day_total.(1:size(datafiles_use,1),df=df_rvs,orders=vcat(36:40,46:50), weights= order_weights))
	day_red = DataFrame(summarize_day_total.(1:size(datafiles_use,1),df=df_rvs,orders=vcat(14:22,24:25), weights= order_weights))
	day_blue = DataFrame(summarize_day_total.(1:size(datafiles_use,1),df=df_rvs,orders=vcat(1,4:6,8:10,12), weights= order_weights))
end;

# ╔═╡ 9d89d84d-f107-4a8d-be18-55af8edb546e
begin
	local plt = plot()
	local pfit = Polynomials.fit(day_comb.t_sha0[days_to_plt],Δrv[days_to_plt],1)
	Δrv_blue = day_blue.rv_sha0.-NaNMath.mean(day_blue.rv_sha0[days_to_plt])
	Δrv_red =  day_red.rv_sha0.-NaNMath.mean(day_red.rv_sha0[days_to_plt])
	Δrv_farred =  day_farred.rv_sha0.-NaNMath.mean(day_farred.rv_sha0[days_to_plt])

	#scatter!(plt,Date.(Dates.julian2datetime.(day_blue.t_sha0[days_to_plt])), Δrv_blue[days_to_plt],color=:blue)
	scatter!(plt,Date.(Dates.julian2datetime.(day_red.t_sha0[days_to_plt])), Δrv_red[days_to_plt].-Δrv_blue[days_to_plt],color=:green, label="Red-Blue")
	scatter!(plt,Date.(Dates.julian2datetime.(day_farred.t_sha0[days_to_plt])), Δrv_farred[days_to_plt].-Δrv_blue[days_to_plt],color=:red, label="FarRed-Red")
	#rms_resid = std(Δrv[days_to_plt].-pfit.(day_comb.t_sha0[days_to_plt]))
	title!(plt,"ΔRV fit at midday") # \nLinear slope: " * string(round(pfit[1],digits=4)) * " m/s/d.  RMS about fit: " * string(round(rms_resid,digits=3)) * "m/s")
	xlabel!(plt,"Day")
	ylabel!(plt,"ΔRV at solar noon (m/s)")
	xlims!(extrema(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt]))))
	ylims!(plt,-6,6)#daily_rv_lo,daily_rv_hi)
	plt
end

# ╔═╡ 833476e3-f1f1-4042-87e0-4c21ef37f7a7
begin
	plot(legend=:none)
	scatter!(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt])),day_comb.rv_slope[days_to_plt],ms=2,color=:black)
	scatter!(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt])),day_blue.rv_slope[days_to_plt],color=:blue)
	scatter!(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt])),day_red.rv_slope[days_to_plt],color=:green)
	scatter!(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt])),day_farred.rv_slope[days_to_plt],color=:red)
	title!("Daily RV slope for Order subsets")
	xlabel!("Day")
	ylabel!("RV Slope (m/s/day)")
	xlims!(extrema(Date.(Dates.julian2datetime.(day_comb.t_sha0[days_to_plt]))))
	ylims!(daily_slope_lo-2,daily_slope_hi+2)
end

# ╔═╡ 31f212b1-dac3-4c77-a340-dd2837fc8807
function plot_rvs_day(obs_idx::Integer)
	local plt = plot()
	sha_mask = abs.(datafiles_use.manifest[obs_idx].solar_hour_angle).<= solar_hour_angle_threshold
	 snr_mask = map(i->datafiles_use.manifest[obs_idx].order_snrs[i][order_idx_to_use_for_snr],1:length(datafiles_use.manifest[obs_idx].order_snrs)) .>= order_snr_threshold
	mask = sha_mask .& snr_mask
	plot!(plt,day_comb.Δt[obs_idx][mask].*24,day_comb.Δt[obs_idx][mask].*day_comb.rv_slope[obs_idx].+day_comb.rv_sha0[obs_idx], label=:none)
	scatter!(plt,day_comb.Δt[obs_idx][mask].*24,	day_comb.rv[obs_idx][mask], ms=1.5,label="Unbinned")
	(t_binned, rvs_binned) = bin_times_and_rvs_max_Δt(times=day_comb.Δt[obs_idx][mask],rvs=day_comb.rv[obs_idx][mask],Δt_threshold= 5/(24*60))
	scatter!(plt,t_binned.*24,	rvs_binned, label="Binned 5min")
	xlabel!(plt,"Δt from solar noon (hr)")
	ylabel!(plt,"RV (m/s)")
	title!(plt,datafiles_use.datestr[obs_idx] * ": RMS= " * string(round(day_comb.rms_rv_binned[obs_idx],digits=3)) * "m/s Slope= " * string(round(day_comb.rv_slope[obs_idx],digits=3)) * "m/s/d" )
	return plt
end

# ╔═╡ 8bb918eb-cdb6-4a15-a344-88ec22df92b2
plot_rvs_day(obs_idx_plt)

# ╔═╡ 4441d217-5f2e-4052-8fda-e11672ae1813
md"""
## Package Setup
"""

# ╔═╡ e0756a3b-998a-4f8d-9ab3-df50422ce14b
md"""
## UI Tweaks
"""

# ╔═╡ 44ae02ac-6f9a-4af5-8901-3428da03af22
PlutoUI.TableOfContents()

# ╔═╡ b05ae27f-13df-4ae6-8189-1c21d86ee2df
html"""<style>
main {
    max-width: 1030px;
    align-self: flex-start;
    margin-left: 20px;
}
"""

# ╔═╡ Cell order:
# ╟─87597e02-19c7-418a-8eb3-ad25695bb5bf
# ╟─ba3a32ee-9664-4279-b27b-0fed992e4d1c
# ╟─f64f9468-25ee-4799-b776-b1cfbe5b7f57
# ╠═c259acbf-da3d-4b73-af2b-8fc22f2548d7
# ╠═15741d2d-5fc2-4dc0-baf1-3ba3a72b5399
# ╠═d6219a57-7c1d-4b23-9b80-ae96f005b762
# ╟─60559545-8f85-4027-8471-4ce3d7cc952a
# ╠═3c4b1028-bcf7-45b0-bba2-dab5a63724d0
# ╠═785a2042-bb12-4a2a-8117-c2cfb396d29c
# ╠═433c64f5-3f0c-4442-a5ea-5d57e14a2038
# ╠═4ae8bce1-5890-4647-94e5-e0dfd14e9eb5
# ╠═188cf9ab-f839-490b-99fa-04611c355bfe
# ╟─0e1f6708-53d8-4cca-acf3-e7b2cac1dd7a
# ╟─de7e3ebb-f9d2-4830-9bec-f034d965ab76
# ╟─b12592bf-ec54-435f-b795-183cf84bcb4c
# ╠═91d634c4-e4ec-49b8-a035-5df0740f9fb0
# ╠═e00fc7cb-05a7-4215-b7f6-5bad647ab090
# ╟─8cb9eca6-2ebe-4faa-b6c6-5ea267e64dc9
# ╠═4c62a111-112c-4b8d-979d-dea84e35232a
# ╟─2d56d9c0-b23a-4399-a86e-77edbf25d8b5
# ╟─1174095b-a339-47b2-9f21-3f5b2c52b4eb
# ╠═343f3226-a6e1-4bf5-a816-f92cb74cf434
# ╠═e9f30692-375c-49fc-bf0d-f045951ae11f
# ╟─21d22e3c-a56a-4c03-a9aa-226acea40c72
# ╟─8bb918eb-cdb6-4a15-a344-88ec22df92b2
# ╟─592b8b00-91a0-4ab4-b26b-9fc358978e9a
# ╟─701ebc68-2809-4544-be85-c4fadae6c391
# ╠═c863a4fc-270b-4473-beaa-d34ba2d719de
# ╠═5bc2c3e1-62dc-4f64-9339-87c6cbca8a8b
# ╠═a79fb3f3-a07c-43b8-9be6-72e7347ff40e
# ╟─ebf32f98-3ba0-4650-83ca-4d7a5b310c05
# ╟─fb17d03b-880c-4be6-bd50-f31bdd3e2fd4
# ╟─bbf59177-36a8-480c-af39-1998b37032d2
# ╟─f7db71f7-9bac-4aad-a636-5c60ae4c259b
# ╟─9dc77012-f845-4a5b-a059-fd13c1308b3d
# ╟─8af9d40f-8d50-4376-82e4-01981a54893f
# ╟─92175df7-b060-4f6f-8acb-265bdbe6713a
# ╟─0d9dd1ea-9ef8-460c-8d95-20a46eb1ebb6
# ╟─f8e4b341-1613-4993-a193-437f2adf4a21
# ╟─cfce959c-c171-4923-b5f9-29666d9c20e4
# ╟─fe0b7b47-39bb-45f8-873e-9a4a07d7d1b3
# ╟─b6da427f-0c72-4de7-bdca-b85688ac9bed
# ╠═b94f222a-0c35-49d7-94f2-a2347efb6b2d
# ╟─e7541e93-8b7a-4c26-901e-d887931b8103
# ╠═e858ca22-5ddd-45aa-8815-c3acc3f9741e
# ╟─1c42a7ce-3e97-48b3-8ddd-89cab6177bfe
# ╟─0b07ac3c-d61e-43ff-acff-f907019835d6
# ╟─473302a3-2541-49a1-8f21-906bc2c9d279
# ╟─34397405-e939-464c-a384-5d6c8fb9ca83
# ╟─9b81f986-ea3c-49f5-be04-b12b392aacd6
# ╟─7b2e6266-6dfe-47b5-8f28-3b9d76f57a42
# ╟─c8621274-00c9-4f67-b285-c624d1346527
# ╟─7d63e3eb-6132-439a-bc68-05e73cd5b068
# ╟─fca2fa7d-99fe-44bd-a34a-1b1552a3ae48
# ╟─72e90532-02da-4d31-88c3-7766393d9f92
# ╠═c487c08a-93a1-42dc-b465-3268be5954df
# ╟─d4e7987c-94bf-4194-8e16-88c2581832d1
# ╟─4c71c3d6-1f9e-4812-9f36-2e65a43c8207
# ╟─b2b81585-2144-4b37-b0a7-286c862bb812
# ╟─0e9b78b0-208f-4b29-96c3-4b4129416bdd
# ╟─9833f763-ece1-4739-a6d1-0f3182d32d5d
# ╟─9d89d84d-f107-4a8d-be18-55af8edb546e
# ╟─833476e3-f1f1-4042-87e0-4c21ef37f7a7
# ╟─2c88f409-6763-41d8-b1cc-168df855728b
# ╟─38bb2810-3a30-4a39-99dd-83518c093be2
# ╠═c23c4063-947e-436a-910b-c43afdb9c244
# ╠═26caf850-c46e-453c-aad9-c7bdb71f858b
# ╠═bcbd90c4-60b1-4c90-aa4b-69b0634e2fcd
# ╠═0b6ec47b-3610-453f-b0c3-a99d5fb83ab4
# ╠═c1bf8bf8-fb05-4731-8140-c0f6332496bd
# ╠═b957ac39-458f-48b6-8683-229d09ecaea4
# ╟─d9877f46-c461-4615-bb0c-0c7541238899
# ╟─93465fa0-83bc-413b-b46e-7ed0d3d0d1a7
# ╟─022f7788-a2d0-4489-af0b-71692f800fdf
# ╟─34bda3ea-4238-4670-9f6b-9207f63047cb
# ╠═72522f83-2b3c-4e40-a1de-c3508d50fdab
# ╠═6fb4b458-6192-4000-8efd-f6fb02488861
# ╟─31f212b1-dac3-4c77-a340-dd2837fc8807
# ╟─4441d217-5f2e-4052-8fda-e11672ae1813
# ╟─42b4e6fe-a897-11eb-2676-f7fd96f35a22
# ╟─e0756a3b-998a-4f8d-9ab3-df50422ce14b
# ╟─44ae02ac-6f9a-4af5-8901-3428da03af22
# ╠═b05ae27f-13df-4ae6-8189-1c21d86ee2df


text:
### Determine the number of flights from each airport, include a list of any airports not used.
reload("nmr")

nmr.runjob(nmr.NMR(1, ["AComp_Passenger_data.csv"], nmr.mapper_parserecordacomp, "acomp.csv"))
nmr.runjob(nmr.NMR(1, ["acomp.csv"], nmr.reducer_numberofflights, "acomp_flights_reduced.csv"))
c = nmr.runcombiner("acomp_flights_reduced.csv", nmr.combiner_parsejson)

@test typeof(c) == Vector{String}
@test length(c) == 17
@test JSON.parse(c[1])[1] == "AMS"
@test JSON.parse(c[1])[2] == 13

ext: jl

max_stars_repo_name: deyandyankov/nmr

max_issues_repo_name: deyandyankov/nmr

max_forks_repo_name: deyandyankov/nmr


text:
# Rename file to offsetarrays?

# If ndim(::OffsetArray) == 1, we can convert to a regular strided pointer and offset.
@inline VectorizationBase.stridedpointer(a::OffsetArrays.OffsetArray{<:Any,1}) = gesp(stridedpointer(parent(a)), (-@inbounds(a.offsets[1]),))

struct OffsetStridedPointer{T, N, P <: VectorizationBase.AbstractStridedPointer{T}} <: VectorizationBase.AbstractStridedPointer{T}
    ptr::P
    offsets::NTuple{N,Int}
end
# if ndim(A::OffsetArray) ≥ 2, then eachindex(A) isa Base.OneTo, index starting at 1.
# but multiple indexing is calculated using offsets, so we need a special type to express this.
@inline function VectorizationBase.stridedpointer(A::OffsetArrays.OffsetArray)
    OffsetStridedPointer(stridedpointer(parent(A)), A.offsets)
end

@inline function VectorizationBase.stridedpointer(
    B::Adjoint{T,A}
) where {T,A<:OffsetArrays.OffsetArray{T}}
    Boff = parent(B)
    OffsetStridedPointer(
        stridedpointer(parent(Boff)'),
        Boff.offsets
    )
end
@inline function Base.transpose(A::OffsetStridedPointer)
    OffsetStridedPointer(
        transpose(A.ptr), A.offsets
    )
end
# Tuple of length == 1, use ind directly.
# @inline VectorizationBase.offset(ptr::OffsetStridedPointer, ind::Tuple{I}) where {I} = VectorizationBase.offset(ptr.ptr, ind)
# Tuple of length > 1, subtract offsets.
# @inline VectorizationBase.offset(ptr::OffsetStridedPointer{<:Any,N}, ind::Tuple) where {N} = VectorizationBase.offset(ptr.ptr, ntuple(n -> ind[n] + ptr.offsets[n], Val{N}()))
@inline VectorizationBase.stride1offset(ptr::OffsetStridedPointer, ind::Tuple{I}) where {I} = VectorizationBase.stride1offset(ptr.ptr, ind)
@inline VectorizationBase.stridedoffset(ptr::OffsetStridedPointer, ind::Tuple{I}) where {I} = VectorizationBase.stridedoffset(ptr.ptr, ind)
# Tuple of length > 1, subtract offsets.
@inline VectorizationBase.stride1offset(ptr::OffsetStridedPointer{<:Any,N}, ind::Tuple) where {N} = VectorizationBase.stride1offset(ptr.ptr, ntuple(n -> vsub(ind[n], ptr.offsets[n]), Val{N}()))
@inline VectorizationBase.stridedoffset(ptr::OffsetStridedPointer{<:Any,N}, ind::Tuple) where {N} = VectorizationBase.stridedoffset(ptr.ptr, ntuple(n -> vsub(ind[n], ptr.offsets[n]), Val{N}()))
@inline Base.similar(p::OffsetStridedPointer, ptr::Ptr) = OffsetStridedPointer(similar(p.ptr, ptr), p.offsets)
@inline Base.pointer(p::OffsetStridedPointer) = pointer(p.ptr)
@inline VectorizationBase.gesp(p::OffsetStridedPointer, i) = similar(p.ptr, gep(p, i))
# @inline VectorizationBase.gesp(p::OffsetStridedPointer, i) = similar(p, gep(p.ptr, i))
# If an OffsetArray is getting indexed by a (loop-)constant value, then this particular vptr object cannot also be eachindexed, so we can safely return a stridedpointer
@inline function VectorizationBase.subsetview(ptr::OffsetStridedPointer{<:Any,N}, ::Val{I}, i) where {I,N}
    subsetview(gesp(ptr.ptr, ntuple(n -> 0 - @inbounds(ptr.offsets[n]), Val{N}())), Val{I}(), i)
end

@inline function VectorizationBase.stride1offset(ptr::OffsetStridedPointer{<:Any,<:Any,<:VectorizationBase.AbstractBitPointer}, ind::Tuple{I}) where {I}
    VectorizationBase.stride1offset(ptr.ptr, (vsub(ind[1], ptr.offsets[1]),))
end
@inline function VectorizationBase.stridedoffset(ptr::OffsetStridedPointer{<:Any,<:Any,<:VectorizationBase.AbstractBitPointer}, ind::Tuple{I}) where {I}
    VectorizationBase.stridedoffset(ptr.ptr, (vsub(ind[1], ptr.offsets[1]),))
end
@inline VectorizationBase.gesp(ptr::VectorizationBase.AbstractBitPointer, i) = OffsetStridedPointer(ptr, vsub.(-1, unwrap.(i)))


text:
# extract items useful for plotting from the history generated by the
# fallback `result` method:
function plotting_report(fields, scales, history)

    n_models = length(history)
    n_parameters = length(fields)

    A = Array{Any}(undef, (n_models, n_parameters))
    measurements = Vector{Float64}(undef, n_models)

    for j in eachindex(history)
        entry = history[j]
        A[j,:] = [recursive_getproperty(entry.model, fld) for fld in fields]
        measurements[j] = entry.measurement[1]
    end

    return plotting=(parameter_names=string.(fields) |> collect,
                     parameter_scales=scales |> collect,
                     parameter_values = A,
                     measurements = measurements)

end

# return a named tuple with some keys removed:
function delete(nt::NamedTuple, target_keys...)
    zipped = tuple(zip(keys(nt), values(nt))...)
    filtered = filter(zipped) do tup
        !(tup[1] in target_keys)
    end
    return (; filtered...)
end

signature(measure) =
    if orientation(measure) == :loss
        1
    elseif orientation(measure) == :score
        -1
    else
        0
    end

max_stars_repo_name: alan-turing-institute/MLJTuning

max_issues_repo_name: alan-turing-institute/MLJTuning

max_forks_repo_name: alan-turing-institute/MLJTuning


text:
@testset "Core definitions and constructors" begin
    @testset "Inner Constructors" begin
        # SVector
        @test SVector{1,Int}((1,)).data === (1,)
        @test SVector{1,Float64}((1,)).data === (1.0,)
        @test_throws Exception SVector{1,Int}()

        @test SMatrix{1,1,Int}((1,)).data === (1,)
        @test SMatrix{1,1,Float64}((1,)).data === (1.0,)
        @test_throws Exception SMatrix{1,Int}()

        #@test MArray{(1,),Int,1,Tuple{Int}}((1,)).data === (1,)
        #@test MArray{(1,),Float64,1,Tuple{Float64}}((1,)).data === (1.0,)

        # default constructors
        #@test try; MArray{(1,),Int,1,Tuple{Int}}(); true; catch; false; end
    end

    @testset "Type parameter errors" begin
        # (not sure what type of exception these should be?)
        @test_throws Exception SVector{1.0,Int}((1,))
        @test_throws Exception SVector{2,Int}((1,))
        @test_throws Exception SVector{1,3}((1,))

        @test_throws Exception SMatrix{1.0,1,Int,1}((1,))
        @test_throws Exception SMatrix{1,1.0,Int,1}((1,))
        @test_throws Exception SMatrix{2,1,Int,1}((1,))
        @test_throws Exception SMatrix{1,2,Int,1}((1,))
        @test_throws Exception SMatrix{1,1,3,1}((1,))
        @test_throws Exception SMatrix{1,1,Int,2}((1,))

    end
#=
    @testset "Outer constructors" begin
        @test SArray{(1,),Int,1}((1,)).data === (1,)
        @test_inferred SArray{(1,),Int,1}((1,))
        @test SArray{(1,),Int}((1,)).data === (1,)
        @test_inferred SArray{(1,),Int}((1,))
        @test SArray{(1,)}((1,)).data === (1,)
        @test_inferred SArray{(1,)}((1,))

        @test MArray{(1,),Int,1}((1,)).data === (1,)
        @test_inferred MArray{(1,),Int,1}((1,))
        @test_inferred MArray{(1,),Int,1}()
        @test MArray{(1,),Int}((1,)).data === (1,)
        @test_inferred MArray{(1,),Int}((1,))
        @test_inferred MArray{(1,),Int}()
        @test MArray{(1,)}((1,)).data === (1,)
        @test_inferred MArray{(1,)}((1,))
    end

    @testset "Constructors that convert types" begin
        @test SArray{(1,),Float64,1}((1,)).data === (1.0,)
        @test_inferred SArray{(1,),Float64,1}((1,))
        @test SArray{(1,),Float64}((1,)).data === (1.0,)
        @test_inferred SArray{(1,),Float64}((1,))

        @test MArray{(1,),Float64,1}((1,)).data === (1.0,)
        @test_inferred MArray{(1,),Float64,1}((1,))
        @test MArray{(1,),Float64}((1,)).data === (1.0,)
        @test_inferred MArray{(1,),Float64}((1,))
    end

    @testset "eltype conversion" begin
        sa_int = SArray{(1,)}((1,))
        ma_int = MArray{(1,)}((1,))

        sa_float = SArray{(1,)}((1.0,))
        ma_float = MArray{(1,)}((1.0,))

        @test convert(SArray{(1,),Float64}, sa_int) === sa_float
        @test_inferred convert(SArray{(1,),Float64}, sa_int)
        @test convert(SArray{(1,),Float64,1}, sa_int) === sa_float
        @test_inferred convert(SArray{(1,),Float64,1}, sa_int)
        @test convert(SArray{(1,),Float64,1,Tuple{Float64}}, sa_int) === sa_float
        @test_inferred convert(SArray{(1,),Float64,1,Tuple{Float64}}, sa_int)

        @test convert(MArray{(1,),Float64}, ma_int) == ma_float
        @test_inferred convert(MArray{(1,),Float64}, ma_int)
        @test convert(MArray{(1,),Float64,1}, ma_int) == ma_float
        @test_inferred convert(MArray{(1,),Float64,1}, ma_int)
        @test convert(MArray{(1,),Float64,1,Tuple{Float64}}, ma_int) == ma_float
        @test_inferred convert(MArray{(1,),Float64,1,Tuple{Float64}}, ma_int)
    end

    @testset "StaticArray conversion" begin
        sa_int = SArray{(1,)}((1,))
        ma_int = MArray{(1,)}((1,))

        sa_float = SArray{(1,)}((1.0,))
        ma_float = MArray{(1,)}((1.0,))

        # SArray -> MArray
        @test convert(MArray, sa_int) == ma_int
        @test_inferred convert(MArray, sa_int)
        @test convert(MArray{(1,)}, sa_int) == ma_int
        @test_inferred convert(MArray{(1,)}, sa_int)
        @test convert(MArray{(1,),Int}, sa_int) == ma_int
        @test_inferred convert(MArray{(1,),Int}, sa_int)
        @test convert(MArray{(1,),Int,1}, sa_int) == ma_int
        @test_inferred convert(MArray{(1,),Int,1}, sa_int)
        @test convert(MArray{(1,),Int,1,Tuple{Int}}, sa_int) == ma_int
        @test_inferred convert(MArray{(1,),Int,1,Tuple{Int}}, sa_int)

        @test convert(MArray{(1,),Float64}, sa_int) == ma_float
        @test_inferred convert(MArray{(1,),Float64}, sa_int)
        @test convert(MArray{(1,),Float64,1}, sa_int) == ma_float
        @test_inferred convert(MArray{(1,),Float64,1}, sa_int)
        @test convert(MArray{(1,),Float64,1,Tuple{Float64}}, sa_int) == ma_float
        @test_inferred convert(MArray{(1,),Float64,1,Tuple{Float64}}, sa_int)

        # MArray -> SArray
        @test convert(SArray, ma_int) === sa_int
        @test_inferred convert(SArray, ma_int)
        @test convert(SArray{(1,)}, ma_int) === sa_int
        @test_inferred convert(SArray{(1,)}, ma_int)
        @test convert(SArray{(1,),Int}, ma_int) === sa_int
        @test_inferred convert(SArray{(1,),Int}, ma_int)
        @test convert(SArray{(1,),Int,1}, ma_int) === sa_int
        @test_inferred convert(SArray{(1,),Int,1}, ma_int)
        @test convert(SArray{(1,),Int,1,Tuple{Int}}, ma_int) === sa_int
        @test_inferred convert(SArray{(1,),Int,1,Tuple{Int}}, ma_int)

        @test convert(SArray{(1,),Float64}, ma_int) === sa_float
        @test_inferred convert(SArray{(1,),Float64}, ma_int)
        @test convert(SArray{(1,),Float64,1}, ma_int) === sa_float
        @test_inferred convert(SArray{(1,),Float64,1}, ma_int)
        @test convert(SArray{(1,),Float64,1,Tuple{Float64}}, ma_int) === sa_float
        @test_inferred convert(SArray{(1,),Float64,1,Tuple{Float64}}, ma_int)
    end

    @testset "AbstractArray conversion" begin
        sa = SArray{(2,2)}((3,4,5,6))
        ma = MArray{(2,2)}((3,4,5,6))
        a = [3 5; 4 6]

        @test convert(SArray{(2,2)}, a) === sa
        @test_inferred convert(SArray{(2,2)}, a)
        @test convert(SArray{(2,2),Int}, a) === sa
        @test_inferred convert(SArray{(2,2),Int}, a)
        @test convert(SArray{(2,2),Int,2}, a) === sa
        @test_inferred convert(SArray{(2,2),Int,2}, a)
        @test convert(SArray{(2,2),Int,2,NTuple{4,Int}}, a) === sa
        @test_inferred convert(SArray{(2,2),Int,2,NTuple{4,Int}}, a)

        @test convert(MArray{(2,2)}, a) == ma
        @test_inferred convert(MArray{(2,2)}, a)
        @test convert(MArray{(2,2),Int}, a) == ma
        @test_inferred convert(MArray{(2,2),Int}, a)
        @test convert(MArray{(2,2),Int,2}, a) == ma
        @test_inferred convert(MArray{(2,2),Int,2}, a)
        @test convert(MArray{(2,2),Int,2,NTuple{4,Int}}, a) == ma
        @test_inferred convert(MArray{(2,2),Int,2,NTuple{4,Int}}, a)

        @test convert(Array, sa) == a
        @test_inferred convert(Array, sa)
        @test convert(Array{Int}, sa) == a
        @test_inferred convert(Array{Int}, sa)
        @test convert(Array{Int,2}, sa) == a
        @test_inferred convert(Array{Int,2}, sa)

        @test convert(Array, ma) == a
        @test_inferred convert(Array, ma)
        @test convert(Array{Int}, ma) == a
        @test_inferred convert(Array{Int}, ma)
        @test convert(Array{Int,2}, ma) == a
        @test_inferred convert(Array{Int,2}, ma)
    end =#
end


text: include("./sop.jl")


text:
module Histograms
import StatsBase

export Hist2d

type Hist1d
    P::Integer
    C::Integer
    N::Integer
    wordlist::Array{Int32,1}
    wordcnt::Array{Float64,1}
    entropy::Float64
end

type Hist2d
    wordcnt::Array{Float64,2}
end

function hist1d{T<:Integer}(X::Array{T,1})
    ntrials = length(X)
    N = float(StatsBase.counts(X, 0:maximum(X)))
    W = int32(sort(unique(X)))
    nW = length(W)
    N = N[N.>0]
    return Hist1d(ntrials,nW,1,W,N,0.)
end

function hist1d{T<:Integer}(X::Array{T,2})
    ntrials = size(X,2)
    mx = maximum(X,2)+1
    M = cumprod(mx)
    M = [1;M[1:end-1]]
    Y = X'*M
    return hist1d(Y)
end

function hist1d{T<:Integer}(X::Array{T,3})
    nbins = size(X,3)
    H = Array(Hist1d,nbins)
    for i=1:nbins
        H[i] = hist1d(squeeze(X[:,:,i],3))
    end
    return H
end

end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(highlight_search),Vector{SubString{String}},Regex})   # time: 0.059157036
    Base.precompile(Tuple{Core.kwftype(typeof(textview)),NamedTuple{(:active_match, :search_regex), Tuple{Int64, Regex}},typeof(textview),String,NTuple{4, Int64}})   # time: 0.058709044
    Base.precompile(Tuple{Core.kwftype(typeof(textview)),NamedTuple{(:active_match, :frozen_lines_at_beginning, :search_regex), Tuple{Int64, Int64, Regex}},typeof(textview),String,NTuple{4, Int64}})   # time: 0.025712194
    Base.precompile(Tuple{Core.kwftype(typeof(textview)),NamedTuple{(:active_match, :frozen_lines_at_beginning, :frozen_columns_at_beginning, :search_regex), Tuple{Int64, Int64, Int64, Regex}},typeof(textview),String,NTuple{4, Int64}})   # time: 0.01779813
    Base.precompile(Tuple{typeof(textview),String,NTuple{4, Int64}})   # time: 0.015094872
    Base.precompile(Tuple{typeof(parse_decoration),String})   # time: 0.014008783
    Base.precompile(Tuple{Core.kwftype(typeof(textview)),NamedTuple{(:active_match, :frozen_lines_at_beginning, :parse_decorations_before_view, :search_regex), Tuple{Int64, Int64, Bool, Regex}},typeof(textview),String,NTuple{4, Int64}})   # time: 0.013855798
    Base.precompile(Tuple{typeof(align_string_per_line),String,Int64,Symbol})   # time: 0.012167775
    Base.precompile(Tuple{Core.kwftype(typeof(textview)),NamedTuple{(:active_match, :frozen_lines_at_beginning, :frozen_columns_at_beginning, :search_regex, :show_ruler), Tuple{Int64, Int64, Int64, Regex, Bool}},typeof(textview),String,NTuple{4, Int64}})   # time: 0.011858754
    Base.precompile(Tuple{typeof(highlight_search),String,Regex})   # time: 0.011790905
    Base.precompile(Tuple{typeof(fit_string_in_field),String,Int64})   # time: 0.010749603
    Base.precompile(Tuple{Core.kwftype(typeof(textview)),NamedTuple{(:frozen_columns_at_beginning, :frozen_lines_at_beginning, :maximum_number_of_columns, :maximum_number_of_lines), NTuple{4, Int64}},typeof(textview),String,NTuple{4, Int64}})   # time: 0.010748933
    Base.precompile(Tuple{Core.kwftype(typeof(textview)),NamedTuple{(:maximum_number_of_lines, :maximum_number_of_columns), Tuple{Int64, Int64}},typeof(textview),String,NTuple{4, Int64}})   # time: 0.009999969
    Base.precompile(Tuple{Core.kwftype(typeof(textview)),NamedTuple{(:frozen_columns_at_beginning, :maximum_number_of_columns), Tuple{Int64, Int64}},typeof(textview),String,NTuple{4, Int64}})   # time: 0.009609414
    Base.precompile(Tuple{typeof(align_string),String,Int64,Symbol})   # time: 0.006738714
    Base.precompile(Tuple{Core.kwftype(typeof(textview)),NamedTuple{(:frozen_columns_at_beginning, :frozen_lines_at_beginning, :title_lines), Tuple{Int64, Int64, Int64}},typeof(textview),Vector{SubString{String}},NTuple{4, Int64}})   # time: 0.00565889
    Base.precompile(Tuple{typeof(get_and_remove_decorations),String})   # time: 0.004928707
    Base.precompile(Tuple{Core.kwftype(typeof(highlight_search)),NamedTuple{(:active_match,), Tuple{Int64}},typeof(highlight_search),String,Regex})   # time: 0.004911768
    Base.precompile(Tuple{Core.kwftype(typeof(textview)),NamedTuple{(:active_match, :search_matches), Tuple{Int64, Dict{Int64, Vector{Tuple{Int64, Int64}}}}},typeof(textview),Vector{SubString{String}},NTuple{4, Int64}})   # time: 0.004447444
    Base.precompile(Tuple{Core.kwftype(typeof(highlight_search)),NamedTuple{(:active_match, :start_line, :end_line), Tuple{Int64, Int64, Int64}},typeof(highlight_search),Vector{SubString{String}},Regex})   # time: 0.004381782
    Base.precompile(Tuple{typeof(get_decorations),String})   # time: 0.00303308
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:foreground, :background, :underline), Tuple{String, String, DecorationState}},Type{Decoration}})   # time: 0.002692545
    Base.precompile(Tuple{Core.kwftype(typeof(fit_string_in_field)),NamedTuple{(:crop_side, :printable_string_width), Tuple{Symbol, Int64}},typeof(fit_string_in_field),String,Int64})   # time: 0.002690417
    Base.precompile(Tuple{Core.kwftype(typeof(highlight_search)),NamedTuple{(:active_match,), Tuple{Int64}},typeof(highlight_search),Vector{SubString{String}},Regex})   # time: 0.002600482
    Base.precompile(Tuple{Core.kwftype(typeof(fit_string_in_field)),NamedTuple{(:add_space_in_continuation_char,), Tuple{Bool}},typeof(fit_string_in_field),String,Int64})   # time: 0.002591723
    Base.precompile(Tuple{Core.kwftype(typeof(fit_string_in_field)),NamedTuple{(:printable_string_width,), Tuple{Int64}},typeof(fit_string_in_field),String,Int64})   # time: 0.00258039
    Base.precompile(Tuple{typeof(left_crop),String,Int64})   # time: 0.00254965
    Base.precompile(Tuple{Core.kwftype(typeof(fit_string_in_field)),NamedTuple{(:crop_side,), Tuple{Symbol}},typeof(fit_string_in_field),String,Int64})   # time: 0.002482469
    Base.precompile(Tuple{Core.kwftype(typeof(fit_string_in_field)),NamedTuple{(:add_space_in_continuation_char, :keep_ansi), Tuple{Bool, Bool}},typeof(fit_string_in_field),String,Int64})   # time: 0.002358202
    Base.precompile(Tuple{typeof(printable_textwidth_per_line),String})   # time: 0.002345092
    Base.precompile(Tuple{Core.kwftype(typeof(align_string)),NamedTuple{(:fill,), Tuple{Bool}},typeof(align_string),String,Int64,Symbol})   # time: 0.00233521
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:foreground, :background, :reversed), Tuple{String, String, DecorationState}},Type{Decoration}})   # time: 0.002195955
    Base.precompile(Tuple{Core.kwftype(typeof(fit_string_in_field)),NamedTuple{(:keep_ansi,), Tuple{Bool}},typeof(fit_string_in_field),String,Int64})   # time: 0.002191823
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:foreground, :background, :bold), Tuple{String, String, DecorationState}},Type{Decoration}})   # time: 0.002106865
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:foreground, :background, :bold, :underline, :reversed), Tuple{String, String, DecorationState, DecorationState, DecorationState}},Type{Decoration}})   # time: 0.002102631
    Base.precompile(Tuple{Core.kwftype(typeof(fit_string_in_field)),NamedTuple{(:add_space_in_continuation_char, :crop_side), Tuple{Bool, Symbol}},typeof(fit_string_in_field),String,Int64})   # time: 0.002086886
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:background,), Tuple{String}},Type{Decoration}})   # time: 0.001959384
    Base.precompile(Tuple{typeof(convert),Type{String},Decoration})   # time: 0.001886508
    Base.precompile(Tuple{Core.kwftype(typeof(right_crop)),NamedTuple{(:printable_string_width,), Tuple{Int64}},typeof(right_crop),String,Int64})   # time: 0.001595281
    Base.precompile(Tuple{typeof(drop_inactive_properties),Decoration})   # time: 0.001566229
end


text:
using ADI
using ProgressLogging
using PSFModels: Gaussian
using PyCall
using SiriusB

fits = pyimport("astropy.io.fits")
pro = pyimport("proplot")

pro.use_style("ggplot")
pro.rc["image.cmap"] = "inferno"
pro.rc["image.origin"] = "lower"
pro.rc["grid"] = false

_label(a::Type{<:PCA}) = "PCA"
_label(a::Type{<:NMF}) = "NMF"
_label(a::Type{<:GreeDS}) = "GreeDS"

_epochs = Dict(
    "2020feb04" => "2020-02-04",
    "2020nov21" => "2020-11-21",
    "2020nov28" => "2020-11-28"
)

function plot_mosaic(res_cubes, angles; epoch, label, ncomps)
    py"""
    import proplot as pro
    fig, axs = pro.subplots(nrows=2, ncols=5, share=3, width="7.5in", space=0)
    """
    i = 0
    flats = collapse.(res_cubes, Ref(angles))
    vmax = mapreduce(maximum, max, flats)
    vmin = mapreduce(minimum, min, flats)

    for (n, flat) in zip(ncomps, flats)
        lab = "ncomp=$n"

        py"""
        m = axs[$i].imshow($flat, vmin=$vmin, vmax=$vmax)
        axs[$i].text(6, 6, $lab, color="w")
        """
        i += 1
    end

    ticks = range(40, 159, length=3)
    tick_labs = @. string(round(SiriusB.pxscale * (ticks - 99.5), digits=1))

    figname = make_filename_friendly(figuredir("reports", "$(epoch)_$(label)_mosaic.pdf"))
    py"""
    axs.format(
        yticks=$ticks,
        yticklabels=$tick_labs,
        ylabel="y [arcsec]",
        xticks=$ticks,
        xticklabels=$tick_labs,
        xlabel="x [arcsec]",
    )
    fig.colorbar(m, loc="r")
    fig.savefig($figname)
    """
end

function plot_results(res_cubes, angles, contrast_curves; epoch, label, fwhm)
    py"""
    import proplot as pro
    fig, axs = pro.subplots([[1, 3, 3], [2, 3, 3]], share=0, width="7.5in")
    cycle = pro.Cycle("viridis_r", $(length(res_cubes)))
    lines = []
    """
    # plot the average STIM map and SLIM probability map
    ticks = range(20, 179, length=5)
    tick_labs = @. string(round(SiriusB.pxscale * (ticks - 99.5), digits=1))

    stimmap, slimmask = slimmap(res_cubes, angles; N=ceil(Int, π/4 * fwhm^2))
    slimprob = stimmap .* slimmask
    py"""
    axs[0].imshow($stimmap, vmin=0, vmax=1)
    axs[0].text(6, 6, "mean STIM", color="w")
    m = axs[1].imshow($slimprob, vmin=0, vmax=1)
    axs[1].text(6, 6, "SLIM map", color="w")
    fig.colorbar(m, loc="l", label="STIM probability")

    axs[0].format(
        ylabel="y [arcsec]",
        yticks=$ticks,
        yticklabels=$tick_labs,
        xticks="null"
    )
    axs[1].format(
        ylabel="y [arcsec]",
        xlabel="x [arcsec]",
        yticks=$ticks,
        yticklabels=$tick_labs,
        xticks=$ticks,
        xticklabels=$tick_labs
    )
    """
    # now plot contrast curves
    for curve in contrast_curves
        dist = curve.distance .* SiriusB.auscale
        py"""
        l = axs[2].plot($dist, $(curve.contrast), cycle=cycle)
        lines.extend(l)
        """
    end
    figname = make_filename_friendly(figuredir("reports", "$(epoch)_$(label)_contrast_curves.pdf"))
    py"""
    axs[2].text(0.25, 5e-4, $(_epochs[epoch]), fontsize=14)
    fig.colorbar(lines, loc="r", label="ncomp", values=range(1, 11))
    axs[2].dualx(lambda x: x * $(SiriusB.parallax), label="separation [arcsec]")
    ax2 = axs[2].alty(
        yticks=[0.09, 0.31, 0.54, 0.77, 1],
        yticklabels=["8", "6", "4", "2", "0"],
        label="Δ mag"
    )

    axs[2].format(
        ylim=(10**(-3.5), 1),
        yscale="log",
        yformatter="log",
        xlabel="projected separation [AU]",
        ylabel="5σ contrast",
        grid=True
    )
    fig.savefig($figname)
    """
end

function make_filename_friendly(filename)
    tmp = lowercase(filename)
    tmp = replace(tmp, " " => "_")
    tmp = replace(tmp, "(" => "-")
    tmp = replace(tmp,  ")" => "")
end

@progress name="epoch" for epoch in ["2020feb04", "2020nov21", "2020nov28"]
    cube = fits.getdata(datadir("epoch_$epoch", "processed", "$(epoch)_sirius-b_cube_calib_registered_crop.fits"))
    angles = fits.getdata(datadir("epoch_$epoch", "processed", "$(epoch)_sirius-b_pa.fits"))

    raw_median = collapse(cube)
    res = gaussian_fit(raw_median, (100.5, 100.5))
    fwhm = res[3]
    psf_model = Gaussian{eltype(cube)}(res[1:2], fwhm)
    @info "FWHM for epoch $epoch is $fwhm px"
    starphot = Metrics.estimate_starphot(cube, fwhm)
    psfphot = Metrics.estimate_starphot(psf_model, fwhm) * exp(res[4])
    @info "star photometry is ≈ $starphot (PSF photometry was $psfphot)"

    # naturally mask out inner FWHM
    av_cube = AnnulusView(cube; inner=fwhm)

    ncomps = 1:10
    @progress name="algorithm" for alg in [PCA, NMF, GreeDS]
        label = _label(alg)
        @withprogress name="ncomp" begin
            i = 1
            N = length(ncomps)
            outputs = map(ncomps) do n
                _alg = alg(n)
                fname = make_filename_friendly(datadir("epoch_$epoch", "residuals", "$(epoch)_sirius-b_residual_$label-$n.fits"))
                res_cube = load_or_produce(fname) do
                    subtract(_alg, av_cube; angles, fwhm)
                end
                fname_cc = make_filename_friendly(datadir("epoch_$epoch", "residuals", "$(epoch)_sirius-b_contrast-curve_$label-$n.csv"))
                cc = load_or_produce(fname_cc) do
                    contrast_curve(_alg, av_cube, angles, psf_model; fwhm, angles, starphot, nbranch=6)
                end
                @logprogress i/N
                i += 1
                res_cube, cc
            end
            res_cubes = map(first, outputs)
            contrast_curves = map(last, outputs)
            plot_mosaic(res_cubes, angles; epoch, label, ncomps)
            plot_results(res_cubes, angles, contrast_curves; epoch, label, fwhm)
        end
    end
end

max_stars_repo_name: mileslucas/sirius-b

max_issues_repo_name: mileslucas/sirius-b

max_forks_repo_name: mileslucas/sirius-b


text:
using JSON

"""
Reads data in from JSON function and treat it like a read only dictionary
"""
type SettingsJSONFile <: SettingsSource
    data::Dict
    filename::AbstractString
    watched::Bool
    file_watcher::Task

    function SettingsJSONFile(filename::AbstractString, watched::Bool=true)
        data = JSON.parsefile(filename)

        # Now get our settings source object so we can update the data
        # asynchroniously if desired
        settings_source = new(data, filename, watched, @async begin end)

        # Now set up the real watcher
        if watched
            settings_source.file_watcher = @async begin
                while settings_source.watched
                    name, status = watch_file(settings_source.filename)
                    if status.changed
                        settings_source.data = JSON.parsefile(settings_source.filename)
                    end
                end
            end
        end

        return settings_source
    end
end

Base.getindex(item::SettingsJSONFile, key) = item.data[key]

Base.keys(item::SettingsJSONFile) = keys(item.data)


text:
module Shapefile
    import Base: read, show, +, /, ./

    export Handle, RGBGradient, LogLinearRGBGradient, LinearRGBGradient

    type Rect{T}
        top::T
        left::T
        bottom::T
        right::T
    end

    abstract ESRIShape

    type NullShape <: ESRIShape
    end

    type Interval{T}
        left::T
        right::T
    end

    type Point{T} <: ESRIShape
        x::T
        y::T
    end

    +{T}(a::Point{T},b::Point{T}) = Point{T}(a.x+b.x,a.y+b.y)
    /{T}(a::Point{T},val::Real) = Point{T}(a.x/val,a.y/val)
    ./{T}(a::Point{T},val::Real) = Point{T}(a.x./val,a.y./val)

    type PointM{T,M} <: ESRIShape
        x::T
        y::T
        m::M # measure
    end

    type PointZ{T,M} <: ESRIShape
        x::T
        y::T
        z::T
        m::M # measure
    end

    type Polyline{T} <: ESRIShape
        MBR::Rect{T}
        parts::Vector{Int32}
        points::Vector{Point{T}}
    end

    type PolylineM{T,M} <: ESRIShape
        MBR::Rect{T}
        parts::Vector{Int32}
        points::Vector{Point{T}}
        measures::Vector{M}
    end

    type PolylineZ{T,M} <: ESRIShape
        MBR::Rect{T}
        parts::Vector{Int32}
        points::Vector{Point{T}}
        zvalues::Vector{T}
        measures::Vector{M}
    end

    type Polygon{T} <: ESRIShape
        MBR::Rect{T}
        parts::Vector{Int32}
        points::Vector{Point{T}}
    end

    show{T}(io::IO,p::Polygon{T}) = print(io,"Polygon(",length(p.points)," ",T," Points)")

    type PolygonM{T,M} <: ESRIShape
        MBR::Rect{T}
        parts::Vector{Int32}
        points::Vector{Point{T}}
        measures::Vector{M}
    end

    type PolygonZ{T,M} <: ESRIShape
        MBR::Rect{T}
        parts::Vector{Int32}
        points::Vector{Point{T}}
        zvalues::Vector{T}
        measures::Vector{M}
    end

    type MultiPoint{T} <: ESRIShape
        MBR::Rect{T}
        points::Vector{Point{T}}
    end

    type MultiPointM{T,M} <: ESRIShape
        MBR::Rect{T}
        points::Vector{Point{T}}
        measures::Vector{M}
    end

    type MultiPointZ{T,M} <: ESRIShape
        MBR::Rect{T}
        points::Vector{Point{T}}
        zvalues::Vector{T}
        measures::Vector{M}
    end

    type MultiPatch{T,M} <: ESRIShape
        MBR::Rect{T}
        parts::Vector{Int32}
        parttypes::Vector{Int32}
        points::Vector{Point{T}}
        zvalues::Vector{T}
        # measures::Vector{M} # (optional)
    end

    type Handle
        code::Int32
        length::Int32
        version::Int32
        shapeType::Int32
        MBR::Rect{Float64}
        zrange::Interval{Float64}
        mrange::Interval{Float64}
        shapes::Vector{ESRIShape}
    end

    function read{T}(io::IO,::Type{Rect{T}})
        minx = read(io,T)
        miny = read(io,T)
        maxx = read(io,T)
        maxy = read(io,T)
        Rect{T}(miny,minx,maxy,maxx)
    end

    read(io::IO,::Type{NullShape}) = NullShape()

    function read{T}(io::IO,::Type{Point{T}})
        x = read(io,T)
        y = read(io,T)
        Point{T}(x,y)
    end

    function read{T,M}(io::IO,::Type{PointM{T,M}})
        x = read(io,T)
        y = read(io,T)
        m = read(io,M)
        PointM{T,M}(x,y,m)
    end

    function read{T,M}(io::IO,::Type{PointZ{T,M}})
        x = read(io,T)
        y = read(io,T)
        z = read(io,T)
        m = read(io,M)
        PointZ{T,M}(x,y,z,m)
    end

    function read{T}(io::IO,::Type{Polyline{T}})
        box = read(io,Rect{T})
        numparts = read(io,Int32)
        numpoints = read(io,Int32)
        parts = Array(Int32,numparts)
        read!(io, parts)
        points = Array(Point{T},numpoints)
        read!(io, points)
        Polyline{T}(box,parts,points)
    end

    function read{T,M}(io::IO,::Type{PolylineM{T,M}})
        box = read(io,Rect{T})
        numparts = read(io,Int32)
        numpoints = read(io,Int32)
        parts = Array(Int32,numparts)
        read!(io, parts)
        points = Array(Point{T},numpoints)
        read!(io, points)
        mrange = Array(M,2)
        read!(io, mrange)
        measures = Array(M,numpoints)
        read!(io, measures)
        PolylineM{T,M}(box,parts,points,measures)
    end

    function read{T,M}(io::IO,::Type{PolylineZ{T,M}})
        box = read(io,Rect{T})
        numparts = read(io,Int32)
        numpoints = read(io,Int32)
        parts = Array(Int32,numparts)
        read!(io, parts)
        points = Array(Point{T},numpoints)
        read!(io, points)
        zrange = Array(T,2)
        read!(io, zrange)
        zvalues = Array(T,numpoints)
        read!(io, zvalues)
        mrange = Array(M,2)
        read!(io, mrange)
        measures = Array(M,numpoints)
        read!(io, measures)
        PolylineZ{T,M}(box,parts,points,zvalues,measures)
    end

    function read{T}(io::IO,::Type{Polygon{T}})
        box = read(io,Rect{Float64})
        numparts = read(io,Int32)
        numpoints = read(io,Int32)
        parts = Array(Int32,numparts)
        read!(io, parts)
        points = Array(Point{T},numpoints)
        read!(io, points)
        Polygon{T}(box,parts,points)
    end

    function read{T,M}(io::IO,::Type{PolygonM{T,M}})
        box = read(io,Rect{Float64})
        numparts = read(io,Int32)
        numpoints = read(io,Int32)
        parts = Array(Int32,numparts)
        read!(io, parts)
        points = Array(Point{T},numpoints)
        read!(io, points)
        mrange = Array(M,2)
        read!(io, mrange)
        measures = Array(M,numpoints)
        read!(io, measures)
        PolygonM{T,M}(box,parts,points,measures)
    end

    function read{T,M}(io::IO,::Type{PolygonZ{T,M}})
        box = read(io,Rect{Float64})
        numparts = read(io,Int32)
        numpoints = read(io,Int32)
        parts = Array(Int32,numparts)
        read!(io, parts)
        points = Array(Point{T},numpoints)
        read!(io, points)
        zrange = Array(T,2)
        read!(io, zrange)
        zvalues = Array(T,numpoints)
        read!(io, zvalues)
        mrange = Array(M,2)
        read!(io, mrange)
        measures = Array(M,numpoints)
        read!(io, measures)
        PolygonZ{T,M}(box,parts,points,zvalues,measures)
    end

    function read{T}(io::IO,::Type{MultiPoint{T}})
        box = read(io,Rect{Float64})
        numpoints = read(io,Int32)
        points = Array(Point{T},numpoints)
        read!(io, points)
        MultiPoint{T}(box,points)
    end

    function read{T,M}(io::IO,::Type{MultiPointM{T,M}})
        box = read(io,Rect{Float64})
        numpoints = read(io,Int32)
        points = Array(Point{T},numpoints)
        read!(io, points)
        mrange = Array(M,2)
        read!(io, mrange)
        measures = Array(M,numpoints)
        read!(io, measures)
        MultiPointM{T,M}(box,points,measures)
    end

    function read{T,M}(io::IO,::Type{MultiPointZ{T,M}})
        box = read(io,Rect{Float64})
        numpoints = read(io,Int32)
        points = Array(Point{T},numpoints)
        read!(io, points)
        zrange = Array(T,2)
        read!(io, zrange)
        zvalues = Array(T,numpoints)
        read!(io, zvalues)
        mrange = Array(M,2)
        read!(io, mrange)
        measures = Array(M,numpoints)
        read!(io, measures)
        MultiPointZ{T,M}(box,points,zvalues,measures)
    end

    function read{T,M}(io::IO,::Type{MultiPatch{T,M}})
        box = read(io,Rect{Float64})
        numparts = read(io,Int32)
        numpoints = read(io,Int32)
        parts = Array(Int32,numparts)
        read!(io, parts)
        parttypes = Array(Int32,numparts)
        read!(io, parttypes)
        points = Array(Point{T},numpoints)
        read!(io, points)
        zrange = Array(T,2)
        read!(io, zrange)
        zvalues = Array(T,numpoints)
        read!(io, zvalues)
        # mrange = Array(M,2)
        # read!(io, mrange)
        # measures = Array(M,numpoints)
        # read!(io, measures)
        MultiPatch{T,M}(box,parts,parttypes,points,zvalues) #,measures)
    end

    function read(io::IO,::Type{ESRIShape})
        num = bswap(read(io,Int32))
        rlength = bswap(read(io,Int32))
        shapeType = read(io,Int32)
        if(shapeType == 0)
            return read(io,NullShape)
        elseif(shapeType == 1)
            return read(io,Point{Float64})
        elseif(shapeType == 3)
            return read(io,Polyline{Float64})
        elseif(shapeType == 5)
            return read(io,Polygon{Float64})
        elseif(shapeType == 8)
            return read(io,MultiPoint{Float64})
        elseif(shapeType == 11)
            return read(io,PointZ{Float64,Float64})
        elseif(shapeType == 13)
            return read(io,PolylineZ{Float64,Float64})
        elseif(shapeType == 15)
            return read(io,PolygonZ{Float64,Float64})
        elseif(shapeType == 18)
            return read(io,MultiPointZ{Float64,Float64})
        elseif(shapeType == 21)
            return read(io,PointM{Float64,Float64})
        elseif(shapeType == 23)
            return read(io,PolylineM{Float64,Float64})
        elseif(shapeType == 25)
            return read(io,PolygonM{Float64,Float64})
        elseif(shapeType == 28)
            return read(io,MultiPointM{Float64,Float64})
        elseif(shapeType == 31)
            return read(io,MultiPatch{Float64,Float64})
        else
            error("Unknown shape type $shapeType")
        end
    end

    function read(io::IO,::Type{Handle})
        code = bswap(read(io,Int32))
        read(io,Int32,5)
        fileSize = bswap(read(io,Int32))
        version = read(io,Int32)
        shapeType = read(io,Int32)
        MBR = read(io,Rect{Float64})
        zmin = read(io,Float64)
        zmax = read(io,Float64)
        mmin = read(io,Float64)
        mmax = read(io,Float64)
        shapes = Array(ESRIShape,0)
        file = Handle(code,fileSize,version,shapeType,MBR,Interval(zmin,zmax),Interval(mmin,mmax),shapes)
        while(!eof(io))
            push!(shapes,read(io,ESRIShape))
        end
        file
    end

    #If Compose.jl is present, define useful interconversion functions
    isdefined(:Compose) && isa(Compose, Module) && include("compose.jl")
end # module


text:
"""
Generic Fermion module
"""
module AbstractFermion
    #export FermionFields,Wx!,Wdagx!,clear!,substitute_fermion!,Dx!

    #export FermionFields,
    #        Wx!,Wdagx!,clear!,substitute_fermion!,Dx!,fermion_shift!,
    #        fermion_shiftB!,add!,set_wing_fermi!,WdagWx!,apply_periodicity

    import ..Actions:FermiActionParam,FermiActionParam_Wilson,
                FermiActionParam_WilsonClover,FermiActionParam_Staggered
    import ..Gaugefields:GaugeFields,GaugeFields_1d,SU3GaugeFields,SU2GaugeFields,SU3GaugeFields_1d,SU2GaugeFields_1d,
                staggered_phase,SUn,SU2,SU3,SUNGaugeFields,SUNGaugeFields_1d,SU
    using Random
    abstract type FermionFields end

    function clear!(a::FermionFields)
        n1,n2,n3,n4,n5,n6 = size(a.f)
        for i6=1:n6
            for i5=1:n5
                for i4=1:n4
                    for i3=1:n3
                        for i2=1:n2
                            @simd for i1=1:n1
                                a.f[i1,i2,i3,i4,i5,i6]= 0
                            end
                        end
                    end
                end
            end
        end
    end

    function add!(c::FermionFields,alpha::Number,a::FermionFields,beta::Number,b::FermionFields) #c = c + alpha*a + beta*b
        n1,n2,n3,n4,n5,n6 = size(a.f)

        for i6=1:n6
            for i5=1:n5
                for i4=1:n4
                    for i3=1:n3
                        for i2=1:n2
                            @simd for i1=1:n1
                                #println(a.f[i1,i2,i3,i4,i5,i6],"\t",b.f[i1,i2,i3,i4,i5,i6] )
                                c.f[i1,i2,i3,i4,i5,i6] += alpha*a.f[i1,i2,i3,i4,i5,i6] +beta*b.f[i1,i2,i3,i4,i5,i6]
                            end
                        end
                    end
                end
            end
        end
        return
    end

    function add!(c::FermionFields,alpha::Number,a::FermionFields) #c = c + alpha*a
        n1,n2,n3,n4,n5,n6 = size(a.f)

        for i6=1:n6
            for i5=1:n5
                for i4=1:n4
                    for i3=1:n3
                        for i2=1:n2
                            @simd for i1=1:n1
                                #println(a.f[i1,i2,i3,i4,i5,i6],"\t",b.f[i1,i2,i3,i4,i5,i6] )
                                c.f[i1,i2,i3,i4,i5,i6] += alpha*a.f[i1,i2,i3,i4,i5,i6]
                            end
                        end
                    end
                end
            end
        end
        return
    end

    function add!(coeff::Number,c::FermionFields,alpha::Number,a::FermionFields) #c = coeff*c + alpha*a
        n1,n2,n3,n4,n5,n6 = size(a.f)

        for i6=1:n6
            for i5=1:n5
                for i4=1:n4
                    for i3=1:n3
                        for i2=1:n2
                            @simd for i1=1:n1
                                #println(a.f[i1,i2,i3,i4,i5,i6],"\t",b.f[i1,i2,i3,i4,i5,i6] )
                                c.f[i1,i2,i3,i4,i5,i6] = coeff*c.f[i1,i2,i3,i4,i5,i6] + alpha*a.f[i1,i2,i3,i4,i5,i6]
                            end
                        end
                    end
                end
            end
        end
        return
    end

    function Wx!(xout::FermionFields,U::Array{GaugeFields,1},
        x::FermionFields,temps::Array{FermionFields,1},fparam::T)  where T <: FermiActionParam
        error("This FermionFields type is not supported yet in function Wx!. The type is $(typeof(x))")
    end

    function Dx!(xout::T,U::Array{G,1},
        x::T,temps::Array{T,1}) where  {T <: FermionFields,G <: GaugeFields}
        error("This FermionFields type is not supported yet in function Dx!. The type is $(typeof(x))")

        return
    end

    function Ddagx!(xout::T,U::Array{G,1},
        x::T,temps::Array{T,1}) where  {T <: FermionFields,G <: GaugeFields}
        error("This FermionFields type is not supported yet in function Ddagx!. The type is $(typeof(x))")

        return
    end

    function Dxplus!(xout::T,ν::Int64,U::Array{G,1},
        x::T,temps::Array{T,1}) where  {T <: FermionFields,G <: GaugeFields}
        #temp = temps[4]
        error("This FermionFields type is not supported yet in function Dx!. The type is $(typeof(x))")

        return
    end

    function fermion_shiftB!(b::T,evensite,u::Array{GaugeFields{SU{NC}},1},μ,a::T)  where {NC,T <: FermionFields}
        error("This FermionFields type is not supported yet in function fermion_shiftB!. The type is $(typeof(x))")
    end

    function Wdagx!(xout::FermionFields,U::Array{GaugeFields,1},
        x::FermionFields,temps::Array{FermionFields,1},fparam::T)  where T <: FermiActionParam
        error("This FermionFields type is not supported yet in function Wdagx!. The type is $(typeof(x))")
    end

    function WdagWx!(xout::T,U::Array{G,1},x::T,temps::Array{T,1},fparam) where {T <: FermionFields,G <: GaugeFields}
        temp = temps[5]
        Wx!(temp,U,x,temps,fparam)

        Wdagx!(xout,U,temp,temps,fparam)

        return
    end

    function WdagWx!(xout::T,U::Array{G,1},x::T,temps::Array{T,1},fparam,indices) where {T <: FermionFields,G <: GaugeFields}
        WdagWx!(xout,U,x,temps,fparam)
        return
    end

    function substitute_fermion!(H,j,x::FermionFields)
        error("substitute_fermion! is not implemented for $(typeof(x))")
    end

    function set_wing_fermi!(a::FermionFields)
        NT = a.NT
        NZ = a.NZ
        NY = a.NY
        NX = a.NX
        NC = a.NC

        #!  X-direction
        for ialpha=1:4
            for it=1:NT
                for iz = 1:NZ
                    for iy=1:NY
                        for k=1:NC
                            a[k,0,iy,iz,it,ialpha] = a.BoundaryCondition[1]*a[k,NX,iy,iz,it,ialpha]
                        end
                    end
                end
            end
        end

        for ialpha=1:4
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for k=1:NC
                            a[k,NX+1,iy,iz,it,ialpha] = a.BoundaryCondition[1]*a[k,1,iy,iz,it,ialpha]
                        end
                    end
                end
            end
        end

        #Y-direction
        for ialpha = 1:4
            for it=1:NT
                for iz=1:NZ
                    for ix=1:NX
                        for k=1:NC
                            a[k,ix,0,iz,it,ialpha] = a.BoundaryCondition[2]*a[k,ix,NY,iz,it,ialpha]
                        end
                    end
                end
            end
        end

        for ialpha=1:4
            for it=1:NT
                for iz=1:NZ
                    for ix=1:NX
                        for k=1:NC
                            a[k,ix,NY+1,iz,it,ialpha] = a.BoundaryCondition[2]*a[k,ix,1,iz,it,ialpha]
                        end
                    end
                end
            end
        end

        for ialpha=1:4
            # Z-direction
            for it=1:NT
                for iy=1:NY
                    for ix=1:NX
                        for k=1:NC
                            a[k,ix,iy,0,it,ialpha] = a.BoundaryCondition[3]*a[k,ix,iy,NZ,it,ialpha]
                            a[k,ix,iy,NZ+1,it,ialpha] = a.BoundaryCondition[3]*a[k,ix,iy,1,it,ialpha]

                        end
                    end
                end
            end

            #T-direction
            for iz=1:NZ
                for iy=1:NY
                    for ix=1:NX
                        for k=1:NC
                            a[k,ix,iy,iz,0,ialpha] = a.BoundaryCondition[4]*a[k,ix,iy,iz,NT,ialpha]
                            a[k,ix,iy,iz,NT+1,ialpha] = a.BoundaryCondition[4]*a[k,ix,iy,iz,1,ialpha]
                        end
                    end
                end
            end

        end

    end

    function fermion_shift!(b::F,u::Array{GaugeFields{SU{NC}},1},μ,a::F) where {NC,F <: FermionFields}
        if μ == 0
            substitute!(b,a)
            return
        end

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        #NC = a.NC

        if μ > 0
            idel = zeros(Int64,4)
            idel[μ] = 1

            n6 = size(a.f)[6]
            for ialpha=1:n6
                for it=1:NT
                    it1 = it + idel[4]
                    for iz=1:NZ
                        iz1 = iz + idel[3]
                        for iy=1:NY
                            iy1 = iy + idel[2]
                            for ix=1:NX
                                ix1 = ix + idel[1]
                                for k1=1:NC
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    @simd for k2=1:NC
                                        b[k1,ix,iy,iz,it,ialpha] += u[μ][k1,k2,ix,iy,iz,it]*a[k2,ix1,iy1,iz1,it1,ialpha]
                                    end
                                end
                            end
                        end
                    end
                end
            end

        elseif μ < 0
            idel = zeros(Int64,4)
            idel[-μ] = 1
            n6 = size(b.f)[6]
            for ialpha =1:n6
                for it=1:NT
                    it1 = it - idel[4]
                    for iz=1:NZ
                        iz1 = iz -idel[3]
                        for iy=1:NY
                            iy1 = iy -idel[2]
                            for ix=1:NX
                                ix1 = ix -idel[1]

                                for k1=1:NC
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    @simd for k2=1:NC
                                        b[k1,ix,iy,iz,it,ialpha] += conj(u[-μ][k2,k1,ix1,iy1,iz1,it1])*a[k2,ix1,iy1,iz1,it1,ialpha]
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end

    end

    function fermion_shift!(b::F,evensite::Bool,u::Array{GaugeFields{SU{NC}},1},μ,a::F) where {NC,F <: FermionFields}
        if μ == 0
            substitute!(b,a)
            return
        end
        ibush = ifelse(evensite,0,1)

        NX = a.NX
        NY = a.NY
        NZ = a.NZ
        NT = a.NT
        #NC = a.NC

        if μ > 0
            idel = zeros(Int64,4)
            idel[μ] = 1

            n6 = size(a.f)[6]
            for ialpha=1:n6
                for it=1:NT
                    it1 = it + idel[4]
                    for iz=1:NZ
                        iz1 = iz + idel[3]
                        for iy=1:NY
                            iy1 = iy + idel[2]
                            xran =1+(1+ibush+iy+iz+it)%2:2:NX
                            for ix in xran
                                ix1 = ix + idel[1]
                                for k1=1:NC
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    @simd for k2=1:NC
                                        b[k1,ix,iy,iz,it,ialpha] += u[μ][k1,k2,ix,iy,iz,it]*a[k2,ix1,iy1,iz1,it1,ialpha]
                                    end
                                end
                            end
                        end
                    end
                end
            end

        elseif μ < 0
            idel = zeros(Int64,4)
            idel[-μ] = 1
            n6 = size(b.f)[6]
            for ialpha =1:n6
                for it=1:NT
                    it1 = it - idel[4]
                    for iz=1:NZ
                        iz1 = iz -idel[3]
                        for iy=1:NY
                            iy1 = iy -idel[2]
                            xran =1+(1+ibush+iy+iz+it)%2:2:NX
                            for ix in xran
                                ix1 = ix -idel[1]

                                for k1=1:NC
                                    b[k1,ix,iy,iz,it,ialpha] = 0
                                    @simd for k2=1:NC
                                        b[k1,ix,iy,iz,it,ialpha] += conj(u[-μ][k2,k1,ix1,iy1,iz1,it1])*a[k2,ix1,iy1,iz1,it1,ialpha]
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end

    end

    function apply_periodicity(ix,iy,iz,it,NX,NY,NZ,NT,BC)

        ixflag1 = (ix < 1)
        ixflag2 = (ix > NX)
        ix += ifelse(ixflag1,NX,0)
        ix += ifelse(ixflag2,-NX,0)

        sign =  ifelse(ixflag1 || ixflag2,BC[1],1)
        #sign =  ifelse(ixflag2,BC[1],1)

        iyflag1 = (iy < 1)
        iyflag2 = (iy > NY)

        iy += ifelse(iyflag1,NY,0)
        iy += ifelse(iyflag2,-NY,0)

        sign = ifelse(iyflag1 || iyflag2,BC[2],1)
        #sign = ifelse(iyflag2,BC[2],1)

        izflag1 = (iz < 1)
        izflag2 = (iz > NZ)

        iz += ifelse(izflag1,NZ,0)
        iz += ifelse(izflag2,-NZ,0)

        sign = ifelse(izflag1 || izflag2,BC[3],1)
        #sign = ifelse(izflag2,BC[3],1)

        itflag1 = (it < 1)
        itflag2 = (it > NT)

        it += ifelse(itflag1,NT,0)
        it += ifelse(itflag2,-NT,0)

        sign = ifelse(itflag1 || itflag2,BC[4],1)
        #sign = ifelse(itflag2,BC[4],1)

        return ix,iy,iz,it,sign
    end

    """
c-------------------------------------------------c
c     Random number function for Gaussian  Noise
    with σ^2 = 1/2
c-------------------------------------------------c
    """
    function gauss_distribution_fermi!(x::FermionFields)
        NC = x.NC
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        n6 = size(x.f)[6]
        σ = sqrt(1/2)

        for ialpha = 1:n6
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            for ic=1:NC

                                x[ic,ix,iy,iz,it,ialpha] = σ*randn()+im*σ*randn()
                            end
                        end
                    end
                end
            end
        end

        set_wing_fermi!(x)

        return
    end

    """
c-------------------------------------------------c
c     Random number function for Gaussian  Noise
    with σ^2 = 1/2
c-------------------------------------------------c
    """
    function gauss_distribution_fermi!(x::FermionFields,randomfunc,σ)
        NC = x.NC
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        n6 = size(x.f)[6]
        #σ = sqrt(1/2)

        for mu = 1:n6
            for ic=1:NC
                for it=1:NT
                    for iz=1:NZ
                        for iy=1:NY
                            for ix=1:NX
                                v1 = sqrt(-log(randomfunc()+1e-10))
                                v2 = 2pi*randomfunc()

                                xr = v1*cos(v2)
                                xi = v1 * sin(v2)

                                x[ic,ix,iy,iz,it,mu] = σ*xr + σ*im*xi
                            end
                        end
                    end
                end
            end
        end

        set_wing_fermi!(x)

        return
    end

    function gauss_distribution_fermi!(x::FermionFields,randomfunc)
        σ = 1
        gauss_distribution_fermi!(x,randomfunc,σ)
    end

    function gauss_distribution_fermi_Z2!(x::FermionFields)
        NC = x.NC
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        n6 = size(x.f)[6]
        #σ = sqrt(1/2)

        for mu = 1:n6
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            for ic=1:NC
                                x[ic,ix,iy,iz,it,mu] = rand([-1,1])
                            end
                        end
                    end
                end
            end
        end

        set_wing_fermi!(x)

        return
    end

    """
c-------------------------------------------------c
c     Random number function Z4  Noise
c     https://arxiv.org/pdf/1611.01193.pdf
c-------------------------------------------------c
    """
    function Z4_distribution_fermi!(x::FermionFields)
        NC = x.NC
        NX = x.NX
        NY = x.NY
        NZ = x.NZ
        NT = x.NT
        n6 = size(x.f)[6]
        θ = 0.0
        N::Int32 = 4
        Ninv = Float64(1/N)
        for ialpha = 1:n6
            for it=1:NT
                for iz=1:NZ
                    for iy=1:NY
                        for ix=1:NX
                            for ic=1:NC
                                θ = Float64(rand(0:N-1))*π*Ninv # r \in [0,π/4,2π/4,3π/4]
                                x[ic,ix,iy,iz,it,ialpha] = cos(θ)+im*sin(θ)
                            end
                        end
                    end
                end
            end
        end

        set_wing_fermi!(x)

        return
    end

end


text:
using JLD2
using Statistics

mutable struct TempPrint <: LoggerKeyData

    function TempPrint(logger_init_info)
        new()
    end
end

function lg_step!(self::TempPrint, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    if cur_step_total % 500 == 0
        println("At step: ", cur_step_total)
    end
end

function lg_episode_end!(self::TempPrint, cur_step_in_episode, cur_step_total)
    println("Finished episode: ", cur_step_in_episode)
end

function save_log(self::TempPrint, save_dict::Dict)
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity