text:
# coding: utf-8

# In[1]:

#merge sort is so called advanced sort
#its like a binary search in a way
#we keep dividing the list into two parts
#until it reaches base case when there are only two elements
#we sort the two elements and combine two sublists together
#we do a traversal on both sublists
#we run through elements one by one
#we compare two lists with each other and insert into the bigger list
#its like shuffling poker deck
#say u got 1 and 2 from two different sublists
#1 is smaller than 2
#we insert 1 back to the bigger list
#the second element is 1.5
#it is still smaller than 2
#we keep inserting elements from a sorted sublist
#until the element is larger than 2
#now we insert 2 into the bigger list
#we keep doing the same procedure until we come back to the whole list
#we do a final sort, its all done

#from my description, u can easily tell that we need recursion
#i assume that is why its slowing down the sorting

# In[2]:

function merge_sort(arr)

    #so this is to check the base case
    #we wanna make sure the base case is a sublist with two elements
    #otherwise we got nothing to sort
    if length(arr)>1

        #this is how we keep splitting the list into two halves
        ind=Int32(floor(length(arr)/2))

        #we define each half as an individual list
        #then we run the same trick on each half recursively
        left=arr[1:ind]
        right=arr[ind+1:end]
        left=merge_sort(left)
        right=merge_sort(right)

        #let me denote ind_left as the index in left list
        #ind_right as the index in right list
        #ind_combined as the index in right and left combined list
        ind_left=1
        ind_right=1
        ind_combined=1

        #now we have three scenarios
        #the first case is both indexes are under list length
        #so we can compare two elements
        #we put the larger one into a bigger list
        #assuming an item from right side goes to a bigger list
        #we increase the right side index by 1
        #and come back to the loop to see which element is bigger
        while (ind_left<=length(left)) && (ind_right<=length(right))
            if left[ind_left]<right[ind_right]
                arr[ind_combined]=left[ind_left]
                ind_left+=1
                ind_combined+=1
            else
                arr[ind_combined]=right[ind_right]
                ind_right+=1
                ind_combined+=1
            end
        end

        #the second case is when one index is outta list length
        #the other isnt
        #it can occur under two scenarios
        #one is that we use // method to divide a list into half
        #for a list with an odd number of elements
        #we definitely have one side bigger than the other
        #the other scenario is more common
        #assume it is an extreme case
        #we have two list, left and right
        #all the elements in left is smaller than the smallest element in right
        #so we have placed all elements in left into a bigger list
        #the index of right list hasnt increased even one bit
        #thats when we use the below section to add elements in right into a bigger list
        while (ind_left<=length(left)) && (ind_right>length(right))
            arr[ind_combined]=left[ind_left]
            ind_left+=1
            ind_combined+=1
        end

        #vice versa
        while (ind_left>length(left)) && (ind_right<=length(right))
            arr[ind_combined]=right[ind_right]
            ind_right+=1
            ind_combined+=1
        end
    end

    return arr

end

# In[3]:

for _ in 1:100

    test_arr=rand(1000)

    if !(merge_sort(test_arr)==sort(test_arr))
        printstyled("Erreur",color=:red)
    end

end

max_stars_repo_path: merge sort.jl

max_stars_repo_name: je-suis-tm/search-and-sort

max_issues_repo_path: merge sort.jl

max_issues_repo_name: HafizSabia/search-and-sort

max_forks_repo_path: merge sort.jl

max_forks_repo_name: HafizSabia/search-and-sort


text:
#run everything on one process
if nprocs() > 1
    rmprocs(workers())
end

#get the correct paths and filenames on all systems
testfolder = dirname(@__FILE__())
importsfile = joinpath(testfolder,"imports.jl")
utilfile = joinpath(testfolder,"testutil.jl")

#import all functionality on the local process
include(importsfile)
include(utilfile)

#make the run repeatable
srand(0)

println()
println("===========================================")
println("+++++++++++++++++++++++++++++++++++++++++++")
println("Running all tests with ",nprocs()," process")
println("+++++++++++++++++++++++++++++++++++++++++++")
println("===========================================")
println()

#all test categories
alltests = [
    "unit",
    "performance",
    "functionality"]

println("==================")
println("Running all tests:")
println("==================")

for t in alltests
    tfile = joinpath(testfolder,string(t,".jl"))
    println("  * $(tfile) *")
    include(tfile)
    println()
end

#run the runner tests on additional processes
addprocs(2)
println()
println("================================================")
println("++++++++++++++++++++++++++++++++++++++++++++++++")
println("Running the multiprocess code tests with ",nprocs()," processes")
println("++++++++++++++++++++++++++++++++++++++++++++++++")
println("================================================")
println()

for i=1:nworkers()
    remotecall_wait(include,workers()[i],importsfile)
    remotecall_wait(include,workers()[i],utilfile)
end

partests = [
    "unit/jobsegments",
    "unit/runners",
    "unit/mhrunner",
    "unit/smhrunner",
    "unit/gmhrunner",
    "unit/multiprocess",
    "functionality/sintest1"]

for t in partests
    tfile = joinpath(testfolder,string(t,".jl"))
    println("  * $(tfile) *")
    try
        include(tfile)
    catch e
        if nprocs() > 1
            rmprocs(workers())
        end
        throw(e)
    end
    println()
end

if nprocs() > 1
    rmprocs(workers())
end
sleep(1.0)

ext: jl


text:
module TestStaniczenkoRewire
    using BioEnergeticFoodWebs
    using Test

    A = [0 0 1 1 ; 0 0 0 1 ; 0 0 0 0 ; 0 0 0 0]
    expected_rg = [0 0 0 0 ; 0 0 0 0 ; 0 0 0 0 ; 0 0 1 0]
    parameters = model_parameters(A, rewire_method = :DO)
    @test BioEnergeticFoodWebs.rewiring_graph(parameters) == expected_rg
    expected_pl = [0 0 0 0 ; 0 0 1 0 ; 0 0 0 0 ; 0 0 0 0]
    @test BioEnergeticFoodWebs.potential_newlinks(A, expected_rg, parameters) == expected_pl
    expected_newA = [0 0 0 0 ; 0 0 1 1 ; 0 0 0 0 ; 0 0 0 0]
    parameters[:extinctions] = [1]
    @test BioEnergeticFoodWebs.Staniczenko_rewire(parameters) == expected_newA

end

module TestStaniczenkoRewire_TRissue
    using BioEnergeticFoodWebs
    using Test

    A = [0 0 0 0 0 0 ;
         0 0 0 0 0 0 ;
         0 1 0 0 0 0 ;
         0 0 1 0 0 0 ;
         1 0 0 0 0 0 ;
         1 0 0 1 0 0 ]
    expected_rg = [0 0 0 1 0 0 ;
                   0 0 0 0 0 0 ;
                   0 0 0 0 0 0 ;
                   0 0 0 0 0 0 ;
                   0 0 0 0 0 0 ;
                   0 0 0 0 0 0 ]
    parameters = model_parameters(A, rewire_method = :DO)
    @test BioEnergeticFoodWebs.rewiring_graph(parameters) == expected_rg
    expected_pl = zero(A)
    @test BioEnergeticFoodWebs.potential_newlinks(A, expected_rg, parameters) == expected_pl
    expected_newA = [0 0 0 0 0 0 ;
                     0 0 0 0 0 0 ;
                     0 1 0 0 0 0 ;
                     0 0 1 0 0 0 ;
                     1 0 0 0 0 0 ;
                     0 0 0 0 0 0 ]
    parameters[:extinctions] = [6]
    @test BioEnergeticFoodWebs.Staniczenko_rewire(parameters) == expected_newA

end

module TestStaniczenkoRewire_exemplepaper
    using BioEnergeticFoodWebs
    using Test

    A = [0 0 0 0 0 0 0 ;
         0 0 0 0 0 0 0 ;
         0 0 0 0 0 0 0 ;
         1 0 0 0 0 0 0 ;
         1 1 0 0 0 0 0 ;
         0 1 1 0 0 0 0 ;
         0 0 0 1 1 1 0 ]
    expected_R = [0 1 0 0 0 0 0 ;
                  1 0 1 0 0 0 0 ;
                  0 0 0 0 0 0 0 ;
                  0 0 0 0 0 0 0 ;
                  0 0 0 0 0 0 0 ;
                  0 0 0 0 0 0 0 ;
                  0 0 0 0 0 0 0 ]
    new_A = [0 0 0 0 0 0 0 ;
             0 0 0 0 0 0 0 ;
             0 0 0 0 0 0 0 ;
             0 0 0 0 0 0 0 ;
             1 1 0 0 0 0 0 ;
             1 1 1 0 0 0 0 ;
             0 0 0 0 1 1 0 ]
    p = model_parameters(A, rewire_method = :DO)
    p[:extinctions] = [4]
    @test BioEnergeticFoodWebs.rewiring_graph(p) == expected_R
    @test BioEnergeticFoodWebs.Staniczenko_rewire(p) == new_A
end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
#####
##### Schedules
#####

"""
    AbstractSchedule{R}

Abstract type for a parameter schedule, which represents a function from
nonnegative integers to numbers of type `R`. Subtypes must implement the
`getindex(s::AbstractSchedule, i::Int)` operator.
"""
abstract type AbstractSchedule{R} end

#####
##### Constant Schedule
#####

struct ConstSchedule{R} <: AbstractSchedule{R}
  value :: R
end

Base.getindex(s::ConstSchedule, i::Int) = s.value

Base.convert(::Type{ConstSchedule{R}}, x::R) where {R <: Number} =
  ConstSchedule(x)

#####
##### Piecewise linar schedule
#####

"""
    PLSchedule{R} <: AbstractSchedule{R}

Type for piecewise linear schedules.

# Constructors

    PLSchedule(cst)

Return a schedule with a constant value `cst`.

    PLSchedule(xs, ys)

Return a piecewise linear schedule such that:
  - For all `i`, `(xs[i], ys[i])` belongs to the schedule's graph.
  - Before `xs[1]`, the schedule has value `ys[1]`.
  - After `xs[end]`, the schedule has value `ys[end]`.
"""
struct PLSchedule{R} <: AbstractSchedule{R}
  # We keep the internal representation simple for JSON serialization
  xs :: Vector{Int}
  ys :: Vector{R}
  function PLSchedule{R}(xs, ys) where R
    @assert !isempty(xs)
    @assert length(xs) == length(ys)
    new{R}(xs, ys)
  end
end

PLSchedule(xs, ys) = PLSchedule{eltype(ys)}(xs, ys)

PLSchedule(cst) = PLSchedule([0], [cst])

function Base.getindex(s::PLSchedule{R}, i::Int) where R
  ptidx = findlast(x -> x <= i, s.xs)
  if isnothing(ptidx)
    # We are before the first point
    return s.ys[1]
  elseif ptidx == length(s.xs)
    # We are past the last point
    return s.ys[end]
  else
    # We are between two points
    x0, y0 = s.xs[ptidx], s.ys[ptidx]
    x1, y1 = s.xs[ptidx+1], s.ys[ptidx+1]
    y = y0 + (y1 - y0) / (x1 - x0) * (i - x0)
    R <: Integer && (y = ceil(Int, y))
    return y
  end
end

function test()
  s = PLSchedule([0, 10, 20], [0, 10, 30])
  xs = [-1, 0, 2, 10, 11, 20, 25]
  ys = [0,  0, 2, 10, 12, 30, 30]
  @assert [s[x] for x in xs] == ys
end

#test()

#####
##### Step function
#####

"""
    StepSchedule{R} <: AbstractSchedule{R}

Type for step function schedules.

# Constructor

    StepSchedule(;start, change_at, values)

Return a schedule that has initial value `start`. For all `i`, the schedule
takes value `values[i]` at step `change_at[i]`.
"""
struct StepSchedule{R} <: AbstractSchedule{R}
  start :: R
  xs :: Vector{Int}
  ys :: Vector{R}
  function StepSchedule{R}(start, xs, ys) where R
    @assert length(xs) == length(ys)
    return new{R}(start, xs, ys)
  end
end

StepSchedule(; start, change_at, values) =
  StepSchedule{typeof(start)}(start, change_at, values)

function Base.getindex(s::StepSchedule, i::Int)
   idx = findlast(x -> x <= i, s.xs)
   isnothing(idx) && (return s.start)
   return s.ys[idx]
end

"""
    CyclicSchedule(base, mid, term; n, xmid=0.45, xback=0.90)

Return the [`PLSchedule`](@ref) that is typically used for cyclic
learning rate scheduling.
"""
function CyclicSchedule(base, mid, term; n, xmid=0.45, xback=0.90)
  nmid  = floor(Int, xmid * n)
  nback = floor(Int, xback * n)
  return PLSchedule([1, nmid, nback, n], [base, mid, base, term])
end


text:
push!(LOAD_PATH, "./")
push!(LOAD_PATH, "../")
push!(LOAD_PATH, "../logs/")
push!(LOAD_PATH, "../../VICAS/src/")

using PyPlot
using DeepQLearning
using Simulation_Env_module
using CAS_QMDP
using Interact

rc("font", family="Times New Roman")
rc("text", usetex=true)

function get_section_id(x::Float64, y::Float64, num_sections::Int64)
	ac_pos_angle = norm_angle(atan2(y, x))
	return floor(ac_pos_angle / (2 * pi / num_sections)) + 1
end

function viz_policy_multi(policy::DQNPolicy;
	CASType::AbstractString="",
	state_formation::Symbol=:closest,
	correctionOffOption::Symbol=:multi,
	num_intruders::Int64=1,
	resolution::Int64=100,
	correction::Bool=true,
	colorbarOption::Bool=false,
	xlabelOption::Bool=true,
	ylabelOption::Bool=false,
	legendOption::Bool=false,
	labelleft::Bool=true,
	labelbottom::Bool=true,
	fontsize::Int64=12,
	titleOption::Bool=true,
    format::AbstractString="pdf",
    figsize::Float64=3.5)

	rc("font", size=fontsize)

	XMIN = -1200
	XMAX = 1200
	YMIN = XMIN
	YMAX = XMAX
	v = 50.0
	viz = figure()

	save_button = button("Save Policy Slice", value=0)
	display(save_button)
	fig_name_input = textbox("Figure Name", value=CASType*"PolicySlice")
	display(fig_name_input)
	print_state_button = button("Print full_state", value=0)
	display(print_state_button)
	int_pos_input = textbox("Intruder Position", value="0 0")
	display(int_pos_input)

	@manipulate for int_1_x = XMIN : 150 : XMAX,
					int_1_y = YMIN : 150 : YMAX,
					int_1_p = 0. : 30. : 360,
					int_2_x = XMIN : 150 : XMAX,
					int_2_y = YMIN : 150 : YMAX,
					int_2_p = 0. : 30. : 360,
					p = 0. : 30. : 360,
					deviation = -180. : 30. : 180.,
					dist_to_dest = 100. : 500. : DEST_RANGE
					# ...

		num_fix_intruders = 0
		lowfi_state_matrix = ones(STATE_DIM, 5) * TERM_VAR # 5 is a defautl max num of intruders

		withfig(viz) do
			int_1_rho = sqrt(int_1_x^2 + int_1_y^2)
			if int_1_rho <= SENSING_RANGE
				lowfi_state_matrix[:, 1] .= [int_1_rho, norm_angle(atan2(int_1_y, int_1_x)), deg2rad(int_1_p), v, v]
				num_fix_intruders += 1
			end
			int_2_rho = sqrt(int_2_x^2 + int_2_y^2)
			if int_2_rho <= SENSING_RANGE
				lowfi_state_matrix[:, 2] .= [int_2_rho, norm_angle(atan2(int_2_y, int_2_x)), deg2rad(int_2_p), v, v]
				num_fix_intruders += 1
			end
			# ...

			function get_heat(x::Float64, y::Float64) # of the moving intruder
    			rho = sqrt(x^2 + y^2)
    			theta = norm_angle(atan2(y, x))
    			if rho <= SENSING_RANGE
    				num_intruders = num_fix_intruders + 1
    				lowfi_state_matrix[:, num_intruders] .= [rho, theta, deg2rad(p), v, v]
    			else # if rho >= SENSING_RANGE
    				num_intruders = num_fix_intruders
    			end

    			if num_intruders > 0
					closest_intruder_ind = indmin(lowfi_state_matrix[IND_rho, :])
    				low_fi_state = vec(lowfi_state_matrix[:, closest_intruder_ind])

    				if correction == true
    					full_state = ones(policy.env.num_sections * STATE_DIM + AUG_STATE_DIM) * TERM_VAR

    					# full_state formualtion method: 4 closest intruders
    					# Sort lowfi_state_matrix by rho, resulting in acsending order
    					if state_formation == :closest
    						temp_lowfi = copy(lowfi_state_matrix)
    						# lowfi_state_matrix = transpose(sortrows(transpose(lowfi_state_matrix), by=x->x[IND_rho]))
    						for i in 1 : policy.env.num_sections
    							# if !is_terminal(vec(lowfi_state_matrix[:, i]))
    							closestInd = indmin(temp_lowfi[IND_rho, :])
    							if temp_lowfi[IND_rho, closestInd] <= SENSING_RANGE
    								full_state[STATE_DIM * (i - 1) + 1 : STATE_DIM * i] = vec(temp_lowfi[:, closestInd])
    								temp_lowfi[IND_rho, closestInd] = Inf
    							end
    						end

    					else # full_state formulation method: closest intruders in each section
    						section_ac_arr = Dict{Int64, Array{Float64,2}}()
    						for isection in 1 : policy.env.num_sections
    							section_ac_arr[isection] = Array{Float64,2}(STATE_DIM, 0)
    						end

    						for i in 1 : num_intruders
    							# if !is_terminal(vec(lowfi_state_matrix[:, i]))
    							if vec(lowfi_state_matrix[:, i])[IND_rho] <= SENSING_RANGE
    								ix = lowfi_state_matrix[IND_rho, i] * cos(lowfi_state_matrix[IND_theta, i])
    								iy = lowfi_state_matrix[IND_rho, i] * sin(lowfi_state_matrix[IND_theta, i])
    								section_id_i = get_section_id(ix, iy, policy.env.num_sections)
    								section_ac_arr[section_id_i] = hcat(section_ac_arr[section_id_i], lowfi_state_matrix[:, i])
    							end
    						end

    						for isection in 1 : policy.env.num_sections
    							if length(section_ac_arr[isection]) > 0
    								closest_intruder_ind = indmin(section_ac_arr[isection][IND_rho, :])
    								full_state[STATE_DIM * (isection - 1) + 1 : STATE_DIM * isection] .=
										vec(section_ac_arr[isection][:, closest_intruder_ind])
    							end
    						end
						end # if state_formation

						# Augment states
    					prev_dist_to_dest = dist_to_dest + v
    					full_state[end - AUG_STATE_DIM + 1 : end] .= [deg2rad(deviation), dist_to_dest, prev_dist_to_dest]

    					normalize_obs!(policy.env, full_state)
    					action = get_action(policy, full_state, low_fi_state) # returns symbol

     				else # if correction == false
     					full_state = ones(policy.env.num_sections * STATE_DIM + AUG_STATE_DIM) * TERM_VAR # dummy variable
     					if correctionOffOption == :multi
     						q_matrix = Matrix{Float64}(length(policy.env.action_space), 0)
     						for i in 1 : num_intruders
     							if vec(lowfi_state_matrix[:, i])[IND_rho] <= SENSING_RANGE
     								state_i = vec(lowfi_state_matrix[:, i])
     								qvals = get_qvals(policy.env.low_fi_policy, state_i)
     								q_matrix = hcat(q_matrix, qvals)
     							end
     						end
     						if size(q_matrix, 2) == 0
     							q_matrix = hcat(q_matrix, terminal_state)
     						end
     						min_qvals = minimum(q_matrix, 2)
							the_ind = indmax(min_qvals)

     						action = Simulation_Env_module.actions(policy.env)[the_ind]

     					elseif correctionOffOption == :closest
							qvals = get_qvals(policy.env.low_fi_policy, low_fi_state)
							action = Simulation_Env_module.actions(policy.env)[indmax(qvals)]
     					end
     				end # if

     				if action == :COC
    					return 2.0, full_state, lowfi_state_matrix
    				end
     				return rad2deg(policy.env.advisory_to_action_dict[action]), full_state, lowfi_state_matrix # get numerical value of action

     			else # if num_intruders == 0, i.e. no intruder
     				return 2.0, full_state, lowfi_state_matrix
     			end
			end # get_heat

			# plot
			action_map = zeros(resolution, resolution)
			x_arr = linspace(XMIN, XMAX, resolution)
			y_arr = linspace(YMIN, YMAX, resolution)

			for j in 1:resolution
				for i in 1:resolution
					action_map[j, i], _, _ = get_heat(x_arr[i], y_arr[j])
				end
			end

			action_map[Int(resolution / 2), Int(resolution / 2)] = 10.0
			action_map[Int(resolution / 2), Int(resolution / 2)+1] = -10.0

			action_map = flipdim(action_map, 1)

			g = imshow(action_map, cmap="jet", extent=(XMIN, XMAX, YMIN, YMAX))

			if colorbarOption == true
				cbar = colorbar(g, fraction=0.088, pad=0.04, ticks=[-10, -5, 0, 2, 5, 10], aspect=10, label="Advisory of Ownship")
				cbar[:set_ticklabels](["SR", "WR", "KEEP", "COC", "WL", "SL"])
			end

			# plot Ownship triangle
			scatter(0, 0 , marker=">", color="green", s=180, label="Ownship")

			# plot intruder arrow
			airspace_dim = XMAX - XMIN
            arrow_len = airspace_dim / 12
            arrow(0.8 * XMAX, 0.8 * YMAX,
                arrow_len * cos(deg2rad(p)), # dx
                arrow_len * sin(deg2rad(p)), # dy
                head_width = airspace_dim / 30,
                width = airspace_dim / 90,
                head_length = airspace_dim / 25,
                head_starts_at_zero = "true",
                facecolor = "red",
                length_includes_head = "true")

            scatter(0.8 * XMAX, 0.8 * XMAX , color="red", s=40, label="Free Intruder")

			# plot existing intruder 1 arrow
			airspace_dim = XMAX - XMIN
			arrow_len = airspace_dim / 12

			# plot fixed intruder 1
			arrow(int_1_x, int_1_y,
				arrow_len * cos(deg2rad(int_1_p)), # dx
				arrow_len * sin(deg2rad(int_1_p)), # dy
				head_width = airspace_dim / 20,
				width = airspace_dim / 70,
				head_length = airspace_dim / 18,
				head_starts_at_zero = "true",
				facecolor = "blue",
				length_includes_head = "true")
			scatter(int_1_x, int_1_y , color="blue", s=60, label="Fixed Intruder 1")

			# plot fixed intruder 2
			arrow(int_2_x, int_2_y,
				arrow_len * cos(deg2rad(int_2_p)), # dx
				arrow_len * sin(deg2rad(int_2_p)), # dy
				head_width = airspace_dim / 20,
				width = airspace_dim / 70,
				head_length = airspace_dim / 18,
				head_starts_at_zero = "true",
				facecolor = "magenta",
				length_includes_head = "true")
			scatter(int_2_x, int_2_y , color="magenta", s=60, label="Fixed Intruder 2")
			# ...

			# plot NMAC circle
			th = linspace(0, 2*pi, 50)
			PyPlot.plot(NMAC_RANGE * cos.(th), NMAC_RANGE * sin.(th), linestyle="--", color="red", label="NMAC Range = 150 m")

			# plot sensing circle
			PyPlot.plot(SENSING_RANGE * cos.(th), SENSING_RANGE * sin.(th), linestyle="--", color="blue", label="Sensing Range = 1000 m")

			if ylabelOption == true
				ylabel("\$y\$ (m)")
			end
			if xlabelOption == true
				xlabel("\$x\$ (m)")
			end

			if titleOption == true
				title(CASType)
			end

    		ax = gca()
    		ax[:set_xlim]([-1200,1200])
    		ax[:set_ylim]([-1200,1200])
    		ax[:tick_params](direction="in", color="gray", labelleft=labelleft, labelbottom=labelbottom)
    		fig = gcf()
			fig[:set_size_inches](figsize, figsize)

			if legendOption == true
				legend(loc="center left", bbox_to_anchor=(1.05, 0.2), fancybox=false, edgecolor="black", shadow=false)
			end

			o_saveButton = observe(save_button)
    		o_figNameTextbox = observe(fig_name_input)
			function save_slice(signal::Int64, filename=o_figNameTextbox[])
				if signal > 0
					savefig(filename * "." * format, bbox_inches="tight", pad_inches=0.01)
				end
			end
			on(save_slice, o_saveButton)

    		o_button = observe(print_state_button)
    		o_textbox = observe(int_pos_input)
			function print_state(signal::Int64, pos_str=o_textbox[])
				pos = [parse(Float64, ps) for ps in split(pos_str)]
				if signal > 0
					_, fs, lowfi_state_matrix = get_heat(pos[1], pos[2])
					println("fs = ", fs)
					println("lowfi_state_matrix = ", lowfi_state_matrix)
					println()
				end
			end
			on(print_state, o_button)

    		return viz

		end # withfig
	end
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
using Module1
using Module2: func1, func2
import Module3
import Module4.func1, Module5, Module6.func2
import Module7: func1, func2

module MyModule
using Module8,
      Module9,
      Module10

using Module11: func1,
                func2,
                func3
end


text:
using ClinicalTrialSampleSize
using Base.Test

# OneSampleMean
# Constructor: OneSampleMean()
@test_throws ErrorException OneSampleMean(-Inf, Inf)

params = (
    (OneSampleMean(3, 2), 7.849, 0.8, 1, 0.05, "two"),
    (OneSampleMean(3, 2.5, true), 72.606, 0.8, 1.5, 0.05, "two"),
    (OneSampleMean(10, 7, true), 25.448, 0.9, 5, 0.05, "one"),
    (OneSampleMean(7, 4, false), 6.569, 0.9, 3, 0.1, "one"),
)

for (test, n, p, std, alpha, side) in params

        p_est = power(test; n = n, std = std, alpha = alpha, side = side)
        n_est = sample_size(test; power = p, std = std, alpha = alpha, side = side)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

# One-sample non-inferiority/superiority for mean
# Constructor: OneSampleMeanInferior/Superior
@test_throws ErrorException OneSampleMeanInferior(-Inf, Inf, -1)
@test_throws ErrorException OneSampleMeanInferior(30, 35, 10)
@test_throws ErrorException OneSampleMeanSuperior(-Inf, Inf, 1)
@test_throws ErrorException OneSampleMeanSuperior(35, 30, -10)

params = (
    (OneSampleMeanInferior(2, 2, -1), 6.183, 0.8, 1, 0.05),
    (OneSampleMeanInferior(1.9, 2.5, -1, true), 11.255, 0.8, 0.5, 0.05),
    (OneSampleMeanSuperior(20, 15, 3), 54.369, 0.85, 5.5, 0.05),
    (OneSampleMeanSuperior(121, 119, 2.3, true), 857.934, 0.95, 3, 0.1),
)

for (test, n, p, std, alpha) in params

        p_est = power(test; n = n, std = std, alpha = alpha)
        n_est = sample_size(test; power = p, std = std, alpha = alpha)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

# One-sample equivalence for mean
# Constructor: OneSampleMeanEqual()
@test_throws ErrorException OneSampleMeanEqual(-Inf, Inf, 1)
@test_throws ErrorException OneSampleMeanEqual(30, 30, -1)

params = (
    (OneSampleMeanEqual(2, 2, 0.3), 95.154, 0.8, 1, 0.05),
    (OneSampleMeanEqual(1.94, 1.95, 0, true), 5947.569, 0.85, 0.25, 0.05),
)

for (test, n, p, std, alpha) in params

        p_est = power(test; n = n, std = std, alpha = alpha)
        n_est = sample_size(test; power = p, std = std, alpha = alpha)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

# One-sample proportion
# Constructor: OneSampleProp()
@test_throws ErrorException OneSampleProp(-1.2, 1.2)

params = (
    (OneSampleProp(0.5, 0.2), 21.802, 0.8, 0.05, "two"),
    (OneSampleProp(0.9, 0.4), 3.232, 0.85, 0.05, "two"),
    (OneSampleProp(0.44, 0.4), 1666.615, 0.95, 0.1, "two"),
    (OneSampleProp(0.7, 0.5), 49.682, 0.9, 0.05, "one"),
    (OneSampleProp(0.32, 0.27), 369.834, 0.8, 0.1, "one"),
)

for (test, n, p, alpha, side) in params

        p_est = power(test; n = n, alpha = alpha, side = side)
        n_est = sample_size(test; power = p, alpha = alpha, side = side)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

# One-sample non-inferiority/superiority proportion
# Constructor: OneSamplePropInferior/Superior
@test_throws ErrorException OneSamplePropInferior(-1.2, 1.2, 0)
@test_throws ErrorException OneSamplePropInferior(0.3, 0.2, 1)
@test_throws ErrorException OneSamplePropSuperior(1.3, 0.3, 0)
@test_throws ErrorException OneSamplePropSuperior(0.4, 0.2, 1.2)

params = (
    (OneSamplePropInferior(0.25, 0.3, -0.06), 11592.29, 0.8, 0.05),
    (OneSamplePropInferior(0.71, 0.8, -0.05), 691.445, 0.85, 0.1),
    (OneSamplePropInferior(0.4, 0.4, -0.02), 5138.308, 0.95, 0.1),
    (OneSamplePropSuperior(0.7, 0.5, 0.1), 129.834, 0.8, 0.05),
    (OneSamplePropSuperior(0.32, 0.27, 0.03), 2922.942, 0.85, 0.1),
    (OneSamplePropSuperior(0.5, 0.44, 0.2), 57.2, 0.9, 0.2),
)

for (test, n, p, alpha) in params

        p_est = power(test; n = n, alpha = alpha)
        n_est = sample_size(test; power = p, alpha = alpha)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

# One-sample equivalence proportion
# Constructor: OneSamplePropEqual
@test_throws ErrorException OneSamplePropEqual(1.9, 1.2, 0.1)
@test_throws ErrorException OneSamplePropEqual(0.2, 0.2, -1)

params = (
    (OneSamplePropEqual(0.35, 0.3, 0.06), 19482.75, 0.8, 0.05),
    (OneSamplePropEqual(0.85, 0.8, 0.1), 377.619, 0.85, 0.1),
    (OneSamplePropEqual(0.4, 0.4, 0.02), 6304.454, 0.95, 0.1),
)

for (test, n, p, alpha) in params

        p_est = power(test; n = n, alpha = alpha)
        n_est = sample_size(test; power = p, alpha = alpha)

        @test abs(p_est - p) < 1e-3
        @test abs(n_est - n) < 1

end # end for

return true


text:
export NumericParser, DateParser, DateTimeParser
import TextParse: Numeric
import Dates
import Dates: DateFormat

NumericParser(x...) = parser(TextParse.Numeric(x...))
DateParser(format::AbstractString...) = DateParser(Dates.DateFormat.(format)...)
DateParser(format::DateFormat...) = sEither(parser.(TextParse.DateTimeToken.(Dates.Date,format))...)
DateTimeParser(format::AbstractString...) = DateTimeParser(Dates.DateFormat.(format)...)
DateTimeParser(format::DateFormat...) = sEither(parser.(TextParse.DateTimeToken.(Dates.DateTime,format))...)

import TextParse: tryparsenext

result_type(x::AbstractToken) = result_type(typeof(x))
result_type(::Type{<:AbstractToken{T}}) where T = T

struct AbstractTokenParser{P<:AbstractToken,T} <: LeafParser{NCodeunitsState{T},T}
    parser::P
    function AbstractTokenParser(p::AbstractToken)
        new{typeof(p), result_type(p)}(p)
    end
end

parser(x::AbstractToken) = AbstractTokenParser(x)

regex_string(::TextParse.Numeric{<:Integer}) = "-?[[:digit:]]+"

print_constructor(io::IO, x::AbstractTokenParser) =
    print(io, x.parser)

_iterate(parser::AbstractTokenParser, sequence, till, before_i, next_i, state) =
    _iterate_token(parser.parser, sequence, till, before_i, next_i, state)

_iterate_token(parser::AbstractToken, sequence, till, before_i, next_i, state) =
    nothing
function _iterate_token(parser::AbstractToken, sequence, till, before_i, next_i, state::Nothing, opts=TextParse.default_opts)
        r,next_i_ = tryparsenext(parser, sequence, next_i, till,opts)
        if isnull(r)
            nothing
        else
            NCodeunitsState(next_i,next_i_,get(r))
        end
end

"""
    TextParse.tryparsenext(x::CombinedParser,str,i,till,opts=TextParse.default_opts)

TextParse.jl integrates with CombinedParsers.jl both ways.

```jldoctest
julia> using TextParse

julia> p = ("Number:" * Repeat(' ') * TextParse.Numeric(Int))[3]
🗄 Sequence[3]
├─ Number\\:
├─ \\ *  |> Repeat
└─ <Int64>
::Int64

julia> parse(p, "Number:    42")
42

julia> TextParse.tryparsenext(p, "Number:    42")
(Nullable{Int64}(42), 14)
```

"""
function TextParse.tryparsenext(x::CombinedParser,str,i,till,opts=TextParse.default_opts)
    s = _iterate(x,str,till,i,nothing)
    if s === nothing
        Nullable{result_type(x)}(),i
    else
        Nullable(get(x,str,till,tuple_pos(s),i,tuple_state(s))),tuple_pos(s)
    end
end

"""
    regex_prefix(x)

Prefix printed in parser tree node.
"""
regex_prefix(x::AbstractToken) = ""
"""
    regex_suffix(x)

Suffix printed in parser tree node.
"""
regex_suffix(x::AbstractToken) = ""
"""
    regex_inner(x::AbstractToken)

Regex representation of `x`.
See [`regex_string`](@ref)
"""
regex_inner(x::AbstractToken) = "$(typeof(x))"
regex_inner(::TextParse.Numeric{T}) where T = "$(T)"


text:
abstract type SVAE end

"""
	kldiv(model::SVAE, κ)

	KL divergence between Von Mises-Fisher and Hyperspherical Uniform distributions
"""
kldiv(model::SVAE, κ) = .- vmfentropy(model.zdim, κ) .+ model.hue

function log_pxexpectedz(m::SVAE, x)
	xgivenz = m.g(zparams(m, x)[1])
	Flux.Tracker.data(log_normal(x, xgivenz))
end

jacobian_encoder(m::SVAE, x) = Flux.Tracker.jacobian(a -> Chain(m.q, m.μzfromhidden)(a), x)

function log_pz_jacobian_encoder(m::SVAE, x)
	if size(x, 2) > 1
		xs = [x[:, i] for i in 1:size(x, 2)]
		return map(x -> log_pz_jacobian_encoder_singleinstance(m, x), xs)
	else
		return log_pz_jacobian_encoder_singleinstance(m, x)
	end
end

function log_pz_jacobian_encoder_singleinstance(m::SVAE, x)
	@assert size(x, 2) == 1
	s = svd(jacobian_encoder(m, x).data)
	d = reduce(+, log.(abs.(s.S))) * 2
	d + log_pz(m, x)
end

jacobian_decoder(m::SVAE, z) = Flux.Tracker.jacobian(a -> m.g(a), z)

function log_pz_jacobian_decoder(m::SVAE, z)
	if size(z, 2) > 1
		zs = [z[:, i] for i in 1:size(z, 2)]
		return map(z -> log_pz_jacobian_decoder_singleinstance(m, z), zs)
	else
		return log_pz_jacobian_decoder_singleinstance(m, z)
	end
end

function log_pz_jacobian_decoder_singleinstance(m::SVAE, z)
	@assert size(z, 2) == 1
	s = svd(jacobian_decoder(m, z).data)
	d = reduce(+, log.(abs.(s.S))) * 2
	-d + log_pz_from_z(m, z)
end

"""
	infer(m::SVAE, x)

	infer latent variables and sample output x
"""
function infer(m::SVAE, x)
	(μz, κz) = zparams(m, x)
	z = samplez(m, μz, κz)
	xgivenz = m.g(z)
	return xgivenz
end

"""
	zparams(model::SVAE, x)

	Computes μ and κ from the hidden layer
"""

function zparams(model::SVAE, x)
	hidden = model.q(x)
	# return model.μzfromhidden(hidden), max.(model.κzfromhidden(hidden), 100)
	return model.μzfromhidden(hidden), model.κzfromhidden(hidden)
end

"""
	zfromx(m::SVAE, x)

	convenience function that returns latent layer based on the input `x`
"""
zfromx(m::SVAE, x) = samplez(m, zparams(m, x)...)

"""
	samplez(m::SVAE, μz, κz)

	samples z layer based on its parameters
"""
function samplez(m::SVAE, μz, κz)
	ω = sampleω(m, κz)
	normal = Normal()
	v = Adapt.adapt(eltype(Flux.Tracker.data(κz)), rand(normal, size(μz, 1) - 1, size(μz, 2)))
	v = normalizecolumns(v)
	z = householderrotation(vcat(ω, sqrt.(1 .- ω .^ 2 .+ (10 * eps(Float32))) .* v), μz)
	return z
end

"""
	householderrotation(zprime, μ)

	uses Householder reflection to rotate the `zprime` vector according to mapping of e1 to μ
"""
function householderrotation(zprime, μ)
	e1 = similar(μ) .= 0
	e1[1, :] .= 1
	u = e1 .- μ
	normalizedu = normalizecolumns(u)
	return zprime .- 2 .* sum(zprime .* normalizedu, dims = 1) .* normalizedu
end

function sampleω(model::SVAE, κ)
	m = model.zdim
	c = @. sqrt(4κ ^ 2 + (m - 1) ^ 2)
	b = @. (-2κ + c) / (m - 1)
	a = @. (m - 1 + 2κ + c) / 4
	d = @. (4 * a * b) / (1 + b) - (m - 1) * log(m - 1)
	ω = rejectionsampling(m, a, b, d, κ)
	return ω
end

function rejectionsampling(m, a, b, d, κ)
	beta = Beta((m - 1) / 2, (m - 1) / 2)
	T = eltype(Flux.Tracker.data(a))
	ϵ, u = Adapt.adapt(T, rand(beta, size(a)...)), rand(T, size(a)...)

	accepted = isaccepted(ϵ, u, m, Flux.data(a), Flux.Tracker.data(b), Flux.data(d))
	it = 0
	while (!all(accepted)) & (it < 10000)
		mask = .! accepted
		ϵ[mask] = Adapt.adapt(T, rand(beta, sum(mask)))
		u[mask] = rand(T, sum(mask))
		accepted[mask] = isaccepted(mask, ϵ, u, m, Flux.data(a), Flux.data(b), Flux.data(d))
		it += 1
	end
	if it >= 10000
		println("Warning - sampler was stopped by 10000 iterations - it did not accept the sample!")
		println(m)
		println(a)
		println(b)
		println(d)
		println(κ)
		println(accepted)
	end
	return @. (1 - (1 + b) * ϵ) / (1 - (1 - b) * ϵ)
end

# function rejectionsampling(m, a, b, d)
# 	beta = Beta((m - 1) / 2, (m - 1) / 2)
# 	T = eltype(Flux.Tracker.data(a))
# 	ϵ, u = Adapt.adapt(T, rand(beta, size(a)...)), rand(T, size(a)...)
#
# 	accepted = isaccepted(ϵ, u, m, Flux.data(a), Flux.Tracker.data(b), Flux.data(d))
# 	while !all(accepted)
# 		mask = .! accepted
# 		ϵ[mask] = Adapt.adapt(T, rand(beta, sum(mask)))
# 		u[mask] = rand(T, sum(mask))
# 		ia = isaccepted(mask, ϵ, u, m, Flux.data(a), Flux.data(b), Flux.data(d))
# 		accepted[mask] = ia
# 	end
# 	return @. (1 - (1 + b) * ϵ) / (1 - (1 - b) * ϵ)
# end

isaccepted(mask, ϵ, u, m:: Int, a, b, d) = isaccepted(ϵ[mask], u[mask], m, a[mask], b[mask], d[mask]);
function isaccepted(ϵ, u, m:: Int, a, b, d)
	ω = @. (1 - (1 + b) * ϵ) / (1 - (1 - b) * ϵ)
	t = @. 2 * a * b / (1 - (1 - b) * ϵ)
	@. (m - 1) * log(t) - t + d >= log(u)
end

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
using LinearAlgebra
using Statistics
using StatsBase

function linear_regression(X, y)
    Q, R = modified_gram_schmidt(X)
    f = transpose(Q) * y
    betas = backsolve(R, f)
    return betas
end

# Helper Functions

function back_substitute(U, y)
    n, _ = size(U)
    x = zeros(n)
    for i in n:1
        x[i] = y[i]
        for j in (i+1): n
            x[i] = x[i] - U[i, j] * x[j]
        end
        x[i] = x[i] / U[i, i]
    end
    return x
end

function back_solve(A, b)
    Q,R = modified_gram_schmidt(A)
    y = tranpose(Q) * b
    x = back_substitute(R, y)
    return x
end

function classical_gram_schmidt(A)
    m, n = size(A)
    R = zeros(n, n)
    Q = zeros(m, n)
    for j in 1:n
        v_j = A[:, j]
        for i in 1:(j-1)
            R[i, j] = dot(Q[:, i], A[:, j])
            v_j = v_j - R[i, j] * Q[:, i]
        end
        R[j, j] = norm(v_j)
        Q[:, j] = v_j / R[j, j]
    end
    return Q, R
end

function modified_gram_schmidt(A)
    m, n = size(A)
    R = zeros(Float64, n, n)
    Q = zeros(Float64, m, n)

    V = zeros(Float64, n, n)
    for i in 1:n
        V[:, i] = A[:, i]
    end

    for i in 1:n
        R[i, i] = norm(V[:, i])
        Q[:, i] = V[:, i] / R[i, i]
        for j in (i+1):n
            R[i, j] = dot(Q[:, i], V[:, j])
            V[:, j] = V[:, j] - R[i, j] * Q[:, i]
        end
    end
    return Q, R
end

function check_stability(Q)
    m, n = size(Q)
    total_dot_error = 0
    for i in 1:m
        for j in 1:n
            if i != j
                total_dot_error += abs(dot(Q[:, i], Q[:, j]))
            end
        end
    end
    return total_dot_error
end

function boostrap_confint(X, y, bootstrap_samples)
    n, d = size(X)
    beta_arr = zeros(bootstrap_samples, d)

    for i in 1:bootstrap_samples
        sample_arr = sample(1:n,n)
        sampled_X = X[sample_arr, :]
        sampled_y = y[sample_arr]
        beta_arr[i, :] = linear_regression(sampled_X, sampled_y)
    end

    betas = linear_regression(X,y)
    confidence_intervals = zeros(n, 2)

    for i in 1:n
        beta_i = beta_arr[i, :]
        beta_std = std(beta_i)
        confidence_intervals[i, 1] = betas[i] - 1.96 * beta_std
        confidence_intervals[i, 2] = betas[i] + 1.96 * beta_std

    end
    return confidence_intervals
end

# Very Simple Stability Analysis

A = [1 2 3; 4 1 6; 7 8 1]
Q, R = classical_gram_schmidt(A)

classical_error = check_stability(Q)
Q, R = modified_gram_schmidt(A)
modified_error = check_stability(Q)

errror_ratios = classical_error / modified_error
println("classical error $classical_error")
println("modified error $modified_error")
println("Error ratios: $errror_ratios")

max_stars_repo_name: CornellDataScience/interpret-ml

max_issues_repo_name: CornellDataScience/interpret-ml

max_forks_repo_name: CornellDataScience/interpret-ml


text:
function get_attitude_provider(obj::FieldStateMapper)
    return jcall(obj, "getAttitudeProvider", AttitudeProvider, ())
end

function get_frame(obj::FieldStateMapper)
    return jcall(obj, "getFrame", Frame, ())
end

function get_mu(obj::FieldStateMapper)
    return jcall(obj, "getMu", RealFieldElement, ())
end

function get_orbit_type(obj::FieldStateMapper)
    return jcall(obj, "getOrbitType", OrbitType, ())
end

function get_position_angle_type(obj::FieldStateMapper)
    return jcall(obj, "getPositionAngleType", PositionAngle, ())
end

function get_reference_date(obj::FieldStateMapper)
    return jcall(obj, "getReferenceDate", FieldAbsoluteDate, ())
end

function map_array_to_state(obj::FieldStateMapper, arg0::FieldAbsoluteDate, arg1::Vector{RealFieldElement}, arg2::Vector{RealFieldElement}, arg3::PropagationType)
    return jcall(obj, "mapArrayToState", FieldSpacecraftState, (FieldAbsoluteDate, Vector{RealFieldElement}, Vector{RealFieldElement}, PropagationType), arg0, arg1, arg2, arg3)
end

function map_array_to_state(obj::FieldStateMapper, arg0::RealFieldElement, arg1::Vector{RealFieldElement}, arg2::Vector{RealFieldElement}, arg3::PropagationType)
    return jcall(obj, "mapArrayToState", FieldSpacecraftState, (RealFieldElement, Vector{RealFieldElement}, Vector{RealFieldElement}, PropagationType), arg0, arg1, arg2, arg3)
end

function map_date_to_double(obj::FieldStateMapper, arg0::FieldAbsoluteDate)
    return jcall(obj, "mapDateToDouble", RealFieldElement, (FieldAbsoluteDate,), arg0)
end

function map_double_to_date(obj::FieldStateMapper, arg0::RealFieldElement)
    return jcall(obj, "mapDoubleToDate", FieldAbsoluteDate, (RealFieldElement,), arg0)
end

function map_double_to_date(obj::FieldStateMapper, arg0::RealFieldElement, arg1::FieldAbsoluteDate)
    return jcall(obj, "mapDoubleToDate", FieldAbsoluteDate, (RealFieldElement, FieldAbsoluteDate), arg0, arg1)
end

function map_state_to_array(obj::FieldStateMapper, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement}, arg2::Vector{RealFieldElement})
    return jcall(obj, "mapStateToArray", void, (FieldSpacecraftState, Vector{RealFieldElement}, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function set_position_angle_type(obj::FieldStateMapper)
    return jcall(obj, "setPositionAngleType", void, ())
end


text:
@testset "find_alpha.jl" begin

end


text:
# Recreate plot with CMF vs. orbital period:

include("../../../src/regress.jl")
using Statistics
using PyPlot
using DelimitedFiles

period = [1.510826, 2.421937,4.049219, 6.101013, 9.207540,12.352446,18.772866]

dir = "../../../data/POSTERIOR_CMF/"
fname = ["T1b_IRON_FRAC.out","T1c_IRON_FRAC.out","T1d_IRON_FRAC.out","T1e_IRON_FRAC.out","T1f_IRON_FRAC.out","T1g_IRON_FRAC.out","T1h_IRON_FRAC.out"]

np = 7; nsamp = 10000

cmf_norm = zeros(np,nsamp)
cmf_mean = zeros(np)
cmf_median = zeros(np)
cmf_sig = zeros(np)
cmf_sig1 = zeros(np)
cmf_sig2 = zeros(np)
for ip=1:np
  data = readdlm(string(dir,fname[ip]))
  cmf_norm[ip,:] .= data[:,2]
  cmf_mean[ip] = mean(cmf_norm[ip,:])
  cmf_median[ip] = median(cmf_norm[ip,:])
  cmf_sig[ip] = std(cmf_norm[ip,:])
  cmf_sig1[ip] = cmf_median[ip]-sort(cmf_norm[ip,:])[1587]
  cmf_sig2[ip] = sort(cmf_norm[ip,:])[8413]-cmf_median[ip]
end

# Now, carry out regressions for each of these:
coeff_samp = zeros(2,nsamp)
fn = zeros(2,np)
fn[1,:] .= 1.0
fn[2,:] .= period ./6
sig = ones(np)
cmf_avg = 0.0
for isamp=1:nsamp
  global cmf_avg += mean(cmf_norm[:,isamp])
  coeff,cov = regress(fn,cmf_norm[:,isamp],sig)
  coeff_samp[:,isamp] = coeff
end

cmf_avg /= nsamp

# Now, compute standard deviation:
cmf_var = 0.0
for isamp=1:nsamp
  global cmf_var += (mean(cmf_norm[:,isamp]) - cmf_avg)^2
end
cmf_rms = sqrt(cmf_var/nsamp)

# Now, draw 1-sigma and 2-sigma confidence intervals:
period_trial = collect(0.0:0.1:20.0)
conf = zeros(5,length(period_trial))
for (i,p) in enumerate(period_trial)
  # Predict samples at this period
  csamp = coeff_samp[1,:] .+ coeff_samp[2,:] .* p /6
  # Now, sort these and pick out 68.3% & 95% confidence limits:
  csamp = sort(csamp)
  conf[:,i] = [csamp[250],csamp[1415],csamp[5000],csamp[8415],csamp[9750]]
end
# Now, plot median, 1-sigma and 2-sigma confidence intervals:

clf()
#plot(period_trial,conf[1,:],color="C1",alpha=0.5,linestyle="--")
errorbar(period,cmf_mean,yerr=cmf_sig,fmt="o")
#errorbar(period,cmf_median,yerr=[cmf_sig1';cmf_sig2'],fmt="o")
#plot(period_trial,conf[5,:],color="C1",alpha=0.5,linestyle="--")
plot(period_trial,conf[3,:],color="C1",linewidth=2)
plot(period_trial,conf[2,:],color="C1",linewidth=2,linestyle=":")
plot(period_trial,conf[4,:],color="C1",linewidth=2,linestyle=":")
#plot(period_trial,cmf_avg  .+ 0.0 .* period_trial,color="C0",linewidth=2)
fill_between([0,20],[cmf_avg-cmf_rms,cmf_avg-cmf_rms],[cmf_avg+cmf_rms,cmf_avg+cmf_rms],alpha=0.2,color="C0")
xlabel("Orbital Period [d]")
ylabel("Iron mass fraction (%)")
# Overplot solar system bodies:
text([16],[28],"Mars")
plot([0,20],[30.0,30.0],"--",color="k")
text([1],[8],"Moon")
plot([0,20],[10.0,10.0],"--",color="k")
text([18],[33],"Earth")
plot([0,20],[32.0,32.0],"--",color="k")
#text([17],[31.5],"Venus")
#plot([0,20],[31.,31.],"--",color="k")
axis([0,20,0,40])
text([15],[20],"Chondrites (L/CI)")
#plot([0,20],[22.0,22.0],"--",color="k")
plot([0,20],[18.7,18.7],"--",color="k")
text([1],[20.6],"b")
text([2.7],[22.6],"c")
text([3.5],[13],"d")
text([6.3],[20],"e")
text([8.7],[15],"f")
text([11.8],[12],"g")
text([18.1],[11.8],"h")

tight_layout()
savefig("../Norm_cmf_vs_period.pdf",bbox_inches = "tight")

println("Coefficients of fit: ",mean(coeff_samp[1,:]),"+-",std(coeff_samp[1,:])," ",mean(coeff_samp[2,:]),"+-",std(coeff_samp[2,:]))


text:
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.

using DSP #for phase unwrapping

"""
    𝐂 = AMFMdemod(𝚿)
    𝐒 = AMFMdemod([𝚿₀,𝚿₁,𝚿₂])
    𝐒 = AMFMdemod(𝐳)

Create a *numerical canonical triplet* 'numTriplet' from a *numerical component* 'numComp'.

Create a *numerical component set* 'numSet' from a *vector of numerical components*.

Create a *numerical component set* 'numSet' from a *numerical AM--FM model* 'numModel'.

# Examples
```@example
using ISA
ψ₀ = AMFMcomp(t->exp(-t^2),t->2.0,0.0)
fs = 16_000
t = 0:1/fs:1
𝚿₀ = numComp( ψ₀(t), fs )
𝐂₀ = AMFMdemod(𝚿₀)
```

```@example
using ISA
ψ₀ = AMFMcomp(t->exp(-t^2),t->200.0,0.0)
ψ₁ = AMFMcomp(t->1.0,t->100*t,0.1)
ψ₂ = AMFMcomp(t->1+0.8*cos(2t),t->100 + 70.5*sin(t),π)
fs = 16_000
t = 0:1/fs:1
𝚿₀ = numComp( ψ₀(t), fs )
𝚿₁ = numComp( ψ₁(t), fs )
𝚿₂ = numComp( ψ₂(t), fs )
𝐳 = numModel([𝚿₀,𝚿₁,𝚿₂])
𝐒 = AMFMdemod(𝐳)
```
"""
function AMFMdemod(𝚿::numComp; derivMethod="center11")::numTriplet
  return numTriplet( 𝚿.Ψ,
                     𝚿.t,
                     𝚿.fs,
                     abs.(𝚿.Ψ),
                     derivApprox(unwrap(angle.(𝚿.Ψ)),fs=𝚿.fs, method=derivMethod),
                     real.(𝚿.Ψ),
                     imag.(𝚿.Ψ),
                     unwrap(angle.(𝚿.Ψ)),
                     )
end
function AMFMdemod(V::Vector{numComp}; derivMethod="center11")::numSet
  return numSet([ AMFMdemod(𝚿;derivMethod=derivMethod) for 𝚿∈V])
end
function AMFMdemod(𝐳::numModel; derivMethod="center11")::numSet
  return AMFMdemod(𝐳.𝚿ₖ; derivMethod=derivMethod)
end
"""
    f′ = derivApprox(f; fs, method)

Numerically approximate the derivative `f′` of a signal `f` with sampling rate `fs` using one of the a numerical differentiation methods listed below.

`forward` - forward difference

`backward` - backward difference

`center3` - 3-pt stencil central difference

`center5` - 5-pt stencil central difference

`center7` - 7-pt stencil central difference

`center9` - 9-pt stencil central difference

`center11` - 11-pt stencil central difference [default]

`center13` - 13-pt stencil central difference

`center15` - 15-pt stencil central difference

**References**:
[Numerical Differentiation](http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/central-differences/), [Finite Difference Coefficients Calculator](https://web.media.mit.edu/~crtaylor/calculator.html)

Called with a single input `f`, is equivalent to `derivApprox(f, fs=1.0, method="center11")`.

A vector of the same length as the input is returned; values at the signal edges that cannot be computed are filled with `NaN`.

# Examples
```@example
using ISA
x = collect(0.0:0.01:10.0)
f = sin.(x)
f′ = derivApprox(f, fs=1.0, method="center11")
```
"""
function derivApprox(x::Vector{Float64}; fs=1.0, method="center11")::Vector{Float64}
    if method == "forward"
        if length(x)<2; error("derivApprox:vector too short for selected method"); end
        x′ = append!([NaN],fs.*diff(x))
    end
    if method == "backward"
        if length(x)<2; error("derivApprox:vector too short for selected method"); end
        x′ = append!(fs.*diff(x),[NaN])
    end
    if method == "center3"
        if length(x)<3; error("derivApprox:vector too short for selected method"); end
        x′ = append!( append!([NaN], fs.*(x[3:end]-x[1:end-2])./2.0 ) ,[NaN])
    end
    if method == "center5"
        if length(x)<5; error("derivApprox:vector too short for selected method"); end
        x′ = append!( append!(NaN*ones(Float64,2), fs.*(x[1:end-4]-8.0.*x[2:end-3]+8.0.*x[4:end-1]-x[5:end])./12.0 ) ,NaN*ones(Float64,2))
    end
    if method == "center7"
        if length(x)<7; error("derivApprox:vector too short for selected method"); end
        x′ = append!( append!(NaN*ones(Float64,3), fs.*(-x[1:end-6]+ 9.0.*x[2:end-5]- 45.0.*x[3:end-4] +45.0.*x[5:end-2] -9.0.*x[6:end-1] + x[7:end]   )./60.0
 ) ,NaN*ones(Float64,3))
    end
    if method == "center9"
        if length(x)<9; error("derivApprox:vector too short for selected method"); end
        x′ = append!( append!(NaN*ones(Float64,4), fs.*(3.0.*x[1:end-8] - 32.0.*x[2:end-7] + 168.0.*x[3:end-6] -672.0.*x[4:end-5] + 672.0.*x[6:end-3] -168.0.*x[7:end-2] + 32.0.*x[8:end-1] -3.0.*x[9:end]  )./840.0
 ) ,NaN*ones(Float64,4))
    end
    if method == "center11"
        if length(x)<11; error("derivApprox:vector too short for selected method"); end
        x′ = append!( append!(NaN*ones(Float64,5), fs.*(-2.0.*x[1:end-10] +25.0.*x[2:end-9] -150.0.*x[3:end-8] +600.0.*x[4:end-7] -2100.0.*x[5:end-6] + 2100.0.*x[7:end-4]  -600.0.*x[8:end-3] +150.0.*x[9:end-2] -25.0.*x[10:end-1] +2.0.*x[11:end]  )./2520.0
 ) ,NaN*ones(Float64,5))
    end
    if method == "center13"
        if length(x)<13; error("derivApprox:vector too short for selected method"); end
        x′ = append!( append!(NaN*ones(Float64,6), fs.*(5.0.*x[1:end-12] - 72.0.*x[2:end-11] + 495.0.*x[3:end-10] -2200.0.*x[4:end-9]   +7425.0.*x[5:end-8]  -23760.0.*x[6:end-7]  +23760.0.*x[8:end-5]   + -7425.0.*x[9:end-4]   + 2200.0.*x[10:end-3] -495.0.*x[11:end-2] + 72.0.*x[12:end-1] -5.0.*x[13:end]  )./27720.0
 ) ,NaN*ones(Float64,6))
    end
    return x′
end


text:
"""

3D represtation of the minimum distance distribution function

"""

using M3GTools ; const m3g = M3GTools

mysim = m3g.init(psf="../gmd_solute/gmdfiles/structure.psf",
                 pdb="../gmd_solute/gmdfiles/structure.pdb",
                 vmd="vmd")

# Selection that defines what the solute is

solute = m3g.select(mysim,"protein")

# File containting the contribution of the solute atoms to the gmd:

gmd_solute_file = "../gmd_solute/gmdfiles/gmd-GMD_ATOM_SOLUTE_CONTRIB.dat"

# Write PDB file containing the 3D respresentation of the GMD:

output = "./gmd3D.pdb"
gmd3D_write(mysim,solute,gmd_solute_file,output,)


text:
using VLConstraintBasedModelGenerationUtilities

# setup path to sequence file -
path_to_gene_sequence_file = "/Users/jeffreyvarner/Desktop/julia_work/VLConstraintBasedModelGenerationUtilities.jl/test/data/G-MZ373340.fasta"

# load -
gene_table = build_gene_table(path_to_gene_sequence_file) |> check

# get the sequence -
gene_seq = gene_table[!,:gene_sequence][1]

# table -
transcription_table = build_transcription_reaction_table("test_gene", gene_seq) |> check


text:
struct EventValues
    events::Vector{Event}
    vals::Vector{Counts}
    time::Counts
end

EventValues(events::Event...) = EventValues(events, zeros(Counts, length(events)), Counts(0))
EventValues(events::NTuple{N, Event}) where N = EventValues(collect(events), zeros(Counts, N), Counts(0))
EventValues(events::Vector{Event}) = EventValues(events, zeros(Counts, length(events)), Counts(0))

struct EventStats
    events::Vector{Event}
    samples::Matrix{Counts}
    time::Vector{Counts}
end

EventStats(events::Vector{Event}) = EventStats(events, zeros(Counts, (0, length(events))),zeros(Counts, 0))

gcscrub() = (GC.gc(); GC.gc(); GC.gc(); GC.gc())

function run_and_collect!(f::Function, evtsets::Vector{EventSet}, vals::Vector{Counts})
    last = 0

    time = UInt64(0)
    for evtset in evtsets
        first = last + 1
        last = last + length(evtset)

        start_counters(evtset)
        time += try
            local t0 = time_ns()
            f()
            (time_ns() - t0)
        finally
            stop_counters!(evtset, view(vals, first:last))
        end
    end
    div(time, length(evtsets))
end

"""
    profile(f::Function, events::Vector{Event}; gcfirst::Bool=true, warmup::Int64=0)

Execute the function `f` once while counting specific `events`.

**Arguments**:

    -`f`: the function to profile
    -`events`: the events to count
    -`gcfirst`: run the gc several times before the execution to reduce gc noise
    -`warmup`: number of times to run the function prior to counting

**Return values**:

`EventValues` containing the events, counts and runtime collected
"""
function profile(f::Function, events::Vector{Event}; gcfirst::Bool=true, warmup::Int64=0, multi::Bool=true)
    gcfirst && gcscrub()

    for i in 1:warmup
      f()
    end

    evtsets = measurement_groups(events)
    if length(evtsets) != 1 && !multi
        throw(ArgumentError("cannot profile the events at the same time"))
    end

    vals = zeros(Counts, length(events))
    time = run_and_collect!(f, evtsets, vals)
    EventValues(collect_events(evtsets), vals, time)
end

profile(f::Function, events::NTuple{N, Event}) where N = profile(f, collect(events))

function _sample(f::Function, evtsets::Vector{EventSet}, counts::Vector{Counts}, samples::Vector{Vector{Counts}}, times::Vector{UInt64};
         max_secs::Real=5., max_epochs::Int64=1000, gcsample::Bool=false)

    start_time = Base.time()
    iters = 1
    while (Base.time() - start_time) < max_secs && iters ≤ max_epochs
        gcsample && gcscrub()

        time = run_and_collect!(f, evtsets, counts)

        push!(samples, copy(counts))
        push!(times, time)
        iters += 1
    end
end

function _sample(f::Function, evtset::EventSet, counts::Vector{Counts}, samples::Vector{Vector{Counts}}, times::Vector{UInt64};
         max_secs::Real=5., max_epochs::Int64=1000, gcsample::Bool=false)

    start_counters(evtset)
    try
        start_time = Base.time()
        iters = 1
        while (Base.time() - start_time) < max_secs && iters ≤ max_epochs
            gcsample && gcscrub()

            reset_counters!(evtset)
            local t0 = time_ns()
            f()
            time = (time_ns() - t0)
            read_counters!(evtset, counts)

            push!(samples, copy(counts))
            push!(times, time)
            iters += 1
        end
    finally
        stop_counters(evtset)
    end
end

"""
    sample(f::Function, events::Vector{Event}; max_secs::Float64=5, max_epochs::Int64=1000, gcsample::Bool=false, warmup::Int64=1)

Execute the function `f` several times, each time counting specific `events`.
Sampling continues until either the maximum number of samples `max_epochs` are collected or the runtime budget `max_secs` is exceeded.

**Arguments**:

    -`f`: the function to profile
    -`events`: the events to count
    -`max_secs`: maximum number of seconds to sample for
    -`gcsample`: run the gc several times before the execution to reduce gc noise
    -`warmup`: number of times to run the function prior to counting

**Return values**:

`EventStats` containing the events, counts and runtime collected
"""
function sample(f::Function, events::Vector{Event}; max_secs::Real=5., max_epochs::Int64=1000, gcsample::Bool=false, warmup::Int64=1, multi::Bool=true)
    num_events = length(events)
    counts = Vector{Counts}(undef, num_events)
    samples = Vector{Counts}[]
    times = UInt64[]

    evtsets = measurement_groups(events)
    if length(evtsets) != 1 && !multi
        throw(ArgumentError("cannot profile the events at the same time"))
    end

    gcscrub()

    for i in 1:warmup
        f()
    end

    if length(evtsets) == 1
        _sample(f, first(evtsets), counts, samples, times; max_secs=max_secs, max_epochs=max_epochs, gcsample=gcsample)
    else
        _sample(f, evtsets, counts, samples, times; max_secs=max_secs, max_epochs=max_epochs, gcsample=gcsample)
    end
    EventStats(collect_events(evtsets), hcat(samples...)', times)
end

sample(f::Function, events::NTuple{N, Event}; kw...) where N = sample(f, collect(events); kw...)

function kwargs(default_events, ex, args...)
    events, ex, params = if isa(ex, Symbol) || (isa(ex, Expr) && (ex.head == :tuple || ex.head == :vect))
      ex, first(args), collect(Iterators.drop(args, 1))
    else
      default_events, ex, collect(args)
    end

    for ex in params
        if isa(ex, Expr) && ex.head == :(=)
            ex.head = :kw
        end
    end
    events, ex, params
end

"""
    profile(ex, args...)

Convience macro for profiling an expression.
Events can be specified as a first argument, otherwise the default events `[BR_INS, BR_MSP, TOT_INS, TOT_CYC]` are counted

Arguments and return values are similar to [`profile`](@ref)

**Example**:
```julia
@profile f(x, y, z) # sampling default events
@profile [PAPI.TOT_INS, PAPI.DP_OPS, native_event] f(x, y, z) gcfirst=false
```
"""
macro profile(ex, args...)
    events, ex, params = kwargs([BR_INS, BR_MSP, TOT_INS, TOT_CYC], ex, args...)
    quote
        profile(() -> $(esc(ex)), Event[$(esc(events))...], $(params...))
    end
end

"""
    sample(ex, args...)

Convience macro for sampling an expression.
Events can be specified as a first argument, otherwise the default events `[BR_INS, BR_MSP, TOT_INS, TOT_CYC]` are counted

Arguments and return values are similar to [`sample`](@ref)

**Example**:
```julia
@sample f(x, y, z) # sampling default events
@sample [PAPI.TOT_INS, PAPI.DP_OPS, native_event] f(x, y, z) max_secs=1
```
"""
macro sample(ex, args...)
    events, ex, params = kwargs([BR_INS, BR_MSP, TOT_INS, TOT_CYC], ex, args...)
    quote
        sample(() -> $(esc(ex)), Event[$(esc(events))...], $(params...))
    end
end


text:
# by lazarusA
# using GLMakie # HIDE
using Legendre, GLMakie
using GeometryBasics, LinearAlgebra, StatsBase
using Makie: get_dim, surface_normals
GLMakie.activate!() # HIDE
# thanks to @jkrumbiegel for the lift
let
    function Y(θ, ϕ, l, m)
        if m < 0
            return (-1)^m * √2 * Nlm(l, abs(m)) * Plm(l, abs(m), cos(θ)) * sin(abs(m)*ϕ)
        elseif m == 0
            return sqrt((2*l+1)/4π)*Plm(l, m, cos(θ))
        else
            return (-1)^m * √2 * Nlm(l, m) * Plm(l, m, cos(θ)) * cos(m*ϕ)
        end
    end
    function getMesh(x,y,z)
        positions = vec(map(CartesianIndices(z)) do i
        GeometryBasics.Point{3, Float32}(
            get_dim(x, i, 1, size(z)),
            get_dim(y, i, 2, size(z)),
            z[i])
        end)
        faces = decompose(GLTriangleFace, Rect2D(0f0, 0f0, 1f0, 1f0), size(z))
        normals = surface_normals(x, y, z)
        vertices = GeometryBasics.meta(positions; normals=normals)
        meshObj = GeometryBasics.Mesh(vertices, faces)
        meshObj
    end
    # Grids of polar and azimuthal angles
    θ = LinRange(0, π, 200)
    ϕ = LinRange(0, 2π, 200)
    x = [sin(θ)*sin(ϕ) for θ in θ, ϕ in ϕ]
    y = [sin(θ)*cos(ϕ) for θ in θ, ϕ in ϕ]
    z = [cos(θ)        for θ in θ, ϕ in ϕ]

    l = Node(4)
    m = Node(1)

    ambient =  Vec3f0(0.75, 0.75, 0.75)
    cmap = (:dodgerblue, :white) # how to include this into menu options?
    with_theme(theme_black()) do
        fig = Figure(resolution = (1400, 800))
        menu = Menu(fig, options = ["Spectral_11", "viridis", "heat", "plasma", "magma", "inferno"])
        Ygrid = lift(l, m) do l, m
            [Y(θ, ϕ, l, m) for θ in θ, ϕ in ϕ]
        end
        Ylm = @lift(abs.($Ygrid))
        Ygrid2 = @lift(vec($Ygrid))

        ax1 = Axis3(fig, aspect = :data, perspectiveness = 0.5, elevation = π/8, azimuth = 2.225π)
        ax2 = Axis3(fig, aspect = :data, perspectiveness = 0.5, elevation = π/8, azimuth = 2.225π)
        pltobj1 = mesh!(ax1, getMesh(x, y, z), color = Ygrid2, colormap = cmap, ambient = ambient)
        pltobj2 = mesh!(ax2, @lift(getMesh($Ylm .* x, $Ylm .* y, $Ylm .* z)), color = Ygrid2,
            colormap = cmap,  ambient = ambient, limits = FRect3D((-0.3,-0.3,-0.6),(0.6,0.6,1.2)))
        #hidespines!(ax1)
        #hidedecorations!(ax1)
        #hidespines!(ax2)
        #hidedecorations!(ax2)
        cbar = Colorbar(fig, pltobj1, label = "Yₗₘ(θ,ϕ)", width = 11, tickalign = 1, tickwidth = 1)
        fig[1,1] = ax1
        fig[1,2] = ax2
        fig[1,3] = cbar
        fig[0,1:2] = Label(fig, @lift("Tesseral Spherical Harmonics l = $($l), m = $($m)"), textsize = 20)
        fig[2, 0] = vgrid!(
            Label(fig, "Colormap", width = nothing),
            menu; tellheight = false, width = 100)
        on(menu.selection) do s
            pltobj1.colormap = s
            pltobj2.colormap = s
        end
        sl = Slider(fig[end+1, 2:3], range = 1:30, startvalue = 2)
        sl2 = Slider(fig[end+1, 2:3], range = @lift(-$(sl.value):1:$(sl.value)))
        connect!(l, sl.value)
        connect!(m, sl2.value)
        tight_ticklabel_spacing!(cbar)
        fig
        display(fig) # HIDE
        save(joinpath(@__DIR__, "output", "tesseralSphericalHSliders.png"), fig, px_per_unit = 2.0) # HIDE
    end
end
using Pkg # HIDE
Pkg.status(["GLMakie", "Legendre","GeometryBasics", "LinearAlgebra", "StatsBase", "Makie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Cairo, Images, Colors

abstract type AbstractCtx end
abstract type AbstractDrawParams end

struct CairoCtx <: AbstractCtx
    params::AbstractDrawParams
    viewer::Cairo.CairoSurfaceBase{UInt32}
end

struct RGBCtx <: AbstractCtx
    params::AbstractDrawParams
    viewer::Cairo.CairoSurfaceBase{UInt32}
end

function render!(env::AbstractEnv, ctx::CairoCtx)
    drawcanvas!(env, CairoContext(ctx.viewer), ctx.params)
    display(ctx.viewer)
end

function render!(env::AbstractEnv, ctx::RGBCtx)
    drawcanvas!(env, CairoContext(ctx.viewer), ctx.params)
    ptr = ccall((:cairo_image_surface_get_data, Cairo.libcairo), Ptr{UInt32}, (Ptr{Nothing},), ctx.viewer.ptr)
    arr = unsafe_wrap(Array, ptr, (ctx.params.screen_width, ctx.params.screen_height))
    rgb_arr = convert.(Float64, channelview(colorview(RGB{N0f8}, permutedims(reinterpret(RGB24, arr), [2, 1]))))
    return rgb_arr
end

@init @require Gtk="4c0ca9eb-093a-5379-98c5-f87ac0bbbf44" begin
	using Gtk

	struct GtkCtx <: AbstractCtx
		params::AbstractDrawParams
		canvas::GtkCanvas
		win::GtkWindowLeaf
	end

	function render!(env::AbstractEnv, ctx::GtkCtx)
		!visible(ctx.win) && visible(ctx.win, true)
		@guarded draw(ctx.canvas) do widget
		    drawcanvas!(env, getgc(ctx.canvas), ctx.params)
		end
		reveal(ctx.canvas, true)
		return ctx.canvas
	end
end


text:
function pull_request_build(::NewPackage,
                            pr::GitHub.PullRequest,
                            current_pr_head_commit_sha::String,
                            registry::GitHub.Repo;
                            auth::GitHub.Authorization,
                            authorized_authors::Vector{String},
                            authorized_authors_special_jll_exceptions::Vector{String},
                            registry_head::String,
                            registry_master::String,
                            suggest_onepointzero::Bool,
                            whoami::String,
                            registry_deps::Vector{<:AbstractString} = String[])::Nothing
    # first check if the PR is open, and the author is authorized - if not, then quit
    # if the PR is open and the author is authorized, then check rules 0 through 10.
    # Rules:
    # 0. A JLL-only author (e.g. `jlbuild`) is not allowed to register non-JLL packages.
    # 1. Only changes a subset of the following files:
    #     - `Registry.toml`,
    #     - `E/Example/Compat.toml`
    #     - `E/Example/Deps.toml`
    #     - `E/Example/Package.toml`
    #     - `E/Example/Versions.toml`
    # 2. TODO: implement this check. When implemented, this check will make sure that the changes to `Registry.toml` only modify the specified package.
    # 3. Normal capitalization
    #     - name should match r"^[A-Z]\w*[a-z]\w*[0-9]?$"
    #     - i.e. starts with a capital letter, ASCII alphanumerics only, contains at least 1 lowercase letter
    # 4. Not too short
    #     - at least five letters
    #     - you can register names shorter than this, but doing so requires someone to approve
    # 5. Standard initial version number - one of 0.0.1, 0.1.0, 1.0.0
    #     - does not apply to JLL packages
    # 6. DISABLED. Repo URL ends with /$name.jl.git where name is the package name. Now that we have support for multiple packages in different subdirectories of a repo, we have disabled this check.
    # 7. Compat for all dependencies
    #     - there should be a [compat] entry for Julia
    #     - all [deps] should also have [compat] entries
    #     - all [compat] entries should have upper bounds
    #     - dependencies that are standard libraries do not need [compat] entries
    #     - dependencies that are JLL packages do not need [compat] entries
    # 8. (only applies to JLL packages) The only dependencies of the package are:
    #     - Pkg
    #     - Libdl
    #     - other JLL packages
    # 9. Version can be installed
    #     - given the proposed changes to the registry, can we resolve and install the new version of the package?
    #     - i.e. can we run `Pkg.add("Foo")`
    # 10. Version can be loaded
    #     - once it's been installed (and built?), can we load the code?
    #     - i.e. can we run `import Foo`
    pkg, version = parse_pull_request_title(NewPackage(), pr)
    this_is_jll_package = is_jll_name(pkg)
    @info("This is a new package pull request", pkg, version, this_is_jll_package)
    pr_author_login = author_login(pr)
    if is_open(pr)
        if pr_author_login in vcat(authorized_authors, authorized_authors_special_jll_exceptions)
            description = "New package. Pending."
            params = Dict("state" => "pending",
                          "context" => "automerge/decision",
                          "description" => description)
            my_retry(() -> GitHub.create_status(registry,
                                                current_pr_head_commit_sha;
                                                auth = auth,
                                                params = params))

            if this_is_jll_package
                if pr_author_login in authorized_authors_special_jll_exceptions
                    this_pr_can_use_special_jll_exceptions = true
                else
                    this_pr_can_use_special_jll_exceptions = false
                end
            else
                this_pr_can_use_special_jll_exceptions = false
            end

            if this_is_jll_package
                g0 = true
                m0 = ""
            else
                if pr_author_login in authorized_authors
                    g0 = true
                    m0 = ""
                else
                    g0 = false
                    m0 = "This package is not a JLL package. The author of this pull request is not authorized to register non-JLL packages."
                end
            end

            g1, m1 = pr_only_changes_allowed_files(NewPackage(),
                                                   registry,
                                                   pr,
                                                   pkg;
                                                   auth = auth)
            g2 = true
            m2 = ""
            if this_pr_can_use_special_jll_exceptions
                g3 = true
                g4 = true
                m3 = ""
                m4 = ""
            else
                g3, m3 = meets_normal_capitalization(pkg)
                g4, m4 = meets_name_length(pkg)
            end
            if this_pr_can_use_special_jll_exceptions
                g5 = true
                m5 = ""
            else
                g5, m5 = meets_standard_initial_version_number(version)
            end

            # g6, m6 = meets_repo_url_requirement(pkg; registry_head = registry_head)
            g6 = true
            m6 = ""

            g7, m7 = meets_compat_for_all_deps(registry_head,
                                               pkg,
                                               version)
            g8_if_jll, m8_if_jll = meets_allowed_jll_nonrecursive_dependencies(registry_head,
                                                                               pkg,
                                                                               version)
            if this_is_jll_package
                g8 = g8_if_jll
                m8 = m8_if_jll
            else
                g8 = true
                m8 = ""
            end

            @info("JLL-only authors cannot register non-JLL packages.",
                  meets_this_guideline = g0,
                  message = m0)
            @info("Only modifies the files that it's allowed to modify",
                  meets_this_guideline = g1,
                  message = m1)
            @info("TODO: implement this check",
                  meets_this_guideline = g2,
                  message = m2)
            @info("Normal capitalization",
                  meets_this_guideline = g3,
                  message = m3)
            @info("Name not too short",
                  meets_this_guideline = g4,
                  message = m4)
            @info("Standard initial version number ",
                  meets_this_guideline = g5,
                  message = m5)
            @info("Repo URL ends with /name.jl.git",
                  meets_this_guideline = g6,
                  message = m6)
            @info("Compat (with upper bound) for all dependencies",
                  meets_this_guideline = g7,
                  message = m7)
            @info("If this is a JLL package, only deps are Pkg, Libdl, and other JLL packages",
                  meets_this_guideline = g8,
                  message = m8)
            g0through8 = Bool[g0,
                              g1,
                              g2,
                              g3,
                              g4,
                              g5,
                              g6,
                              g7,
                              g8]
            if !all(g0through8)
                description = "New package. Failed."
                params = Dict("state" => "failure",
                              "context" => "automerge/decision",
                              "description" => description)
                my_retry(() -> GitHub.create_status(registry,
                                                    current_pr_head_commit_sha;
                                                    auth = auth,
                                                    params = params))
            end
            g9, m9 = meets_version_can_be_pkg_added(registry_head,
                                                    pkg,
                                                    version;
                                                    registry_deps = registry_deps)
            @info("Version can be `Pkg.add`ed",
                  meets_this_guideline = g9,
                  message = m9)
            g10, m10 = meets_version_can_be_imported(registry_head,
                                                     pkg,
                                                     version;
                                                     registry_deps = registry_deps)
            @info("Version can be `import`ed",
                  meets_this_guideline = g9,
                  message = m9)
            g0through10 = Bool[g0,
                               g1,
                               g2,
                               g3,
                               g4,
                               g5,
                               g6,
                               g7,
                               g8,
                               g9,
                               g10]
            allmessages0through10 = String[m0,
                                           m1,
                                           m2,
                                           m3,
                                           m4,
                                           m5,
                                           m6,
                                           m7,
                                           m8,
                                           m9,
                                           m10]
            if all(g0through10) # success
                description = "New package. Approved. name=\"$(pkg)\". sha=\"$(current_pr_head_commit_sha)\""
                params = Dict("state" => "success",
                              "context" => "automerge/decision",
                              "description" => description)
                my_retry(() -> GitHub.create_status(registry,
                                                    current_pr_head_commit_sha;
                                                    auth = auth,
                                                    params = params))
                this_pr_comment_pass = comment_text_pass(NewPackage(),
                                                         suggest_onepointzero,
                                                         version)
                my_retry(() -> update_automerge_comment!(registry,
                                                         pr;
                                                         auth = auth,
                                                         body = this_pr_comment_pass,
                                                         whoami = whoami))
                return nothing
            else # failure
                description = "New package. Failed."
                params = Dict("state" => "failure",
                              "context" => "automerge/decision",
                              "description" => description)
                my_retry(() -> GitHub.create_status(registry,
                                                    current_pr_head_commit_sha;
                                                    auth = auth,
                                                    params = params))
                failingmessages0through10 = allmessages0through10[.!g0through10]
                this_pr_comment_fail = comment_text_fail(NewPackage(),
                                                         failingmessages0through10,
                                                         suggest_onepointzero,
                                                         version)
                my_retry(() -> update_automerge_comment!(registry,
                                                         pr;
                                                         body = this_pr_comment_fail,
                                                         auth = auth,
                                                         whoami = whoami))
                throw(AutoMergeGuidelinesNotMet("The automerge guidelines were not met."))
            end
        else
            throw(AutoMergeAuthorNotAuthorized("Author $(pr_author_login) is not authorized to automerge. Exiting..."))
        end
    else
        throw(AutoMergePullRequestNotOpen("The pull request is not open. Exiting..."))
    end
end