text:
# loadMetropolis.jl
# Julia 0.4 code to generate metropolis
# Authored by Arthur J Delarue

MAX_VELOCITIES = [1.,1.,1.]#[0.1,0.3,0.5]
MEAN_VELOCITIES = [0.03,0.15,0.3]
OUTSIDE_NODES = [(1,1),(1,2),(1,3),(1,4),(1,5),(1,6),(1,7),(1,8),(2,8),(3,8),(4,8),(5,8),(6,8),(7,8),(8,8),(8,7),(8,6),(8,5),(8,4),(8,3),(8,2),(8,1),(7,1),(6,1),(5,1),(4,1),(3,1),(2,1),(1,1)]
# OUTSIDE_NODES = [(1,1),(1,2),(1,3),(2,3),(3,3),(3,2),(3,1),(2,1),(1,1)]
WIDTH = 8
NSUB = 8
TURN_COST = 2.0

using TaxiSimulation
using Distributions
using HDF5, JLD

#Gives the index of vertex, given coords and city number
function coordToLoc(i::Int, j::Int, c::Int, width::Int, nSub::Int)
	subWidth = floor(Int, width/2)
    if c ==0
        return j + (i-1)*width
    else
        return width^2 + (c-1)*subWidth^2 + j + (i-1)*subWidth
    end
end

function load_metropolis_graph(;from_scratch=false)
	path = "/Users/arthurdelarue/.julia/v0.4/TaxiSimulation"
	if from_scratch
		# Create city
		city = Metropolis(WIDTH, NSUB)

		function drawNetworkToText(pb::TaxiProblem, name::String = "graph")
		 	stdin, proc = open(`neato -Tplain -o $(path)/outputs/$(name).txt`, "w")
		 	TaxiSimulation.to_dot(pb,stdin)
		 	close(stdin)
		end

		# GraphViz magic to find coordinates for nodes
		drawNetwork(city)
		GraphViz = Coordinates[]
		indices = Int64[]
		drawNetworkToText(city, "test1")
		fileExists = false
		while (!fileExists)
			sleep(1)
			fileExists = isfile("$(path)/outputs/test1.txt")
		end
		lines = readlines(open("$(path)/outputs/test1.txt"))
		rm("$(path)/outputs/test1.txt")
		index = 2
		while(split(lines[index])[1] == "node")
			push!(indices, convert(Int64, float(split(lines[index])[2])))
			nodeC = Coordinates(float(split(lines[index])[3]), float(split(lines[index])[4]))
			push!(GraphViz, nodeC)
			index += 1
		end

		GraphVizNodes = copy(GraphViz)
		for i = 1:length(GraphViz)
			GraphVizNodes[indices[i]] = GraphViz[i]
		end

		originalNodes = GraphVizNodes
		distances = spzeros(TaxiSimulation.nv(city.network),TaxiSimulation.nv(city.network))

		for i in TaxiSimulation.vertices(city.network), j in TaxiSimulation.out_neighbors(city.network, i)
			distances[i,j] = norm([originalNodes[i].x - originalNodes[j].x, originalNodes[i].y - originalNodes[j].y])
		end

		adjacencyList = Array{Int}[]
		coordinates = zeros(TaxiSimulation.nv(city.network), 2)
		for i = 1:TaxiSimulation.nv(city.network)
			push!(adjacencyList, TaxiSimulation.out_neighbors(city.network, i))
			coordinates[i,1] = originalNodes[i].x
			coordinates[i,2] = originalNodes[i].y
		end

		graph = city.network
		positions = originalNodes

		save("Inputs/input-graph.jld", "adjList", adjacencyList)
		save("Inputs/input-positions.jld", "coordinates", coordinates)
		save("Inputs/input-distances.jld", "distances", distances)
	else
		# load graph
		adjacencyList = load("Inputs/input-graph.jld", "adjList")
		graph = TaxiSimulation.DiGraph(length(adjacencyList))
		for i = 1:TaxiSimulation.nv(graph), j in adjacencyList[i]
			TaxiSimulation.add_edge!(graph, i, j)
		end
		# load coordinates
		coordinates = load("Inputs/input-positions.jld", "coordinates")
		positions = TaxiSimulation.Coordinates[]
		for i = 1:TaxiSimulation.nv(graph)
			push!(positions, TaxiSimulation.Coordinates(coordinates[i,1], coordinates[i,2]))
		end
		# load distances
		distances = load("Inputs/input-distances.jld", "distances")
	end
	return graph, positions, distances
end

function define_travel_times(;from_scratch = false, turn_cost::Float64 = TURN_COST)
	input_speed = 0
	function gen_speed(index::Int)
		input_speed += 0.001 #MEAN_VELOCITIES[index]
		speed = -1
		while speed < 0 || speed > 9.9#MAX_VELOCITIES[index]
			speed = rand(Normal(input_speed, input_speed/3))#15))
		end
		return input_speed
	end
	if from_scratch
		graph, positions, distances = load_metropolis_graph(from_scratch=false)

		minRoadTime = spzeros(TaxiSimulation.nv(graph),TaxiSimulation.nv(graph))
		meanRoadTime = spzeros(TaxiSimulation.nv(graph),TaxiSimulation.nv(graph))
		# Set normal roads
		for i in TaxiSimulation.vertices(graph), j in TaxiSimulation.out_neighbors(graph, i)
			minRoadTime[i,j] = distances[i,j]/MAX_VELOCITIES[2]
			meanRoadTime[i,j] = distances[i,j]/gen_speed(2)
		end
		# Set fast roads on outside of main city
		for k = 1:(length(OUTSIDE_NODES)-1)
			i1,j1 = OUTSIDE_NODES[k]
			node1 = coordToLoc(i1,j1,0,WIDTH,NSUB)
			i2,j2 = OUTSIDE_NODES[k+1]
			node2 = coordToLoc(i2,j2,0,WIDTH,NSUB)
			minRoadTime[node1,node2] = distances[node1,node2]/MAX_VELOCITIES[3]
			minRoadTime[node2,node1] = distances[node2,node1]/MAX_VELOCITIES[3]
			meanRoadTime[node1,node2] = distances[node1,node2]/gen_speed(3)
			meanRoadTime[node2,node1] = distances[node2,node1]/gen_speed(3)
		end
		# Set slow roads between suburbs and city
		for k = 1:NSUB
			suburbNode = coordToLoc(1,1,k,WIDTH,NSUB)
			cityNode = minimum(TaxiSimulation.out_neighbors(graph, suburbNode))
			minRoadTime[suburbNode, cityNode] = distances[suburbNode, cityNode]/MAX_VELOCITIES[1]
			minRoadTime[cityNode, suburbNode] = distances[cityNode, suburbNode]/MAX_VELOCITIES[1]
			meanRoadTime[suburbNode, cityNode] = distances[suburbNode, cityNode]/gen_speed(1)
			meanRoadTime[cityNode, suburbNode] = distances[cityNode, suburbNode]/gen_speed(1)
		end

		cityPaths = realPaths(graph, meanRoadTime, meanRoadTime, positions, turn_cost, turn_cost)
		traveltime = cityPaths.traveltime

		save("Inputs/input-realTimes.jld", "realTimes", cityPaths.traveltime)
		save("Inputs/input-speedLimits.jld", "speedLimits", minRoadTime)
		save("Inputs/input-meanTimes.jld", "meanTimes", meanRoadTime)
	else
		traveltime = load("Inputs/input-realTimes.jld", "realTimes")
		minRoadTime = load("Inputs/input-speedLimits.jld", "speedLimits")
		meanRoadTime = load("Inputs/input-meanTimes.jld", "meanTimes")
	end
	return traveltime, minRoadTime, meanRoadTime
end

function generate_rides(prob = 0.7)
	graph, positions, distances = load_metropolis_graph(from_scratch=false)
	traveltime, minRoadTime, meanRoadTime = define_travel_times(from_scratch=false)

	out = [copy(TaxiSimulation.out_neighbors(graph, i)) for i = TaxiSimulation.vertices(graph)]
	out2 = [vcat([copy(TaxiSimulation.out_neighbors(graph, j)) for j in out[i]]) for i = TaxiSimulation.vertices(graph)]

	cityToCityDist = Categorical([prob, (1-prob)/6, (1-prob)/6, (1-prob)/6, (1-prob)/6, (1-prob)/6, (1-prob)/6])
	cityToSuburbDist = Categorical([prob, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10])
	suburbToCityDist = Categorical([prob, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10, (1-prob)/10])
	suburbToSuburbDist = Categorical([prob, (1-prob)/4, (1-prob)/4, (1-prob)/4, (1-prob)/4])

	numRides = zeros(Int, (TaxiSimulation.nv(graph), TaxiSimulation.nv(graph)))
	travelTimes = zeros(TaxiSimulation.nv(graph), TaxiSimulation.nv(graph))
	for i = 1:TaxiSimulation.nv(graph), j = 1:TaxiSimulation.nv(graph)
		if i != j && !(j in out[i]) && !(j in out2[i])
			# city to city
			if i < 65 && j < 65
				numRides[i,j] = rand(cityToCityDist) - 1
			# city to suburb
			elseif i < 65 && j >= 65
				numRides[i,j] = rand(cityToSuburbDist) - 1
			# suburb to city
			elseif i >= 65 && j < 65
				numRides[i,j] = rand(suburbToCityDist) - 1
			# suburb to suburb
			else
				numRides[i,j] = rand(suburbToSuburbDist) - 1
			end
			if numRides[i,j] > 0
				travelTimes[i,j] = mean(rand(Normal(traveltime[i,j], sqrt(traveltime[i,j])), numRides[i,j]))
				while travelTimes[i,j] < 0
					travelTimes[i,j] = mean(rand(Normal(traveltime[i,j], sqrt(traveltime[i,j])), numRides[i,j]))
				end
			end
		end
	end

	save("Inputs/input-travelTimes-$(prob).jld", "travelTimes", travelTimes)
	save("Inputs/input-numRides-$(prob).jld", "numRides", numRides)
	return travelTimes, numRides
end

function create_metropolis(graph::TaxiSimulation.DiGraph, roadTime::AbstractArray{Float64,2})
	city = Metropolis(WIDTH, NSUB)
	city.network = graph
	city.roadTime = roadTime
	return city
end

function updatedMetropolis(;prob::Float64=0.7, dir = "", max_rounds = 20, num_iter = -1)
	graph, positions, distances = load_metropolis_graph()
	if dir == ""
		path = "Outputs/minr1_i$(max_rounds)_metropolis_$(prob)"
	else
		path = "Outputs/$dir"
	end
	if num_iter < 0
		if isfile("$(path)/end.jld")
			num_iter = load("$(path)/end.jld", "num_iter")
		else
			num_iter = max_rounds
		end
	end
	roadTime = load("$(path)/metropolis-times-$(num_iter).jld", "times")
	return create_metropolis(graph, roadTime)
end

function showRealMetropolis()
	g, p, d = load_metropolis_graph()
	tt, mt, at = define_travel_times()
	visualize(create_metropolis(g, at))
end

function showStartingMetropolis()
	g, p, d = load_metropolis_graph()
	tt, mt, at = define_travel_times()
	visualize(create_metropolis(g, mt))
end

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
# combination of prep_save.jl and align_prep.jl.
# Loads from NMR experiment instead of intermediate BSON file.

import NMRDataSetup
import NMRSpectraSimulator
import NMRSpecifyRegions

include("../src/NMRCalibrate.jl")
import .NMRCalibrate
#import NMRCalibrate

using LinearAlgebra
using FFTW

import PyPlot
#import PlotlyJS
#using Plots; plotly()

import BSON
import JSON

import Statistics
import Random
Random.seed!(25)

##### global constants.
#cs_config_path = "/home/roy/Documents/repo/NMRData/src/input/reduced_cs_config.txt"
SH_config_path = "/home/roy/Documents/repo/NMRData/input/SH_configs/select_compounds_SH_configs_reduce.json"
surrogate_config_path = "/home/roy/Documents/repo/NMRData/input/surrogate_configs/select_compounds_SH_configs.json"
fit_config_path = "/home/roy/Documents/repo/NMRData/input/fit_configs/align_700MHz_type1_select_compounds.json"

# get mapping from molecule names to their spin system info json files.
H_params_path = "/home/roy/Documents/repo/NMRData/input/coupling_info"
dict_compound_to_filename = JSON.parsefile("/home/roy/Documents/repo/NMRData/input/compound_mapping/select_compounds.json")

# specify the NMR experiment folder
#experiment_full_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/misc/bmse000297_ethanol/"
experiment_full_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/NRC/misc/glucose/Sep-25-2018"
#experiment_full_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/BMRB/similar_settings/BMRB-700-20mM/L-Serine"
#experiment_full_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/BMRB/similar_settings/BMRB-500-0.5mM/L-Serine"
#experiment_full_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/NRC/NRC_4_amino_acid_mixture_Jan_2022/1"
#experiment_full_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/BMRB/similar_settings/BMRB-700-20mM/L-Isoleucine"
#experiment_full_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/BMRB/similar_settings/BMRB-700-20mM/L-Glutamine"
#experiment_full_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/BMRB/similar_settings/BMRB-500-0.5mM/L-Leucine"
#experiment_full_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/misc/bmse000795_2_DSS"
#experiment_full_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/misc/gissmo_DSS"
#experiment_full_path = "/home/roy/MEGAsync/outputs/NMR/experiments/misc/bmse000900_L-Phenylalanine"

#experiment_full_path = "/home/roy/MEGAsync/outputs/NMR/experiments/misc/bmse000900_L-Phenylalanine"
#experiment_full_path = "/home/roy/Documents/repo/NMRData/experiments_1D1H/BMRB/glucose-600-bmse000855_1"

# specify where the calibration results should be saved for this experiment.
#project_name = "ethanol"
project_name = "NRC-glucose-2018"
#project_name = "Serine-BMRB-700-20mM-force-mag-eq"
#project_name = "Serine-BMRB-700-20mM-2000entry"
#project_name = "Serine-BMRB-700-20mM-mod"
#project_name = "Serine-BMRB-500-0.5mM-mod" #project_name = "Serine-BMRB-500-0.5mM"

 #project_name = "Serine-glucose-NRC-Jan2022"
#project_name = "NRC-Jan2022-serine-glucose-dss"

#project_name = "Isoleucine-BMRB-700-20mM"
#project_name = "Glutamine-BMRB-700-20mM"
#project_name = "Leucine-BMRB-500-0.5mM" # skipped.
#project_name = "bmse000795_2_DSS" # skipped.
#project_name = "gissmo_DSS"
#project_name = "Phenylalanine-GISSMO-600-bmse000900"

#project_name = "Phenylalanine-GISSMO-600-bmse000900-poster"
#project_name = "glucose-600-bmse000855_1-poster"

#molecule_names = ["Ethanol";]
molecule_names = ["alpha-D-Glucose"; "beta-D-Glucose";]
#molecule_names = ["D-Glucose - 2000 study";] #molecule_names = ["D-(+)-Glucose";]
#molecule_names = ["D-(+)-Glucose";]
#molecule_names = ["L-Serine";]
#molecule_names = ["L-Serine - 2000 study";]
#molecule_names = ["L-Serine - mod";]
#molecule_names = ["L-Serine - 2000 study";] # molecule_names = ["L-Serine";]

#molecule_names = ["L-Serine - 2000 study"; "D-Glucose - 2000 study";]
#molecule_names = ["L-Serine"; "D-(+)-Glucose"]
#molecule_names = ["L-Serine"; "D-(+)-Glucose"; "L-Isoleucine"]
#molecule_names = ["L-Serine"; "D-(+)-Glucose"; "DSS"] # for debugging code with singlets.
#molecule_names = ["L-Serine"; "D-(+)-Glucose"; "L-Glutamine"; "L-Isoleucine"; "L-Leucine"; "DSS"]

#molecule_names = ["L-Serine"; "D-(+)-Glucose"; "L-Leucine"; ]
#molecule_names = ["L-Serine"; "D-(+)-Glucose"; "DSS"; ] # really need plots that highlights regions and compounds.

#molecule_names = ["L-Isoleucine";]
#molecule_names = ["L-Glutamine";]
#molecule_names = ["L-Leucine";]
#molecule_names = ["DSS";]
#molecule_names = ["DSS";]
#molecule_names = ["L-Phenylalanine";]

#molecule_names = ["L-Phenylalanine";]
#molecule_names = ["D-(+)-Glucose";]

#molecule_names = ["L-Phenylalanine, 500 MHz"; "L-Phenylalanine"]

project_base_folder = "/home/roy/MEGAsync/outputs/NMR/align"
#project_base_folder = "/home/roy/MEGAsync/outputs/NMR/calibrate" # old.
project_folder = joinpath(project_base_folder, project_name)
isdir(project_folder) || mkpath(project_folder)

##### user input for loading data.

## for loading the experiment.
solvent_ppm_guess = 4.7
solvent_window_ppm = 0.1

# start the first entry at the frequency corresponding to offset_ppm.
# This is a method to deal with wrap-around frequency.
offset_ppm = 0.3

## for surrogate construction.

#w = [1.0;] # relative concentration.

Δcs_max_scalar_default = 0.2 # In units of ppm. interpolation border that is added to the lowest and highest resonance frequency component of the mixture being simulated.

dummy_SSFID = NMRSpectraSimulator.SpinSysParamsType1(0.0) # level 2 model.

unique_cs_atol = 1e-6 # for assessing chemical equivalence.

## for determining cost function positions
region_min_dist = 0.1 # in ppm.

#####

##### load experiment, normalize data.
s_t, S, hz2ppmfunc, ppm2hzfunc, ν_0ppm, fs, SW, α_0ppm, β_0ppm, λ_0ppm, Ω_0ppm,
    results_0ppm, dic, α_solvent, β_solvent, λ_solvent, Ω_solvent,
    results_solvent = NMRDataSetup.loadspectrum(experiment_full_path;
    solvent_ppm = solvent_ppm_guess,
    solvent_window_ppm = solvent_window_ppm)
#

# start the first entry at the frequency corresponding to offset_Hz, which we set to 0.3 ppm.

offset_Hz = ν_0ppm - (ppm2hzfunc(offset_ppm)-ppm2hzfunc(0.0))
DFT_s = fft(s_t)
U_DFT, U_y, U_inds = NMRDataSetup.getwraparoundDFTfreqs(length(s_t), fs, offset_Hz)

S_U = DFT_s[U_inds]
P_y = hz2ppmfunc.(U_y)

val, ind = findmin( abs.(U_y .- ν_0ppm) )
Z = abs(S_U[ind])

y = S_U ./ Z
#####

####### compute surrogate.
λ0 = λ_0ppm

# get a surrogate where K_{n,i} is encouraged to be no larger than `early_exit_part_size`.
println("Timing: mag equivalence")
@time MEs = NMRSpectraSimulator.getmageqinfomixture(molecule_names,
    H_params_path,
    dict_compound_to_filename;
    unique_cs_atol = unique_cs_atol)

println("Timing: setupmixtureproxies()")
@time mixture_params = NMRSpectraSimulator.setupmixtureSH(molecule_names,
    H_params_path, dict_compound_to_filename, fs, SW,
    ν_0ppm;
    MEs = MEs,
    config_path = SH_config_path)
As = mixture_params

ΩS_ppm = NMRSpectraSimulator.getPsnospininfo(As, hz2ppmfunc)
ΩS_ppm_sorted = sort(NMRSpectraSimulator.combinevectors(ΩS_ppm))

u_offset = 0.5 # in ppm.
u_min = ppm2hzfunc(ΩS_ppm_sorted[1] - u_offset)
u_max = ppm2hzfunc(ΩS_ppm_sorted[end] + u_offset)

println("Timing: fitproxies!()")
dummy_SSFID = NMRSpectraSimulator.SpinSysParamsType1(0.0)
Bs = NMRSpectraSimulator.fitproxies(As, dummy_SSFID, λ0;
    names = molecule_names,
    config_path = surrogate_config_path,
    Δcs_max_scalar_default = Δcs_max_scalar_default,
    u_min = u_min,
    u_max = u_max)

####### end mixture proxy.

##### prepare fit positions.
### prepare positions.
hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)
P_y = hz2ppmfunc.(U_y)

ΩS0 = NMRSpectraSimulator.getΩS(As)
ΩS0_ppm = NMRSpectraSimulator.getPs(ΩS0, hz2ppmfunc)

Δsys_cs, y_cost_all, U_cost_all, P_cost_all, exp_info, cost_inds, cost_inds_set,
    λ_lbs, λ_ubs, κs_β_orderings,
    κs_β_DOFs = NMRCalibrate.prepareoptim(fit_config_path, molecule_names,
    hz2ppmfunc, U_y, y, As; region_min_dist = region_min_dist)
#
a_setp, b_setp, minxs,
    rets = NMRCalibrate.setupitpab(0.1, 10, 0.7; optim_algorithm = :LN_BOBYQA)

# new normalization.
Z = maximum( maximum(abs.(y[cost_inds_set[r]])) for r = 1:length(cost_inds_set) )
y = y ./ Z

# include("inner_kappa.jl")
#@assert 1==3

# for poster.
#include("pkg_align_single_region.jl")
#include("phenylalanine.jl")

#include("pkg_calibrate.jl")
include("pkg_quantify.jl")

using Plots; plotly()

# save data plot.
file_name = "data.html"
plots_save_path = joinpath(project_folder, file_name)

canvas_size = (1000, 400)
plot_obj = Plots.plot( P_y,
    real.(y),
    title = "data spectrum",
    label = "data",
    seriestype = :line,
    ticks = :native,
    xlims = (P_y[1],P_y[end]),
    hover = P_y,
    linewidth = 4,
    xlabel = "ppm",
    ylabel = "real part of spectrum",
    xflip = true,
    size = canvas_size)
Plots.savefig(plot_obj, plots_save_path)

# julia> minxs
# 3-element Vector{Vector{Float64}}:
#  [-0.3182551682670491, 0.049183131629439616, 0.024203062792075736]
#  [0.9948724914072531, 0.06932430981123272, 0.0839330953494708]
#  [0.060496655154809946, 0.10515911754875247, 0.11996321361619477]
#
# julia> w_BLS
# ERROR: UndefVarError: w_BLS not defined
#
# julia> ws
# 3-element Vector{Vector{Float64}}:
#  [75.8591754032393, 0.21482769137077787]
#  [47.263485918404186, 0.33989540682369757]
#  [0.26376149754144756, 0.34040957238909]

ext: jl


text:
function ngram_char(sequence::String ,n::Int64)
  sequence = split(replace(sequence, " ", "_"), "")
  ngrams = String[]
  for i in 1:length(sequence)-n+1
    push!(ngrams, join(sequence[i:i+n-1]))
  end
  ngrams
end

function ngram_array(sequence ,n::Int64)
  ngrams = String[]
  for i in 1:length(sequence)-n+1
    push!(ngrams, join(sequence[i:i+n-1], " "))
  end
  ngrams
end


text:
"""

Julia client to get price data from flickelectric.co.nz

"""

module Flick

using HTTP
using JSON
using Dates: Day, Minute, Second, DateTime
using ExpiringCaches

const SERVER = "https://api.flick.energy"
const MOBILE_ENDPOINT = SERVER * "/customer/mobile_provider"
const PRICE_URL = MOBILE_ENDPOINT * "/price"
const AUTH_URL = SERVER * "/identity/oauth/token"

const DEFAULT_CLIENT_ID = "le37iwi3qctbduh39fvnpevt1m2uuvz"
const DEFAULT_CLIENT_SECRET = "ignwy9ztnst3azswww66y9vd9zt6qnt"

struct Price
    value::Float64
    unit::String
    issued_at::DateTime
    expires_at::DateTime
end

function get_config(filename)
    JSON.parse(open(filename))
end

"""
    get_auth_token()

Using a user's login credentials, and an application id/secret pair from
Flick, retrieve a authentication token for subsequent API calls.

FYI: Tokens appear to expire in 2 months/60 days.
"""
function get_auth_token(config::Dict)
    get_auth_token(
        get(config, "client_id", DEFAULT_CLIENT_ID),
        get(config, "client_secret", DEFAULT_CLIENT_SECRET),
        config["username"],
        config["password"]
    )
end
get_auth_token(username::String, password::String) = get_auth_token(DEFAULT_CLIENT_ID, DEFAULT_CLIENT_SECRET, username, password)

ExpiringCaches.@cacheable Day(30) function get_auth_token(client_id::String, client_secret::String, username::String, password::String)::String
    headers = ["Content-Type" => "application/x-www-form-urlencoded"]
    payload = Dict(
        "grant_type" => "password",
        "client_id" => client_id,
        "client_secret" => client_secret,
        "username" => username,
        "password" => password
    )
    @debug "⚡ Making authentication call to flick"
    r = HTTP.request("POST", AUTH_URL, headers, HTTP.escapeuri(payload))
    @debug "⚡ Authentication response with status code $(r.status)"

    auth_token = JSON.parse(String(r.body))
    @debug "⚡ Authentication response\n$(auth_token)"

    # We only need the "id_token" string.
    auth_token["id_token"]
end

"""
    get_price_detail()

Returns a Dict from the raw JSON returned by the Flick mobile API.

Contains price broken down into components e.g. the spot price + fees.
"""
function get_price_detail(auth_token)
    flick_api_call(auth_token, PRICE_URL)
end

"""
    flick_api_call(auth_token, url)

Add the Authorization header to a call and parses the assumed JSON response.
"""
function flick_api_call(auth_token, url)::Dict
    @debug "⚡ Making call to flick api @ URL " * url
    headers = [
        "Authorization" => "Bearer $(auth_token)"
    ]
    r = HTTP.request("GET", url, headers)
    @info "⚡ Response received from flick API with status $(r.status)"

    JSON.parse(String(r.body))
end

"""
    get_current_price()

Calls are cached in memory for a couple of minutes to avoid hitting the flick api too often.
"""

ExpiringCaches.@cacheable Minute(2) function get_current_price(auth_token::String)::Price
    flick_price_detail = get_price_detail(auth_token)
    @debug "Price detail:\n$(flick_price_detail)"
    @show flick_price_detail["needle"]
    price_value = parse(Float64, flick_price_detail["needle"]["price"])
    valid_from = DateTime(flick_price_detail["needle"]["start_at"], "y-m-dTHH:MM:SS.sZ")
    valid_until = DateTime(flick_price_detail["needle"]["end_at"], "y-m-dTHH:MM:SS.sZ")

    Price(price_value, flick_price_detail["needle"]["unit_code"], valid_from, valid_to)
end

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

__precompile__()
module LocalModule

using Distributed
export myapp

function myapp()
    p = addprocs(1)
    @everywhere p println("Taking over the world...")
    rmprocs(p)
    nothing
end

end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    GeoDataFrame(data, coordnames)

A dataframe object `data` with additional metadata for tracking
the columns `coordnames` that represent spatial coordinates.

## Examples

If the data was already loaded in a normal DataFrame `data`,
and there exists columns named `x`, `y` and `z`, wrap the
data and specify the column names:

```julia
julia> GeoDataFrame(data, [:x,:y,:z])
```

Alternatively, load the data directly into a `GeoDataFrame` object
by using the method [`readgeotable`](@ref).

### Notes

This type is a lightweight wrapper over Julia's DataFrame types.
No additional storage is required other than a vector of symbols
with the columns names representing spatial coordinates.

"""
struct GeoDataFrame{T<:Real,N,DF<:AbstractDataFrame} <: AbstractSpatialData{T,N}
  data::DF
  coordnames::Vector{Symbol}

  function GeoDataFrame{T,N,DF}(data, coordnames) where {T<:Real,N,DF<:AbstractDataFrame}
    new(data, coordnames)
  end
end

function GeoDataFrame(data, coordnames)
  @assert coordnames ⊆ names(data) "coordnames must contain valid column names"
  Ts = Missings.T.(eltypes(data[coordnames]))
  T  = promote_type(Ts...)
  N  = length(coordnames)
  DF = typeof(data)

  GeoDataFrame{T,N,DF}(data, coordnames)
end

function coordinates(geodata::GeoDataFrame)
  rawdata = geodata.data
  cnames = geodata.coordnames
  ctypes = Missings.T.(eltypes(rawdata[cnames]))

  Dict(var => T for (var,T) in zip(cnames,ctypes))
end

function variables(geodata::GeoDataFrame)
  rawdata = geodata.data
  cnames = geodata.coordnames
  vnames = [var for var in names(rawdata) if var ∉ cnames]
  vtypes = Missings.T.(eltypes(rawdata[vnames]))

  Dict(var => T for (var,T) in zip(vnames,vtypes))
end

npoints(geodata::GeoDataFrame) = nrow(geodata.data)

function coordinates!(buff::AbstractVector, geodata::GeoDataFrame, ind::Int)
  rawdata = geodata.data
  cnames = geodata.coordnames

  for (i, cname) in enumerate(cnames)
    @inbounds buff[i] = rawdata[ind,cname]
  end
end

value(geodata::GeoDataFrame, ind::Int, var::Symbol) = geodata.data[ind,var]

Base.view(geodata::GeoDataFrame, inds::AbstractVector{Int}) =
  GeoDataFrame(view(geodata.data, inds, :), geodata.coordnames)

# ------------
# IO methods
# ------------
function Base.show(io::IO, geodata::GeoDataFrame)
  dims = join(size(geodata.data), "×")
  cnames = join(geodata.coordnames, ", ", " and ")
  print(io, "$dims GeoDataFrame ($cnames)")
end

function Base.show(io::IO, ::MIME"text/plain", geodata::GeoDataFrame)
  println(io, geodata)
  show(io, geodata.data, allcols=true, rowlabel=:Row, summary=false)
end

function Base.show(io::IO, ::MIME"text/html", geodata::GeoDataFrame)
  println(io, geodata)
  show(io, MIME"text/html"(), geodata.data)
end


text:
using LoopVectorization, ProgressMeter
const LOOPVECBENCHDIR = joinpath(pkgdir(LoopVectorization), "benchmark")
include(joinpath(LOOPVECBENCHDIR, "benchmarkflops.jl"))


text:
using JLD
using Plots
using LaTeXStrings
#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(1.75)

lightblue = colorant"rgb(63%,74%,78%)"
orange = colorant"rgb(85%,55%,13%)"
silver = colorant"rgb(69%,67%,66%)"
rust = colorant"rgb(72%,26%,6%)"

#Extracting stored varibles
ld = load("./out/jld/BEM_Plot.jld")
for key in keys(ld)
  symb = Symbol(key)
  data = ld[key]
  @eval $(symb) = $(data)
end
accplot = plot(xlabel=L"$q$",
               ylabel=L"$\log_{10}(\mathrm{RMSE})$",
               # legend=:bottomleft,
               ylim=(-2.5,0.5))
plot!(accplot,
      (3 : 8),
      log10.(sqrt.(accuracy_exact)),
      label=L"\mathrm{dense}",
      markershape=:circle,
      markersize=15,
      color=silver,
      linestyle=:solid,
      linewidth=5)

plot!(accplot,
      3 : 8,
      log10.(sqrt.(accuracy_predict[:, 1])),
      label=L"\rho = 1.0",
      markershape=:dtriangle,
      markersize=15,
      color=rust,
      linestyle=:dash,
      linewidth=5)

plot!(accplot,
      3 : 8,
      log10.(sqrt.(accuracy_predict[:, 2])),
      label=L"\rho = 2.0",
      markershape=:utriangle,
      markersize=15,
      color=orange,
      linestyle=:dash,
      linewidth=5)

plot!(accplot,
      3 : 8,
      log10.(sqrt.(accuracy_predict[:, 3])),
      label=L"\rho = 3.0",
      markershape=:diamond,
      markersize=15,
      color=:lightblue,
      linestyle=:dot,
      linewidth=5)

savefig(accplot, "./out/plots/BEM_accplot.pdf")

nnzplot = plot(xlabel=L"$q$",
               ylabel=L"$\log_{10}\left(\mathrm{Nonzeros}\right)$",
               legend=:topleft,
               #ylim=(-2.5,0.5)
               )
plot!(nnzplot,
      (3 : 8),
      log10.((nnz_exact)),
      label=L"\mathrm{dense}",
      markershape=:circle,
      markersize=15,
      color=silver,
      linestyle=:solid,
      linewidth=5)

plot!(nnzplot,
      3 : 8,
      log10.((nnz_predict[:, 1])),
      label=L"\rho = 1.0",
      markershape=:dtriangle,
      markersize=15,
      color=rust,
      linestyle=:dash,
      linewidth=5)

plot!(nnzplot,
      3 : 8,
      log10.((nnz_predict[:, 2])),
      label=L"\rho = 2.0",
      markershape=:utriangle,
      markersize=15,
      color=orange,
      linestyle=:dash,
      linewidth=5)

plot!(nnzplot,
      3 : 8,
      log10.((nnz_predict[:, 3])),
      label=L"\rho = 3.0",
      markershape=:diamond,
      markersize=15,
      color=:lightblue,
      linestyle=:dot,
      linewidth=5)

savefig(nnzplot, "./out/plots/BEM_nnzplot.pdf")

timeplot = plot(xlabel=L"$q$",
                ylabel=L"$\log_{10}\left(\mathrm{time}\ \mathrm{in}\  \mathrm{seconds}\right)$",
                # legend=:bottomleft,
                #ylim=(-2.5,0.5)
                )
plot!(timeplot,
      (3 : 8),
      log10.((timing_exact)),
      label=L"\mathrm{dense}",
      markershape=:circle,
      markersize=15,
      color=silver,
      linestyle=:solid,
      linewidth=5)

plot!(timeplot,
      3 : 8,
      log10.(sqrt.(timing_predict[:, 1])),
      label=L"\rho = 1.0",
      markershape=:dtriangle,
      markersize=15,
      color=rust,
      linestyle=:dash,
      linewidth=5)

plot!(timeplot,
      3 : 8,
      log10.(sqrt.(timing_predict[:, 2])),
      label=L"\rho = 2.0",
      markershape=:utriangle,
      markersize=15,
      color=orange,
      linestyle=:dash,
      linewidth=5)

plot!(timeplot,
      3 : 8,
      log10.(sqrt.(timing_predict[:, 3])),
      label=L"\rho = 3.0",
      markershape=:diamond,
      markersize=15,
      color=:lightblue,
      linestyle=:dot,
      linewidth=5)

savefig(timeplot, "./out/plots/BEM_timeplot.pdf")


text:
module ExtendedVariantModeExtendedVariantFadeTrigger

using ..Ahorn, Maple

@mapdef Trigger "ExtendedVariantMode/ExtendedVariantFadeTrigger" ExtendedVariantFadeTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	variantChange::String="Gravity", valueA::Integer=10, valueB::Integer=10, positionMode::String="LeftToRight", revertOnDeath::Bool=true)

const placements = Ahorn.PlacementDict(
	"Extended Variant Fade Trigger (Extended Variant Mode)" => Ahorn.EntityPlacement(
		ExtendedVariantFadeTrigger,
		"rectangle"
	)
)

Ahorn.editingOptions(trigger::ExtendedVariantFadeTrigger) = Dict{String, Any}(
	"variantChange" => Ahorn.ExtendedVariantDictionary.Variants,
	"positionMode" => Maple.trigger_position_modes
)

end


text:
using GMT
using SeisIO: read_data, u2d, t_win, d2u
using SeisPDF
using DelimitedFiles
using Dates

const one_hour_length = 3600
const one_hour_step = 1800
const fifteen_minute_length = 900
const fifteen_minute_step = 225
const smooth_width_factor = 1.25
const μs = 1e-6 # some kind of nonsense, the time in SeisIO object is μs

function range!(freq, sampling_rate)
    n = size(freq, 1);
    delta = 1.0 / sampling_rate
    total_duration = delta * n

    for i in 1:n
        freq[i] = i / total_duration
    end
end

function slide_ind(startslide::AbstractFloat,endslide::AbstractFloat,fs::AbstractFloat, ideal_starttime::AbstractFloat)
    starttime = ideal_starttime
    startind = convert(Int,round((startslide - starttime) * fs)) + 1
    endind = convert(Int,round((endslide - starttime) * fs)) + 1
    return startind,endind
end

"""
    nearest_start_end(S::DateTime,E::DateTime, cc_len::Int, cc_step::Int)
Return best possible start, end times for given starttime `S` and endtime `E`.
"""
function nearest_start_end(S::DateTime, E::DateTime, fs::Float64, cc_len::Real, cc_step::Real)
    ideal_start = DateTime(Date(S)) # midnight of same day
    starts = Array(ideal_start:Second(cc_step):E)
    ends = starts .+ Second(cc_len) .- Millisecond(convert(Int,1. / fs * 1e3))
    return d2u(starts[findfirst(x -> x >= S, starts)]), d2u(ends[findlast(x -> x <= E,ends)])
end

input_trace_file = ARGS[1]
response_file = ARGS[2]

S = read_data("mseed", input_trace_file, memmap=true)
println(S.t)
data = S.x[1]
fs = S.fs[1]
data_length = Int(86400 * fs) # Not a good practice
response = read_resp_from_sacpz(response_file, fs, length(data))

#demean!(data)
#detrend!(data)

# 1-hour long segment
#one_hour_starttime = S.t[1][1, 2] * 1.0
one_hour_starttime = S.t[1][1, 2] * μs
one_hour_endtime = one_hour_starttime + data_length / fs - 1 / fs
startslide, endslide = nearest_start_end(u2d(one_hour_starttime), u2d(one_hour_endtime), fs, one_hour_length, one_hour_step)
#println(t_win(S.t[1], S.fs[1]))
#println(one_hour_endtime - one_hour_starttime)
#print(endslide - startslide)
println("$(u2d(one_hour_starttime)), $(u2d(one_hour_endtime))")
println("$(u2d(startslide)), $(u2d(endslide))")
startidx, endidx = slide_ind(startslide, endslide, fs, one_hour_starttime)
println("$startidx, $endidx")
slides_of_one_hour, starts_of_one_hour = slide(@view(data[startidx:endidx]), one_hour_length, one_hour_step, fs, startslide, endslide)
#println(size(starts_of_one_hour, 1))

# PSD mean of all 1-hour segments
_, _, center_periods = get_freqs_and_periods(fs, fifteen_minute_length, smooth_width_factor)
psd_results_mean = Array{Float64, 2}(undef, size(slides_of_one_hour, 2), size(center_periods, 1))
#psd_results_mean = Array{Float64, 2}(undef, 1, size(center_periods, 1))

for i in 1:size(slides_of_one_hour, 2)
#for i in 1:1
    #println("One hour summation $i")
    # 15-minute long segment
    fifteen_minute_starttime = starts_of_one_hour[i]
    fifteen_minute_endtime = fifteen_minute_starttime + length(slides_of_one_hour[:, i]) / fs - 1 / fs
    fifteen_minute_startslide, fifteen_minute_endslide = nearest_start_end(u2d(fifteen_minute_starttime), u2d(fifteen_minute_endtime), fs, fifteen_minute_length, fifteen_minute_step)
    fifteen_minute_startidx, fifteen_minute_endidx = slide_ind(fifteen_minute_startslide, fifteen_minute_endslide, fs, fifteen_minute_starttime)
    slides_of_fifteen_minute, starts_of_fifteen_minute = slide(@view(slides_of_one_hour[fifteen_minute_startidx:fifteen_minute_endidx, i]),
                                                               fifteen_minute_length, fifteen_minute_step, fs, fifteen_minute_starttime, fifteen_minute_endtime)
    #println(size(slides_of_fifteen_minute))
    #println(size(starts_of_fifteen_minute))

    psd_15min_fake = Array{Float64, 2}(undef, Int(fifteen_minute_length * fs), size(slides_of_fifteen_minute, 2))
    #psd_15min_fake = Array{Float64, 2}(undef, Int(fifteen_minute_length * fs), 1)
    #println(size(psd_15min_fake))

    for j in 1:size(slides_of_fifteen_minute, 2)
    #for j in 1:1
        #println("15 minutes summation $j")
        #println(u2d(starts_of_fifteen_minute[j]))
        # Deep copy
        trace = deepcopy(slides_of_fifteen_minute[:, j])
        #for i in 1:10
        #    print("$(trace[i]) ")
        #end
        #println()

        demean!(trace)
        detrend!(trace)

        # Taper the signal
        cosine_taper!(trace, length(trace), 0.05)

        # FFT
        fft_result = compute_fft(trace)
        # Remove instrument response
        remove_response!(fft_result, response)

        # Band-pass filter for preventing overamplification
        freqs = Array{Float32}(undef, length(trace))
        f1, f2, f3, f4 = 0.005, 0.05, 48.0, 50.0
        range!(freqs, fs)
        taper = sac_cosine_taper(freqs, f1, f2, f3, f4, fs)
        for idx in 1:length(trace)
            fft_result[idx] *= taper[idx]
        end

        # Calculate PSD
        psd = calculate_psd(fft_result, fs)

        psd_15min_fake[:, j] = deepcopy(psd)
    end
    #println("===")

    psd_results, _ = summarize_psd(transpose(psd_15min_fake), fs, smooth_width_factor)
    psd_results_mean[i, :] = psd_results[1, :]
end

# Get PDF of this 1-hour slide
#psd_results_mean = reshape(psd_results[1, :], 1, :)
pdf_mean_1_hour = summarize_pdf(psd_results_mean)
#println(pdf_mean_1_hour)
"""open("seispdf_pdf_out.txt", "w") do io
    writedlm(io, pdf_mean_1_hour)
end
open("seispdf_center_periods.txt", "w") do io
    writedlm(io, center_periods)
end
"""
pdf_mean_1_hour = reverse(pdf_mean_1_hour, dims=1)
#imshow(pdf_mean_1_hour, x=periods, y=powers; proj=:logx)
imshow(pdf_mean_1_hour)

# Plot PDF of this 1-hour slide
#period_max = log10(maximum(center_periods))
#period_min = log10(minimum(center_periods))
#center_periods_interval_in_logscale = log10(center_periods[2]) - log10(center_periods[1])
#
## Create netCDF grid
#pdf_mean_1_hour = reshape(pdf_mean_1_hour, :)
#pdf_mean_1_hour_grid = xyz2grd(pdf_mean_1_hour, R="$period_min/$period_max/-200/-50", I="$center_periods_interval_in_logscale/1", Z="TLA", V=true)
#
#
#pdf_mean_1_hour_min = minimum(pdf_mean_1_hour)
#pdf_mean_1_hour_max = maximum(pdf_mean_1_hour)
#g_cpt = makecpt(color=:rainbow, T="$pdf_mean_1_hour_min/$pdf_mean_1_hour_max")
#periods = collect(pdf_mean_1_hour_min:pdf_mean_1_hour_max:center_periods_interval_in_logscale)
#powers = collect(-200:-50:1)
#
##grdview(pdf_mean_1_hour_grid, J="X6i/5i", frame=(xlabel="log10(Period)", ylabel="Power [10log10(m**2/sec**4/Hz)] [dB]", axes=:WSne), color=g_cpt, S=100, Q="s", N=0, V=true, Y="4.0", show=true)
##grdview(pdf_mean_1_hour_grid, J="X6i/5i", frame=(xlabel="log10(Period)", ylabel="Power [10log10(m**2/sec**4/Hz)] [dB]", axes=:WSne, ), color=g_cpt , V=true, show=true, colorbar=true)
#grdview(pdf_mean_1_hour_grid, J="X6i/5i",
#        xaxis=(annot=:auto, ticks=:auto, label="log10(Period)",),
#        yaxis=(annot=:auto, ticks=:auto, label="Power [10log10(m**2/sec**4/Hz)] [dB]",),
#        color=g_cpt, V=true, S=100, Q="s", Y=4.0,
#        colorbar=true, show=true)
#imshow(pdf_mean_1_hour, x=periods, y=powers; proj=:logx)
#imshow(pdf_mean_1_hour, x=periods, y=powers)
#imshow(pdf_mean_1_hour, proj=:log)
#colorbar!(g_cpt, B="0.02", D="6.15i/2.5i/5.0i/0.25i", V=true, show=true)


text:
# --- Imports

import Distributions
using Distributions: MvNormal, Uniform, ContinuousUnivariateDistribution, cdf, pdf
using HypothesisTests
using LinearAlgebra
using Plots
using Random

# --- Generate sample of vectors drawn from a uniform distribution on a unit circle

n = 2
num_samples_2d = 50000

dist = Uniform(-π, π)
num_vectors = num_samples_2d
theta = rand(dist, num_vectors)
vectors = transpose(hcat(cos.(theta), sin.(theta)));

# --- Plot histogram of angles

thetas = map(i -> atan(vectors[:, i][2], vectors[:, i][1]), 1:num_vectors)

num_hist_bins = 25
hist_bins = range(-π, π; length=num_hist_bins)
hist = histogram(thetas; bins=hist_bins, normalize=true)
plt = plot(hist)

pdf_2d = 0.5 / π
println("pdf(θ): $(pdf_2d)")

# Display plot
display(plt)

# --- Perform test for uniformity on the circle

# Analytical formula for distribution over angle from x-axis
struct ThetaDistribution2D <: ContinuousUnivariateDistribution end
Distributions.pdf(dist::ThetaDistribution2D, x::Real) = 0.5 / π
Distributions.cdf(dist::ThetaDistribution2D, x::Real) = 0.5 * (x + π) / π

theta_dist = ThetaDistribution2D()

# Perform Anderson-Darling Test
test_results = OneSampleADTest(thetas, theta_dist)
println(test_results)

test_results = ExactOneSampleKSTest(thetas, theta_dist)
println(test_results)

test_results = ApproximateOneSampleKSTest(thetas, theta_dist)
println(test_results)

# --- Generate sample of vectors drawn from a uniform distribution on a unit circle

n = 2
num_samples_2d = 50000

# Generate vectors
dist = MvNormal(zeros(n), ones(n))
num_vectors = num_samples_2d
vectors = rand(dist, num_vectors)
for i = 1:num_vectors
    vectors[:, i] /= norm(vectors[:, i])
end

# --- Plot histogram of angles

thetas = map(i -> atan(norm(vectors[:, i][2]), vectors[:, i][1]), 1:num_vectors)

num_hist_bins = 25
hist_bins = range(0, π; length=num_hist_bins)
hist = histogram(thetas; bins=hist_bins, normalize=true)
plt = plot(hist)

pdf_2d = 1 / π
println("pdf(θ): $(pdf_2d)")

# Display plot
display(plt)

# --- Perform test for uniformity on the circle

# Analytical formula for distribution over angle from x-axis
struct ThetaDistribution2D <: ContinuousUnivariateDistribution end
Distributions.pdf(dist::ThetaDistribution2D, x::Real) = 1 / π
Distributions.cdf(dist::ThetaDistribution2D, x::Real) = x / π

theta_dist = ThetaDistribution2D()

# Perform Anderson-Darling Test
test_results = OneSampleADTest(thetas, theta_dist)
println(test_results)

test_results = ExactOneSampleKSTest(thetas, theta_dist)
println(test_results)

test_results = ApproximateOneSampleKSTest(thetas, theta_dist)
println(test_results)

# --- Generate sample of vectors drawn from a uniform distribution on a unit sphere

n = 3
num_samples_3d = 50000

# Generate vectors
dist = MvNormal(zeros(n), ones(n))
num_vectors = num_samples_3d
vectors = rand(dist, num_vectors)
for i = 1:num_vectors
    vectors[:, i] /= norm(vectors[:, i])
end

# --- Plot histogram of angles

thetas = map(i -> atan(norm(vectors[:, i][2:end]), vectors[:, i][1]), 1:num_vectors)

num_hist_bins = 50
hist_bins = range(0, π; length=num_hist_bins)
hist = histogram(thetas; bins=hist_bins, normalize=true)
plt = plot(hist)

pdf_3d(x) = 0.5 * sin(x)
plt = plot!(pdf_3d, 0, π)

# Display plot
display(plt)

# --- Perform test for uniformity on the sphere

# Analytical formula for distribution over angle from x-axis
struct ThetaDistribution3D <: ContinuousUnivariateDistribution end
Distributions.cdf(dist::ThetaDistribution3D, x::Real) = 0.5 * (1 - cos(x))

theta_dist = ThetaDistribution3D()

# Perform Anderson-Darling Test
test_results = OneSampleADTest(thetas, theta_dist)
println(test_results)

test_results = ExactOneSampleKSTest(thetas, theta_dist)
println(test_results)

test_results = ApproximateOneSampleKSTest(thetas, theta_dist)
println(test_results)

# --- Generate sample of vectors drawn from a uniform distribution on a unit sphere

n = 20
num_samples = 50000

# Generate vectors
dist = MvNormal(zeros(n), ones(n))
num_vectors = num_samples
vectors = rand(dist, num_vectors)
for i = 1:num_vectors
    vectors[:, i] /= norm(vectors[:, i])
end

# --- Plot histogram of angles

thetas = map(i -> atan(norm(vectors[:, i][2:end]), vectors[:, i][1]), 1:num_vectors)

num_hist_bins = 50
hist_bins = range(0, π; length=num_hist_bins)
hist = histogram(thetas; bins=hist_bins, normalize=true)
plt = plot(hist)

pdf_nd(x) = (2^(n-2) / π / binomial(n-2, (n-2) ÷ 2)) * sin(x)^(n-2)
plt = plot!(pdf_nd, 0, π)

# Display plot
display(plt)

# --- Perform test for uniformity on the hypersphere

# Analytical formula for distribution over angle from x-axis
struct ThetaDistributionND <: ContinuousUnivariateDistribution end

function Distributions.cdf(dist::ThetaDistributionND, x::Real)
    k = n-2
    sgn = (-1)^(k÷2)
    value = 0
    for j = 0:(k÷2 - 1)
        value += (-1)^j * binomial(k, j) * sin((k - 2*j) * x) / (k - 2*j)
    end
    value *= 2
    value += sgn * binomial(k, k÷2) * x  # contribution from i = k÷2 term
    value *= sgn / π / binomial(k, k÷2)

    return value
end

theta_dist = ThetaDistributionND()

# Perform Anderson-Darling Test
test_results = OneSampleADTest(thetas, theta_dist)
println(test_results)

test_results = ExactOneSampleKSTest(thetas, theta_dist)
println(test_results)

test_results = ApproximateOneSampleKSTest(thetas, theta_dist)
println(test_results)

# --- Generate sample of vectors drawn from a uniform distribution on a unit sphere

n = 45
num_samples = 50000

# Generate vectors
dist = MvNormal(zeros(n), ones(n))
num_vectors = num_samples
vectors = rand(dist, num_vectors)
for i = 1:num_vectors
    vectors[:, i] /= norm(vectors[:, i])
end

# --- Plot histogram of angles

thetas = map(i -> atan(norm(vectors[:, i][2:end]), vectors[:, i][1]), 1:num_vectors)

num_hist_bins = 50
hist_bins = range(0, π; length=num_hist_bins)
hist = histogram(thetas; bins=hist_bins, normalize=true)
plt = plot(hist)

pdf_nd(x) = (2^(n-2) / π / binomial(n-2, (n-2) ÷ 2)) * sin(x)^(n-2)
plt = plot!(pdf_nd, 0, π)

# Display plot
display(plt)

# --- Perform test for uniformity on the hypersphere

# Analytical formula for distribution over angle from x-axis
struct ThetaDistributionND <: ContinuousUnivariateDistribution end

function Distributions.cdf(dist::ThetaDistributionND, x::Real)
    k = n-2
    value = 0
    norm_ = 0
    for j = 0:(k-1)÷2
        coef = (-1)^j * binomial(k, j) / (k - 2*j)
        value += coef * (1 - cos((k - 2*j) * x))
        norm_ += coef
    end
    value /= 2 * norm_

    return value
end

theta_dist = ThetaDistributionND()

# Perform Anderson-Darling Test
test_results = OneSampleADTest(thetas, theta_dist)
println(test_results)

test_results = ExactOneSampleKSTest(thetas, theta_dist)
println(test_results)

test_results = ApproximateOneSampleKSTest(thetas, theta_dist)
println(test_results)

max_stars_repo_name: ktchu/RESEARCH-math-exploration

max_issues_repo_name: ktchu/RESEARCH-math-exploration

max_forks_repo_name: ktchu/RESEARCH-math-exploration


text:
export LPHRep

MOI.Utilities.@model(_MOIModel,
                     (), (MOI.EqualTo, MOI.LessThan,), (), (),
                     (), (MOI.ScalarAffineFunction,), (), ())
# We need the `SingleVariable` constraints to be bridged so we should say that
# they are not supported. We notably exclude `Integer` as we just ignore
# integrality constraints. Binary constraint should be bridged to integrality
# once https://github.com/jump-dev/MathOptInterface.jl/issues/704 is done.
function MOI.supports_constraint(
    ::_MOIModel{T}, ::Type{MOI.SingleVariable},
    ::Type{<:Union{MOI.EqualTo{T}, MOI.GreaterThan{T}, MOI.LessThan{T},
                   MOI.Interval{T}, MOI.ZeroOne}}) where T
    return false
end

mutable struct LPHRep{T} <: HRepresentation{T}
    model::_MOIModel{T}
    hyperplane_indices::Union{Nothing,
                              Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},
                                                         MOI.EqualTo{T}}}}
    halfspace_indices::Union{Nothing,
                             Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},
                                                        MOI.LessThan{T}}}}
end
function LPHRep(model::_MOIModel{T}) where T
    return LPHRep(model, nothing, nothing)
end
function LPHRep(model::MOI.ModelLike, T::Type = Float64)
    _model = _MOIModel{T}()
    bridged = MOI.Bridges.LazyBridgeOptimizer(_model)
    # Only enable constraint bridges that don't create variables and don't add
    # any variable bridge so that there is an identity mapping betwenen
    # variables of `model` and polyhedra dimensions.
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.GreaterToLessBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.LessToGreaterBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.NonnegToNonposBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.NonposToNonnegBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.ScalarizeBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.VectorizeBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.ScalarFunctionizeBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.VectorFunctionizeBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.SplitIntervalBridge{T})
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.NormInfinityBridge{T})
    # This one creates variables so the user need to consider the polyhedra as the
    # feasible set of the extended formulation.
    MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.NormOneBridge{T})
    MOI.copy_to(bridged, model)
    return LPHRep(_model)
end
# returns `Int64` so need to convert for 32-bit system
FullDim(rep::LPHRep) = convert(Int, MOI.get(rep.model, MOI.NumberOfVariables()))

hvectortype(::Type{LPHRep{T}}) where {T} = SparseVector{T, Int}
similar_type(::Type{LPHRep{S}}, ::FullDim, ::Type{T}) where {S, T} = LPHRep{T}
fulltype(::Type{LPHRep{T}}) where {T} = LPHRep{T}

LPHRep(h::HRep{T}) where {T} = LPHRep{T}(h)
LPHRep{T}(h::HRep) where {T} = convert(LPHRep{T}, h)

supports_names(::Type{<:LPHRep}) = true
function dimension_names(h::LPHRep)
    if MOI.VariableName() in MOI.get(h.model, MOI.ListOfVariableAttributesSet())
        names = [
            MOI.get(h.model, MOI.VariableName(), MOI.VariableIndex(i))
            for i in 1:fulldim(h)
        ]
    else
        return nothing
    end
end

function _constrain_in_func_set(h::HyperPlane, vars, T)
    return _dot(h.a, vars, T), MOI.EqualTo{T}(h.β)
end
function _constrain_in_func_set(h::HalfSpace, vars, T)
    return _dot(h.a, vars, T), MOI.LessThan{T}(h.β)
end
function _constrain_in(model, hs::HIt, vars, T)
    for h in hs
        MOI.add_constraint(model, _constrain_in_func_set(h, vars, T)...)
    end
end

function LPHRep{T}(d::FullDim,
                   hyperplanes::ElemIt{<:HyperPlane{T}},
                   halfspaces::ElemIt{<:HalfSpace{T}};
                   dimension_names = nothing) where {T}
    model = _MOIModel{T}()
    vars = MOI.add_variables(model, fulldim(d))
    if dimension_names !== nothing
        if length(dimension_names) !== fulldim(d)
            throw(DimensionMismatch("Length of dimension_names ($(length(dimension_names))) does not match the full dimension of the polyhedron ($(fulldim(d)))."))
        end
        for (i, name) in enumerate(dimension_names)
            if !isempty(name)
                MOI.set(model, MOI.VariableName(), MOI.VariableIndex(i), name)
            end
        end
    end
    _constrain_in(model, hyperplanes, vars, T)
    _constrain_in(model, halfspaces, vars, T)
    return LPHRep(model)
end

function Base.copy(lp::LPHRep{T}) where {T}
    model = _MOIModel{T}()
    MOI.copy_to(model, lp.model)
    return LPHRep(model)
end

function constraint_indices(rep::LPHRep{T},
                            ::Union{HyperPlaneIndex, HyperPlaneIndices}) where T
    if rep.hyperplane_indices === nothing
        rep.hyperplane_indices = MOI.get(
            rep.model,
            MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},
                                        MOI.EqualTo{T}}())
    end
    return rep.hyperplane_indices
end
function MOI.add_constraint(rep::LPHRep{T},
                            func::MOI.ScalarAffineFunction{T},
                            set::MOI.EqualTo{T}) where T
    rep.hyperplane_indices = nothing
    return MOI.add_constraint(rep.model, func, set)
end
function MOI.delete(rep::LPHRep{T},
                    ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}) where T
    rep.hyperplane_indices = nothing
    return MOI.delete(rep.model, ci)
end

function constraint_indices(rep::LPHRep{T},
                            ::Union{HalfSpaceIndex, HalfSpaceIndices}) where T
    if rep.halfspace_indices === nothing
        rep.halfspace_indices = MOI.get(
            rep.model,
            MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},
                                        MOI.LessThan{T}}())
    end
    return rep.halfspace_indices
end
function MOI.add_constraint(rep::LPHRep{T},
                            func::MOI.ScalarAffineFunction{T},
                            set::MOI.LessThan{T}) where T
    rep.halfspace_indices = nothing
    return MOI.add_constraint(rep.model, func, set)
end
function MOI.delete(rep::LPHRep{T},
                    ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}) where T
    rep.halfspace_indices = nothing
    return MOI.delete(rep.model, ci)
end

function Base.length(idxs::HIndices{T, LPHRep{T}}) where {T}
    return length(constraint_indices(idxs.rep, idxs))
end
function Base.isempty(idxs::HIndices{T, LPHRep{T}}) where {T}
    return isempty(constraint_indices(idxs.rep, idxs))
end
function startindex(idxs::HIndices{T, LPHRep{T}}) where {T}
    if isempty(idxs)
        return nothing
    else
        return eltype(idxs)(1)
    end
end
function Base.get(rep::LPHRep{T}, idx::HIndex{T}) where {T}
    ci = constraint_indices(rep, idx)[idx.value]
    func = MOI.get(rep.model, MOI.ConstraintFunction(), ci)::MOI.ScalarAffineFunction{T}
    # MOI uses `Int64` but `SparseArrays` uses `Int32` by default so `Int64` will create
    # issues with, e.g. preimages with `spzeros(d, n)`, etc...
    indices = Int[t.variable_index.value for t in func.terms]
    values = [t.coefficient for t in func.terms]
    a = sparsevec(indices, values, FullDim(rep))
    set = MOI.get(rep.model, MOI.ConstraintSet(), ci)
    β = MOI.constant(set) - func.constant
    if idx isa HyperPlaneIndex
        return HyperPlane(a, β)
    else
        return HalfSpace(a, β)
    end
end
function nextindex(rep::LPHRep{T}, idx::HIndex{T}) where {T}
    if idx.value >= length(constraint_indices(rep, idx))
        return nothing
    else
        return typeof(idx)(idx.value + 1)
    end
end

function Base.isvalid(lp::LPHRep{T}, idx::HIndex{T}) where {T}
    return 1 ≤ idx.value ≤ length(constraint_indices(rep, idx))
end

dualtype(::Type{LPHRep{T}}, ::Type{AT}) where {T, AT} = dualtype(Intersection{T, AT, Int}, AT)
dualfullspace(h::LPHRep, d::FullDim, ::Type{T}, ::Type{AT}) where {T, AT} = dualfullspace(Intersection{T, AT, Int}, d, T, AT)


text:
@recipe function f(sol::DGSolution; tidx = size(sol.t,1), vars=nothing)
    xguide --> "x"
    yguide --> "u"
    labels = String[]
    for i in 1:size(sol.u[tidx],2)
      push!(labels,"u$i")
    end
    yvector = sol.u[tidx]
    if vars != nothing
      yvector = sol.u[tidx][:,vars]
      labels = labels[vars]
    end
    if typeof(labels) <: String
      label --> labels
    else
      label --> reshape(labels,1,length(labels))
    end
    @series begin
      seriestype  :=  :path
      sol.nodes, yvector
    end
end

max_stars_repo_name: Paulms/DiscontinuousGalerkin

max_issues_repo_name: Paulms/DiscontinuousGalerkin

max_forks_repo_name: Paulms/DiscontinuousGalerkin


text:
# Compare line drawing performance of Matplotlib vs. GR
#
# The results depend strongly on which backend is used in GR and/or Matplotlib.
# Depending on the selection, however, improvements by a factor of 50 are possible.
# The best results can be achieved with the Qt driver (gksqt).
#
# The increase in drawing speed through the GR Matplotlib backend is comparatively
# small, since most of the time is spent in Matplotlib/Python.

ENV["GKSwstype"] = "gksqt"
ENV["MPLBACKEND"]="module://gr.matplotlib.backend_gr"

import PyPlot
import GR

function mpl()
    x = collect(0:0.01:2*pi)

    PyPlot.plot(x, sin.(x))
    PyPlot.ion()
    t = time_ns()
    for i = 1:100
        PyPlot.cla()
        PyPlot.plot(x .+ i / 10.0, sin.(x .+ i / 10.0))
        if haskey(ENV, "MPLBACKEND")
           PyPlot.show()
        else
           PyPlot.pause(0.0001)
        end
    end

    fps = round(Int64, 100 / (1e-9 * (time_ns() - t)))
end

function gr()
    x = collect(0:0.01:2*pi)

    GR.plot(x, sin.(x))
    t = time_ns()
    for i = 1:100
        GR.plot(x .+ i / 10.0, sin.(x .+ i / 10.0))
    end

    fps = round(Int64, 100 / (1e-9 * (time_ns() - t)))
end

using Printf

for i in 1:10
    fps_mpl = mpl()
    fps_gr = gr()
    speedup = float(fps_gr) / fps_mpl

    @printf("MPL: %4d fps  GR: %4d fps  => speedup: %6.1f\r",
            fps_mpl, fps_gr, speedup)
end

@printf("\n")


text:
#feynmantheory.jl

"""
Implementation of Feynman's original variational technique applied to the Polaron model.

See Feynman 1955:
http://dx.doi.org/10.1103/PhysRev.97.660
"""

using ..PolaronMakie

# Equation 31: The <|X(t) - X(s)|^{-1}> * exp(-|t-w|) effective action.
A_integrand(v, w, τ) = (w^2 * τ + (v^2 - w^2) / v * (1 - exp(-v * τ)))^(-0.5) * exp(-τ)

A(v, w, α) = π^(-0.5) * α * v * QuadGK.quadgk(τ -> A_integrand(v, w, τ), 0, Inf)[1]

# Equation 33: Lowest Free energy E = -B - A where B = -3/(4v)*(v-w)^2.
E(v, w, α) = (3 / (4 * v)) * (v - w)^2 - A(v, w, α)

"""
feynman_variation(α::Float64; v = 7.0, w = 6.0)

    Calculate v and w variational Polaron parameters that minimise the free energy, for the supplied α Frohlich coupling.
    This version uses the original athermal action (Feynman 1955).
    Returns v, w.
"""
function feynman_variation(α; v = 7.0, w = 6.0)
    initial = [v, w]
    lower = [0.0, 0.0]
    upper = [100.0, 100.0]

    f(x) = E(x[1], x[2], α)
    od = Optim.OnceDifferentiable(f, initial; autodiff = :forward)
    solution = Optim.optimize(od, lower, upper, initial, Fminbox(BFGS()))
    v, w = Optim.minimizer(solution)

    return v, w
end


text:
@testset "output" begin
  using LightGraphs
  result = driver(4,2)
  for g in keys(result)
      output(result[g],nv(g),ne(g))
  end
end


text:
module AOC2021
    using Memoize: @memoize
    using MetaGraphsNext, Graphs, SimpleWeightedGraphs
    export solve, FullInput, TestInput
# Write your package code here.

    abstract type Puzzle end
    abstract type FullInput <: Puzzle end
    abstract type TestInput <: Puzzle end

    read(::Type{Puzzle}, path::String) = begin
        open(path) do f
            readlines(f)
        end
    end
    read(::Type{FullInput}, num::Int, part::Int) = read(Puzzle, "$(@__DIR__)/../data/puzzle/$(num)-$(part).txt")
    read(::Type{FullInput}, num::Int) = read(Puzzle, "$(@__DIR__)/../data/puzzle/$(num).txt")
    read(::Type{TestInput}, num::Int, part::Int) = read(Puzzle, "$(@__DIR__)/../data/puzzle_test/$(num)-$(part).txt")
    read(::Type{TestInput}, num::Int) = read(Puzzle, "$(@__DIR__)/../data/puzzle_test/$(num).txt")

    solve(::Type{T}, number::Int, part::Int; only_data=false, kwargs...) where T <: Puzzle = read(T,number) |> input-> solve(Val(number), input) |> x-> only_data ? (return x) : solve( Val(number), Val(part);input=x, kwargs...)
    solve(::Type{T}, number::Int, part::Int, ppart::Int; only_data=false, kwargs...) where T <: Puzzle = read(T,number, ppart) |> input-> solve(Val(number), input) |> x-> only_data ? (return x) : solve( Val(number), Val(part);input=x, kwargs...)

    solve(::Val{1}, input) = shift -> map(x->parse(Int,x), input) |> nums->
        sum(1 for i in (shift+1):length(nums) if nums[i] > nums[i-shift])

    solve(::Val{1}, ::Val{1}; input) = input(1)
    solve(::Val{1}, ::Val{2}; input) = input(3)

    solve(::Val{2}, input) = begin
        global forward = 0
        global depth = 0
        global aim = 0
        for l in input
            d, v = split(l)
            v = parse(Int,v)
            if d == "forward"
                forward += v
                depth += v*aim
            elseif d == "down"
                aim += v
            else
                aim -= v
            end
        end
        (forward=forward, depth=depth, aim=aim)
    end
    solve(::Val{2}, ::Val{1}; input) =  input.forward*input.aim
    solve(::Val{2}, ::Val{2}; input) =  input.forward*input.depth

    solve(::Val{3}, input) = input .|> collect .|> (x->parse.(Int,x)) |> x-> hcat(x...) |> transpose

    sol_3_final_processing(vals) = vals .|>
        (x-> [c > 0 ? '1' : '0' for c in x]) .|>
        join |>
        x-> parse.(Int, x; base=2) |>
        x->reduce(*,x)

    solve(::Val{3}, ::Val{1}; input) = input |>
        x->sum(x;dims=1) |>
        x-> (x.>(size(input)[1]/2)) |>
        x-> [x, .!x] |>
        sol_3_final_processing

    solve(::Val{3}, ::Val{2}; input) = begin
        matr=input
        rows, cols = size(matr)
        ind_selector(m, range, c, select_ones) = begin
            r = length(range)
            if r == 1
                return m[range,:]
            else
                s = sum(m[range,c])
                range = select_ones(s,range) ? (first(range):(first(range) + s-1)) : (first(range)+s):last(range)
                ind_selector(m, range, c+1, select_ones)
            end
        end

        matr |>
            m->sortslices(m, dims=1, rev=true) |>
            m->[ind_selector(m, 1:rows, 1, (s,r)-> 2*s >= length(r)), ind_selector(m, 1:rows, 1, (s,r)-> 2*s < length(r))] |>
            sol_3_final_processing
    end

    struct Board
        b::Matrix{Int}
        found::Matrix{Bool}
        Board(lines) = begin
            board = lines .|> split .|> (x-> parse.(Int,x)) |> x-> hcat(x...)
            new(board, zeros(Bool, size(board)))
        end
    end
    bingo!(n::Int, b::Board) = begin
        b.found[b.b .== n] .= true
        max(sum(b.found;dims=1)...,sum(b.found;dims=2)...) == 5 ?  sum(b.b[.!b.found])*n : nothing
    end

    bingos!(nums::Vector{Int}, boards::Vector{Board}; get_first=true) = begin
        for n in nums
            boards_result = bingo!.(n,boards)
            nothings = isnothing.(boards_result)
            if get_first && !all(nothings)
                res = boards_result[.!nothings]
                length(res) == 1 || throw(DomainError("Too many?"))
                return res[1]
            elseif !get_first
                if all(.!nothings)
                    length(boards_result) == 1 || throw(DomainError("TooMany?"))
                    return boards_result[1]
                else
                    boards = boards[nothings]
                end
            end
        end
    end
    solve(::Val{4}, input) = (nums=parse.(Int,split(first(input),',')), boards=[Board(input[(3:7).+((i-1)*6)]) for i in 1:div(length(input),6)])
    solve(::Val{4}, ::Val{1}; input) = bingos!(input[1], input[2])
    solve(::Val{4}, ::Val{2}; input) =  bingos!(input[1], input[2]; get_first=false)

    struct Grid
        g::Matrix{Int}
        isdiag::Bool
        Grid(line) = begin
            ps = line |> x-> map(strip, split(x, "->")) |> ps -> split.(ps,",") .|> ps -> parse.(Int,ps)
            x1 = ps[1][1]
            x2 = ps[2][1]
            y1 = ps[1][2]
            y2 = ps[2][2]
            max_len = max(abs(x2-x1),abs(y2-y1))

            s(x) = x == 0 ? 1 : sign(x)
            r(a1,a2) = a1 == a2 ? [a1 for _ in 1:max_len+1] : a1:s(a2-a1):a2
            m = zeros(Int,max(x1,x2)+1,max(y1,y2)+1)
            for p in zip(r(x1,x2), r(y1,y2))
                m[p[1]+1,p[2]+1] = 1
            end
            new(m, x1==x2 || y1 == y2)
        end
    end
    data5(line) = begin
        ps = line |> x-> map(strip, split(x, "->")) |> ps -> split.(ps,",") .|> ps -> parse.(Int,ps)
        x1 = ps[1][1]
        x2 = ps[2][1]
        y1 = ps[1][2]
        y2 = ps[2][2]
        max_len = max(abs(x2-x1),abs(y2-y1))

        s(x) = x == 0 ? 1 : sign(x)
        r(a1,a2) = a1 == a2 ? [a1 for _ in 1:max_len+1] : a1:s(a2-a1):a2

        (x1==x2 || y1 == y2, Dict((p[1]+1,p[2]+1) => 1 for p in zip(r(x1,x2), r(y1,y2))))
    end
    solve(::Val{5}, input) = data5.(input)
    solve(::Val{5}, ::Val{1}; input, use_diag=false) = begin
        grids = use_diag ? [x[2] for x in input] : [x[2] for x in input if x[1]]
        mergewith(+,grids...) |> values |> (x-> x.>1) |> sum
    end

    solve(::Val{5}, ::Val{2}; input) = solve(Val(5), Val(1); input, use_diag=true)

    solve(::Val{6}, lines) = lines |>
        l->(length(l) == 1 ? l[1] : throw(DomainError(l, "Should only be one line"))) |>
        s-> split(s,",") |>
        n->parse.(Int,n)

    """ Got some help. I solved 6.1 but 6.2 was too deeply nested so I built the cache first by hand
    which was a bit harder than I had hoped. Saw a solution using the @memoize function, and what a difference.
    Hides so much boilerplate and "JustWorks"
    """
    @memoize fish(counter, days) = begin
        days == 0 && return 1
        counter == 0 && return fish(8,days-1) + fish(6,days-1)
        return fish(counter-1, days-1)
    end

    solve(::Val{6}, ::Val{1}; input, days=80) = begin
        # input |> counter |> (d->[v*fish(k, days) for (k,v) in d]) |> sum
        # No need for the counter since fishes are memoized
        fish1(k) = fish(k,days)
        input .|> fish1 |> sum
    end
    solve(::Val{6}, ::Val{2}; input, days=256) = solve(Val(6), Val(1); input=input,days=days)

    solve(::Val{7}, lines) = parse.(Int, split(lines[1],","))
    solve(::Val{7},::Val{1}; input, red_func=identity) = begin
        global res = 100000000
        for i in UnitRange(extrema(input)...)
            lens = abs.(i.-input)
            s = sum([red_func(l) for l in lens])
            s < res && (res = s)
        end
        res
    end
    solve(::Val{7},::Val{2}; input) = solve(Val(7),Val(1); input=input, red_func=(x->convert(Int,(x*(x+1)/2))))

    solve(::Val{8}, lines) = lines .|> (l -> split(l,"|") .|> strip .|> split)
    count_uniqs_8(res_nums) = res_nums .|> length .|> (x-> x ∈ [2,4, 3, 7]) |> sum
    solve(::Val{8}, ::Val{1}; input) = input .|> (x->count_uniqs_8(x[2])) |> sum
    parse_input_decode_output(nums, vals) = begin
        nums = Set.(nums)
        num_to_str = Dict{Int,Set{Char}}()
        decode = Dict{Char, Char}()

        lengths = length.(nums)

        for (i,(l,n)) in enumerate(zip(lengths,nums))
            l == 2 && (num_to_str[1] = n)
            l == 3 && (num_to_str[7] = n)
            l == 4 && (num_to_str[4] = n)
            l == 7 && (num_to_str[8] = n)
        end
        eight = num_to_str[8]
        decode[setdiff(num_to_str[7],num_to_str[1]) |> pop!] = 'a'
        eg = setdiff(eight,num_to_str[4])
        for c in eg
            if setdiff(eight,c) in nums
                decode[c] = 'e'
                num_to_str[9] = setdiff(eight, c)
            else
                decode[c] = 'g'
            end
        end

        bd = setdiff(num_to_str[4],num_to_str[1])
        for c in bd
            if setdiff(eight,c) in nums
                decode[c] = 'd'
                num_to_str[0] = setdiff(eight, c)
            else
                decode[c] = 'b'
            end
        end
        nums = setdiff(nums,values(num_to_str))
        num_to_str[6] = [n for n in nums if length(n) == 6][1]
        decode[setdiff(eight, num_to_str[6]) |> pop!] = 'c'
        decode[setdiff('a':'g',keys(decode)) |> pop!] = 'f'

        encode = Dict((v,k) for (k,v) in decode)
        num_to_str[3] = setdiff(setdiff(eight,encode['b']),encode['e'])
        num_to_str[5] = setdiff(num_to_str[6], encode['e'])
        num_to_str[2] = setdiff(setdiff(eight,encode['b']),encode['f'])

        str_to_num = Dict((v,k) for (k,v) in num_to_str)

        Set.(vals) .|> (s-> str_to_num[s]) |> join |> x->parse(Int,x)
    end
    solve(::Val{8}, ::Val{2}; input) = input .|> (x->bool_decode(x...)) |> sum

    bool_decode(nums, vals) = begin
        s1(x) = sum(x,dims=1) |> vec
        model = 'a':'g'
        bm = hcat((nums .|> (n -> model .∈ n))...)
        vm = hcat((vals .|> (n -> model .∈ n))...)
        s = sum(bm, dims=1)

        # Index of number in bm
        found = Dict(
            i=>findfirst(isequal(c),s)[2] for (i,c) in zip([1,4,7,8],[2,4,3,7])
        )

        bf(n) = bm[:,found[n]]
        bi(ns) = bm[:,ns]

        ind_069 = (bm .⊻ bf(8)) |> s1 |> x->findall(isequal(1),x)
        found[6] = bi(ind_069) .| bf(7) |> s1 .|> isequal(7) |> x-> popat!(ind_069, findfirst(x))
        found[0] = bi(ind_069) .| bf(4) |> s1 .|> isequal(7) |> x-> popat!(ind_069, findfirst(x))
        found[9] = ind_069[1]

        ind_235 = filter(x->x ∉ values(found), 1:10)
        found[2] = bi(ind_235) .| bf(4) |> s1 .|> isequal(7) |> x->popat!(ind_235, findfirst(x))
        found[5] = bi(ind_235) .| bf(2) |> s1 .|> isequal(7) |> x->popat!(ind_235, findfirst(x))
        found[3] = ind_235[1]
        new_b = bi([found[k] for k in 0:9])
        decode = Dict(v=>k for (k,v) in found)
        eachcol(vm) .|> (x-> x.==bm) .|> (x->all(x,dims=1)) .|> vec .|> findfirst .|> (x->decode[x]) |> join |> x->parse(Int,x)

    end
    pints(x) = parse.(Int,x)
    solve(::Val{9}, lines) = begin
        input = lines .|> collect .|> pints |> x->hcat(x...)
        r,c = size(input)
        m = zeros(Int, r+2,c+2)
        ma = maximum(input) + 10000
        m[:,1] .= ma
        m[:,end] .= ma
        m[1,:] .= ma
        m[end,:] .= ma

        m[2:end-1,2:end-1] = input[:]
        m
    end
    get_low_points(m) = begin
        r,c = size(m)
        check(i,j) =  (([m[i,j-1], m[i,j+1], m[i,j], m[i-1,j], m[i+1,j]] .- m[i,j]) .> 0) |> sum |> isequal(4)
        res = zeros(Bool,r,c)
        res[2:end-1, 2:end-1] .= [check(i,j) for i in 2:r-1, j in 2:c-1]
        res
    end
    solve(::Val{9}, ::Val{1}; input) = input |> get_low_points |> x->sum(input[x].+1)
    valley_size_9(c,valleys) = begin
        shifts = [CartesianIndex(0,1), CartesianIndex(0,-1), CartesianIndex(1,0), CartesianIndex(-1,0)]
        to_test = Set{CartesianIndex{2}}([c])
        done = Set{CartesianIndex{2}}()
        global count = 0
        while !isempty(to_test)
            t = pop!(to_test)
            push!(done,t)
            count += valleys[t]
            new_points = setdiff([t+s for s in shifts if valleys[t+s] != 0], done)
            !isempty(new_points) && push!(to_test, new_points...)
        end
        count

    end
    solve(::Val{9}, ::Val{2}; input) = begin
        m = get_low_points(input)
        valleys = input .< 9
        vs(x) = valley_size_9(x,valleys)
        carts = m |> findall .|> vs |> sort |> reverse |> x->reduce(*,x[1:3])

    end
    error_map_10(c) = begin
        c == ')' && return 3
        c == ']' && return 57
        c == '}' && return 1197
        c == '>' && return 25137
    end
    abstract type Res end
    abstract type Incomplete <: Res end
    abstract type Corrupt <: Res end
    struct Res10{T<:Res}
        val::BigInt

    end
    close_open() = Dict(
        ')' => '(',
        ']' => '[',
        '}' => '{',
        '>' => '<'
    )

    close_10(open) = begin
        open_close = Dict((v,k) for (k,v) in close_open())
        res = Int[]
        close_map = Dict(
            ')' => 1,
            ']' => 2,
            '}' => 3,
            '>' => 4
        )
        global s = 0
        while !isempty(open)
            o = pop!(open)
            c = open_close[o]

            s = 5*s + close_map[c]
        end
        s
    end
    solve(::Val{10}, lines) = lines
    parse_10(line) = begin
        co = close_open()

        opens = Char[]
        for c in line
            if c in values(co)
                push!(opens, c)
                continue
            else
                v = pop!(opens)
                if v == co[c]
                    continue
                else
                    return Res10{Corrupt}(error_map_10(c))
                end
            end
        end
        return Res10{Incomplete}(close_10(opens))
    end
    solve(::Val{10}, ::Val{1}; input) = parse_10.(input) |> x-> reduce(+,[r.val for r in x if isa(r,Res10{Corrupt})])
    solve(::Val{10}, ::Val{2}; input) = parse_10.(input) |> x-> sort([r.val for r in x if isa(r,Res10{Incomplete})]) |> x-> x[convert(Int,(length(x)+1)/2)]

    solve(::Val{11}, lines) = lines .|> collect .|> pints |> x->hcat(x...)'
    flashing_11!(input) = begin
        r,c = size(input)
        view_input(ci) = CartesianIndex(max(ci[1]-1,1),max(ci[2]-1,1)):CartesianIndex(min(ci[1]+1,r), min(ci[2]+1,c))
        input .+= 1
        flashing = Set{CartesianIndex{2}}(findall(isequal(10), input))
        has_flashed = Set{CartesianIndex{2}}()
        while !isempty(flashing)
            f = pop!(flashing)
            push!(has_flashed, f)

            bumps = setdiff(view_input(f), flashing)

            input[bumps] .+= 1

            new_flashing = findall(isequal(10), input)

            !isempty(new_flashing) && push!(flashing, new_flashing...)
        end
        input[input.>9] .= 0
        length(has_flashed)
    end
    solve(::Val{11}, ::Val{1}; input, flashes=100) = sum(flashing_11!(input) for _ in 1:flashes)
    solve(::Val{11}, ::Val{2}; input, flashes = 1_000_000) = begin
        l = length(input)
        for i in 1:flashes
            (l == flashing_11!(input)) && return i
        end
    end

    generate_cave_system(lines) = begin
        mg = MetaGraph(Graph(), VertexMeta = Int, EdgeMeta = Bool)

        for l in lines
            n1,n2 = split(l,'-') .|> Symbol
            for p ∈ (n1,n2)
                haskey(mg, p) && continue
                mg[p] = (lowercase(string(p)) == string(p) ? 1 : 1000000)
            end
            mg[n1, n2] = true
        end
        mg
    end

    find_end(mg::MetaGraph, l::Symbol, history::Vector{Symbol}) = begin
        history = deepcopy(history)
        push!(history, l)
        n = code_for(mg,l)
        if l == :end
            return Set([history])
        end
        nl = [label_for(mg, nb) for nb in neighbors(mg,n)]
        isempty(nl) && return Set([push!(history, :fail)])

        mg[l] -= 1
        mg[l] == 0 && rem_vertex!(mg, n)

        union([find_end(deepcopy(mg), nb, history) for nb in nl]...)

    end
    get_all_paths(input) = find_end(input, :start, Symbol[]) |> x -> filter(y->y[end] != :fail,x)
    solve(::Val{12}, lines) = generate_cave_system(lines)
    solve(::Val{12}, ::Val{1}; input) = get_all_paths(input) |> length

    solve(::Val{12}, ::Val{2}; input) = begin
        small_caves = [label_for(input, i) for i in 1:nv(input) if input[label_for(input,i)] == 1 && label_for(input,i) ∉[:start,:end]]
        inc_small(ip, c) = begin
            ip = deepcopy(ip)
            ip[c] += 1
            ip
        end
        caves = union([get_all_paths(inc_small(input, c)) for c in small_caves]...)
        length(caves)
    end

    parse_lines_13(lines) = begin
        coords = Vector{Tuple{Int,Int}}()
        folders = Vector{Tuple{Symbol,Int}}()

        parse_folders_line(l) = begin
            h,t = split(l,'=')
            (Symbol(h[end]),parse(Int,t))
        end
        parse_folders(lns) = [parse_folders_line(l) for l in lns]
        for (i,line) in enumerate(lines)
            if isempty(line)
                push!(folders, parse_folders(lines[i+1:end])...)
                break
            end
            push!(coords,parse.(Int,split(line,',')) |> Tuple)
        end

        (coords, folders)
    end
    mirror(cs,xs,::Val{:y}) = Set([(i, j < xs ? j : 2*xs -j) for (i,j) in cs])
    mirror(cs,ys,::Val{:x}) = Set([(i < ys ? i : 2*ys - i, j) for (i,j) in cs])

    solve(::Val{13}, lines) = parse_lines_13(lines)

    solve(::Val{13}, ::Val{1}; input) = input |> x->mirror(x[1],x[2][1][2], Val(x[2][1][1])) |> length
    solve(::Val{13}, ::Val{2}; input) = begin
        coords, foldrs = input
        global couts = Set([(i+1,j+1) for (i,j) in coords])
        foldrs = [(s,l+1) for (s,l) in foldrs]

        for f in foldrs
            couts = mirror(couts,f[2],Val(f[1]))
        end
        couts
    end

    with_vectors(input, iters=10) = begin
        bs, mapping = input
        global bs = bs
        for _ in 1:iters
            global bs_tmp = Vector{Char}(undef, 2*length(bs)-1)
            for (i,e) in enumerate(bs[1:end-1])
                n = bs[i+1]
                bs_tmp[2*i-1] = e
                bs_tmp[2*i+1] = n
                bs_tmp[2*i] = mapping[(e,n)]
            end
            bs = bs_tmp
        end
        res = Dict{Char,Int}()
        for i in bs
            res[i] = get(res, i, 0) + 1
        end
        maximum(values(res)) - minimum(values(res))
    end
    with_loop(input, iters) = begin
        chars, mapping = input
        pairs = Dict((k,0) for k in keys(mapping))

        for i in 1:length(chars)-1
            pairs[(input.chars[i], input.chars[i+1])] += 1
        end
        for i in 1:iters
            ps = [(v, k,(k[1],mapping[k]),(mapping[k],k[2])) for (k,v) in pairs if v > 0]
            for (v, o,n1,n2) in ps
                pairs[o] -= v
                pairs[n1] += v
                pairs[n2] += v
            end
        end
        res = Dict{Char,Int}()
        for ((k,_),v) in pairs
            res[k] = get(res, k, 0) + v
        end
        res[chars[end]] = get(res, chars[end], 0) + 1
        res
    end
    solve(::Val{14}, lines) = begin
        @assert(isempty(lines[2]))
        basestr = lines[1]
        totuple(x,y) = Tuple(x),only(y)
        (chars=collect(basestr), mapping=Dict(totuple(strip.(split(l,"->"))...) for l in lines[3:end]))
    end
    solve(::Val{14}, ::Val{1}; input, iters=10) = begin
        res = with_loop(input, iters)
        l,h = extrema(values(res))
        h-l
    end
    solve(::Val{14},::Val{2}; input, iters=40) = solve(Val(14),Val(1); input, iters)

    solve(::Val{15}, lines) = begin
        risk_score = zeros(Int, length(lines), length(lines[1]))
        for (i,l) in enumerate(lines)
            for (j,c) in enumerate(l)
                risk_score[i,j] = parse(Int,c)
            end
        end
        risk_score
    end

    prep_dijkstra(input) = begin
        r,c = size(input)
        g = SimpleWeightedDiGraph(r*c)

        ind_to_node = Dict{Tuple{Int,Int},Int}()
        node_to_ind = Dict{Int, Tuple{Int,Int}}()
        get_inds(i,j) = Set([
            (i, min(c, j+1)),
            (i, max(1, j-1)),
            (min(r, i+1), j),
            (max(1, i-1), j)
            ])
        for i in 1:r, j in 1:c
            node = (j+c*(i-1))
            ind_to_node[(i,j)] = node
            node_to_ind[node] = (i,j)
        end
        for (n,ind) in node_to_ind
            for ind2 in get_inds(ind...)
                n2 = ind_to_node[ind2]
                add_edge!(g, n2, n, input[ind2...])
            end
        end
        ps = dijkstra_shortest_paths(g,1).dists[end]
        ps - input[1,1] + input[end,end]
    end
    prep_dijkstra_large(input) = begin
        r,c = size(input)
        m = zeros(Int, 5*r, 5*c)
        for i in 1:5, j in 1:5
            m[1+r*(i-1):r*i,1 + c*(j-1):c*j] = ((input .+ (i+j-3)) .% 9) .+1
        end
        prep_dijkstra(m)
    end
    solve(::Val{15}, ::Val{1}; input) = prep_dijkstra(input)
    solve(::Val{15}, ::Val{2}; input) = prep_dijkstra_large(input)
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct SingleTimeSeriesMetadata <: StaticTimeSeriesMetadata
        name::String
        resolution::Dates.Period
        initial_timestamp::Dates.DateTime
        time_series_uuid::UUIDs.UUID
        length::Int
        scaling_factor_multiplier::Union{Nothing, Function}
        internal::InfrastructureSystemsInternal
    end

A TimeSeries Data object in contigous form.

# Arguments
- `name::String`: user-defined name
- `resolution::Dates.Period`
- `initial_timestamp::Dates.DateTime`: time series availability time
- `time_series_uuid::UUIDs.UUID`: reference to time series data
- `length::Int`: length of this time series
- `scaling_factor_multiplier::Union{Nothing, Function}`: Applicable when the time series data are scaling factors. Called on the associated component to convert the values.
- `internal::InfrastructureSystemsInternal`
"""
mutable struct SingleTimeSeriesMetadata <: StaticTimeSeriesMetadata
    "user-defined name"
    name::String
    resolution::Dates.Period
    "time series availability time"
    initial_timestamp::Dates.DateTime
    "reference to time series data"
    time_series_uuid::UUIDs.UUID
    "length of this time series"
    length::Int
    "Applicable when the time series data are scaling factors. Called on the associated component to convert the values."
    scaling_factor_multiplier::Union{Nothing, Function}
    internal::InfrastructureSystemsInternal
end

function SingleTimeSeriesMetadata(name, resolution, initial_timestamp, time_series_uuid, length, scaling_factor_multiplier=nothing, )
    SingleTimeSeriesMetadata(name, resolution, initial_timestamp, time_series_uuid, length, scaling_factor_multiplier, InfrastructureSystemsInternal(), )
end

function SingleTimeSeriesMetadata(; name, resolution, initial_timestamp, time_series_uuid, length, scaling_factor_multiplier=nothing, internal=InfrastructureSystemsInternal(), )
    SingleTimeSeriesMetadata(name, resolution, initial_timestamp, time_series_uuid, length, scaling_factor_multiplier, internal, )
end

"""Get [`SingleTimeSeriesMetadata`](@ref) `name`."""
get_name(value::SingleTimeSeriesMetadata) = value.name
"""Get [`SingleTimeSeriesMetadata`](@ref) `resolution`."""
get_resolution(value::SingleTimeSeriesMetadata) = value.resolution
"""Get [`SingleTimeSeriesMetadata`](@ref) `initial_timestamp`."""
get_initial_timestamp(value::SingleTimeSeriesMetadata) = value.initial_timestamp
"""Get [`SingleTimeSeriesMetadata`](@ref) `time_series_uuid`."""
get_time_series_uuid(value::SingleTimeSeriesMetadata) = value.time_series_uuid
"""Get [`SingleTimeSeriesMetadata`](@ref) `length`."""
get_length(value::SingleTimeSeriesMetadata) = value.length
"""Get [`SingleTimeSeriesMetadata`](@ref) `scaling_factor_multiplier`."""
get_scaling_factor_multiplier(value::SingleTimeSeriesMetadata) = value.scaling_factor_multiplier
"""Get [`SingleTimeSeriesMetadata`](@ref) `internal`."""
get_internal(value::SingleTimeSeriesMetadata) = value.internal

"""Set [`SingleTimeSeriesMetadata`](@ref) `name`."""
set_name!(value::SingleTimeSeriesMetadata, val) = value.name = val
"""Set [`SingleTimeSeriesMetadata`](@ref) `resolution`."""
set_resolution!(value::SingleTimeSeriesMetadata, val) = value.resolution = val
"""Set [`SingleTimeSeriesMetadata`](@ref) `initial_timestamp`."""
set_initial_timestamp!(value::SingleTimeSeriesMetadata, val) = value.initial_timestamp = val
"""Set [`SingleTimeSeriesMetadata`](@ref) `time_series_uuid`."""
set_time_series_uuid!(value::SingleTimeSeriesMetadata, val) = value.time_series_uuid = val
"""Set [`SingleTimeSeriesMetadata`](@ref) `length`."""
set_length!(value::SingleTimeSeriesMetadata, val) = value.length = val
"""Set [`SingleTimeSeriesMetadata`](@ref) `scaling_factor_multiplier`."""
set_scaling_factor_multiplier!(value::SingleTimeSeriesMetadata, val) = value.scaling_factor_multiplier = val
"""Set [`SingleTimeSeriesMetadata`](@ref) `internal`."""
set_internal!(value::SingleTimeSeriesMetadata, val) = value.internal = val

get_horizon(val::SingleTimeSeriesMetadata) = get_length(val)


text:
# Makie.jl recipes

# if Makie is loaded, use Makie by default
PLOTTING_BACKEND[] = :Makie

# encase the recipe in a module to avoid namespace pollution
module MakieRecipe

using ..Makie

# define an internal plotting recipe for shape values
@recipe(ParetoShapePlot, pareto_shape) do scene
    l_theme = Makie.default_theme(scene, Makie.LineSegments)
    return Makie.Attributes(;
        Makie.default_theme(scene, Makie.Scatter)...,
        linewidth=l_theme.linewidth,
        linecolor=(:grey, 0.7),
        showlines=false,
    )
end

end # module

using .MakieRecipe: ParetoShapePlot

# main plotting function
function _paretoshapeplot(backend::Val{:Makie}, args...; attributes...)
    return Makie.plot(ParetoShapePlot, args...; attributes...)
end
function _paretoshapeplot!(backend::Val{:Makie}, args...; attributes...)
    return Makie.plot!(ParetoShapePlot, args...; attributes...)
end

function Makie.plot!(p::ParetoShapePlot)
    attrs = p.attributes
    points = p[1]
    lattr = (:linewidth, :linecolor, :showlines)
    linewidth, linecolor, showlines = getindex.(p, lattr)
    xminmax = lift(points) do pts
        xmin, xmax = extrema(first, pts)
        xspan = xmax - xmin
        return [xmin - 0.01 * xspan, xmax + 0.01 * xspan]
    end
    thresholds = (0, 0.5, 0.7, 1)
    linestyles = (:dot, :dashdot, :dash, :solid)
    for (thresh, linestyle) in zip(thresholds, linestyles)
        Makie.linesegments!(
            p,
            xminmax,
            fill(thresh, 2);
            linestyle=linestyle,
            linewidth=linewidth,
            color=linecolor,
            visible=showlines,
        )
    end
    ks = filter(∉(lattr), keys(attrs))
    return Makie.scatter!(p, p[1]; (k => attrs[k] for k in ks)...)
end

# designate default plot type for PSISResult
Makie.plottype(::PSISResult) = ParetoShapePlot

# convert PSISResult to its shape values when plotted
function Makie.convert_arguments(P::Makie.PlotFunc, result::PSISResult)
    return Makie.convert_arguments(P, as_array(missing_to_nan(pareto_shape(result))))
end
function Makie.convert_arguments(
    ::Type{<:ParetoShapePlot}, pareto_shape::Union{Real,AbstractVector{<:Real}}
)
    return Makie.convert_arguments(Makie.Scatter, as_array(pareto_shape))
end

# set default labels when an Axis is available
function Makie.plot!(
    ax::Makie.Axis, P::Type{<:ParetoShapePlot}, attributes::Makie.Attributes, args...
)
    plt = invoke(
        Makie.plot!,
        Tuple{typeof(ax),Makie.PlotFunc,typeof(attributes),typeof.(args)...},
        ax,
        P,
        attributes,
        args...,
    )
    _set_default_makie_axis_labels!(ax)
    return plt
end
# dispatch on cases like plot(::PSISResult)
function Makie.plot!(
    ax::Makie.Axis, P::Type{Any}, attributes::Makie.Attributes, result::PSISResult
)
    plt = Makie.plot!(ax, ParetoShapePlot, attributes, result)
    _set_default_makie_axis_labels!(ax)
    return plt
end

function _set_default_makie_axis_labels!(ax)
    if isempty(ax.xlabel[])
        ax.xlabel = "Parameter index"
    end
    if isempty(ax.ylabel[])
        ax.ylabel = "Pareto shape"
    end
    return ax
end


text:
using LinearAlgebra
using SparseArrays
using Plots

"This code solves the 2D Laplace's equation with nonzero Dirichlet and Neumann
boundary conditions using a finite difference method."

m = 100 # number of points
x = LinRange(-1,1,m+2)
y = LinRange(-1,1,m+2)
xint = x[2:end-1]
yint = y[2:end-1]
h = x[2]-x[1]

# ordering of unknowns
# u14 u24 u34 u44
# u13 u23 u33 u43
# u12 u22 u32 u42
# u11 u21 u31 u41

id(i,j,m) = i + (j-1)*m # return global index from i,j
function FD_matrix_2D(h,m)
    A = spzeros(m*m,m*m)
    for i = 1:m, j = 1:m # loop thru indices
        A[id(i,j,m),id(i,j,m)] = 4/h^2
        if i > 1 # avoids leftmost line of nodes
            A[id(i,j,m),id(i-1,j,m)] = -1/h^2 # x-derivative
        end
        if i < m # avoids rightmost line of nodes
            A[id(i,j,m),id(i+1,j,m)] = -1/h^2
        end
        if j > 1
            A[id(i,j,m),id(i,j-1,m)] = -1/h^2 # y-derivative
        end
        if j < m
            A[id(i,j,m),id(i,j+1,m)] = -1/h^2
        end
    end
    return A
end

A = FD_matrix_2D(h,m)

f(x,y) = exp(-10*((x-.25)^2 + (y-.75)^2))
# f(x,y) = 1.0
b = vec(f.(xint,yint'))

# solve Au = f(xi,yj)
u = A\vec(b)
contourf(xint,yint,u)


text:
using Base.Test

include("secret-handshake.jl")

@testset "wink for 1" begin
    @test secret_handshake(1) == ["wink"]
end

@testset "double blink for 10" begin
    @test secret_handshake(2) == ["double blink"]
end

@testset "close your eyes for 100" begin
    @test secret_handshake(4) == ["close your eyes"]
end

@testset "jump for 1000" begin
    @test secret_handshake(8) == ["jump"]
end

@testset "combine two actions" begin
    @test secret_handshake(3) == ["wink", "double blink"]
end

@testset "reverse two actions" begin
    @test secret_handshake(19) == ["double blink", "wink"]
end

@testset "reversing one action gives the same action" begin
    @test secret_handshake(24) == ["jump"]
end

@testset "reversing no actions still gives no actions" begin
    @test secret_handshake(16) == []
end

@testset "all possible actions" begin
    @test secret_handshake(15) == ["wink", "double blink", "close your eyes", "jump"]
end

@testset "reverse all possible actions" begin
    @test secret_handshake(31) == ["jump", "close your eyes", "double blink", "wink"]
end

@testset "do nothing for zero" begin
    @test secret_handshake(0) == []
end

@testset "do nothing if lower 5 bits not set" begin
    @test secret_handshake(32) == []
end


text:
@memoize function ξ(Pk,ℓ,n; q="best") # Eq 18 in SVM + Eq 3 in 2Fast.
    N = 16*1024
    kmino = exp(-25)
    kmaxo = exp(25)
    r0o = 1/kmaxo
    rmini = exp(-11)#exp(-14.555)
    rmaxi = exp(11)#exp(11)

    if q == "auto"
        if ℓ == 2 && n == -2
            q = 0
        else
            q = max(.5, 1+n)
        end
    end

    if q == "best"
        q = best_q(Pk,ℓ,n, kmino, kmaxo)
        if ℓ == 0
            q += .4
        end
        #if ℓ == 1
        #    if n == 0
        #        q -= 1
        #    end
        #end
    end

    if q == "show"
        q = show_best_q(Pk,ℓ,n, kmino, kmaxo)
        if ℓ == 0
            println("Xi l0 correction")
            q += .4
        end
        #if ℓ == 1
        #    if n == 0
        #        println("Xi l1n0 correction")
        #        q -= 1
        #    end
        #end
    end

    if q == "henry"
        q = henry_q(Pk,ℓ,n, kmino, kmaxo, r=kmaxo)#r=rmaxi)
    end

    r, xi = xicalc(Pk, ℓ, -n, N=N, kmin=kmino, kmax=kmaxo, r0=r0o; q=q)
    return r, xi .* (r .^ (-n))
end

max_stars_repo_name: JosephTomlinson/GeneralBiasPk

max_issues_repo_name: JosephTomlinson/GeneralBiasPk

max_forks_repo_name: JosephTomlinson/GeneralBiasPk


text:
#using Compat.Test
#using Compat

using ViscousFlow
using Test
##using TestSetExtensions

#@test isempty(detect_ambiguities(ViscousFlow))
include("pointforce.jl")
include("fields.jl")
include("points.jl")
include("layers.jl")
include("saddle.jl")
include("timemarching.jl")
include("systems.jl")

#@testset ExtendedTestSet "All tests" begin
#    @includetests ARGS
#end

#if isempty(ARGS)
#    include("../docs/make.jl")
#end


text:
abstract type PoissonSolver{dType} end

"""
Solves the Poisson equation for periodic boundary conditions

```
function solve!(p::PoissonSolver, x::AbstractVector, w::AbstractVector) end
```

"""
function solve! end

"""
Evaluates the charge density field
"""
function eval_density end

"""
Evaluates the electrostatic potential
"""
function eval_potential end

"""
Evaluates the electric field
"""
function eval_field end

function eval_density(p::PoissonSolver{DT}, x::AbstractVector{DT}) where {DT}
    [eval_density(p,x_) for x_ in x]
end

function eval_potential(p::PoissonSolver{DT}, x::AbstractVector{DT}) where {DT}
    [eval_potential(p,x_) for x_ in x]
end

function eval_field(p::PoissonSolver{DT}, x::AbstractVector{DT}) where {DT}
    [eval_field(p,x_) for x_ in x]
end

function eval_density!(y::AbstractVector{DT}, p::PoissonSolver{DT}, x::AbstractVector{DT}) where {DT}
    for i in eachindex(x,y)
        y[i] = eval_density(p,x[i])
    end
    return y
end

function eval_potential!(y::AbstractVector{DT}, p::PoissonSolver{DT}, x::AbstractVector{DT}) where {DT}
    for i in eachindex(x,y)
        y[i] = eval_potential(p,x[i])
    end
    return y
end

function eval_field!(y::AbstractVector{DT}, p::PoissonSolver{DT}, x::AbstractVector{DT}) where {DT}
    for i in eachindex(x,y)
        y[i] = eval_field(p,x[i])
    end
    return y
end


text:
using Brownian
using Test

include("test_brownian_motion.jl")
include("test_fbm.jl")


text:
# Variables which can be multiple values
NotEnv = Union{Float64, Bool, AbstractBool, Array, Int64, Interval}
struct EnvVar{K,V}
  worlds::Dict{K,NotEnv}

  EnvVar() = new(Dict{K,V}())
  EnvVar(worlds::Dict{K,V}) = new(worlds)
end
EnvVar() = EnvVar{Any,Any}()
singleton{T}(x::T) = (s = Set{T}();push!(s,x);s)
const noconstraints = Set{Symbol}()

typealias IntervalEnvDict Dict{Set{Symbol},Interval}
typealias IntervalEnvVar EnvVar{Set{Symbol}, Interval}

function intervalenvvar{T<:Real}(x::T,y::T)
  EnvVar{Set{Symbol},Interval}(IntervalEnvDict([noconstraints => Interval(x,y)]))
end

unitinterval(::Type{EnvVar}) = intervalenvvar(0.,1.)

function getindex(e::EnvVar, i::Int64)
  ret = EnvVar()
  for world in e.worlds
    val = world[2][i]
    if isa(val, EnvVar)
      for valworld in val.worlds
        ret.worlds[valworld[1]] = valworld[2]
      end
    else
      ret.worlds[world[1]] = world[2][i]
    end
  end
  ret
end

function getindex(e::EnvVar, i::Int64, j::Int64)
  ret = EnvVar()
  for world in e.worlds
    val = world[2][i,j]
    if isa(val, EnvVar)
      for valworld in val.worlds
        ret.worlds[valworld[1]] = valworld[2]
      end
    else
      ret.worlds[world[1]] = world[2][i,j]
    end
  end
  ret
end

# function getindex(e::EnvVar, i::Int64, j::Int64)
#   ret = EnvVar()
#   for world in e.worlds
#     ret.worlds[world[1]] = world[2][i,j]
#   end
#   ret
# end

# Set functions
function subsumes(a::AbstractBool,e::EnvVar)
  doessubsume = true
  for world in e.worlds
    doessubsume = doessubsume & subsumes(a,world[2])
  end
  doessubsume
end

function convert(::Type{Vector{EnvVar}}, b::HyperBox)
  x = [intervalenvvar(b.intervals[1,i],b.intervals[2,i]) for i = 1:num_dims(b)]
end

function isintersect(e::EnvVar, a::AbstractBool)
  doesisintersect = false # only has to isintersect with one
  for world in e.worlds
    doesisintersect = doesisintersect | isintersect(a,world[2])
  end
  doesisintersect
end

function isintersect(a::AbstractBool, e::EnvVar)
  doesisintersect = false # only has to isintersect with one
  for world in e.worlds
    doesisintersect = doesisintersect | isintersect(a,world[2])
  end
  doesisintersect
end

ConcreteValue = Union{Float64, Int64, Bool}

#FIXME, IVE GOT IN and SIZE AS ASSOCIATE BUT ITS NOT
for op = (:+, :-, :*, :>, :>=, :<=, :<, :&, :|, :in, :/, :size)
  @eval begin
    function ($op)(x::EnvVar, y::EnvVar)
      ret = EnvVar() #FIXME: MAKE TYPE STABLE
      for xworld in x.worlds
        xworldid = xworld[1]
        if haskey(y.worlds,xworldid)     # If we share the same world
          ret.worlds[xworldid] = ($op)(xworld[2],y.worlds[xworldid])
        else
          for yworld in y.worlds
            conjworld = union(xworldid,yworld[1])
            ret.worlds[conjworld] = ($op)(xworld[2],yworld[2])
          end
        end
      end
      ret
    end
  end

  @eval begin
    function ($op)(x::EnvVar, y::ConcreteValue)
      ret = EnvVar() #FIXME: MAKE TYPE STABLE
      for xworld in x.worlds
        xworldid = xworld[1]
        ret.worlds[xworldid] = ($op)(xworld[2],y)
      end
      ret
    end
  end

  @eval begin
    function ($op)(y::ConcreteValue, x::EnvVar)
      ret = EnvVar() #FIXME: MAKE TYPE STABLE
      for xworld in x.worlds
        xworldid = xworld[1]
        ret.worlds[xworldid] = ($op)(y,xworld[2])
      end
      ret
    end
  end
end

for op = (:sqr, :sqrt, :inv)
  @eval begin
    function ($op)(x::EnvVar)
      ret = EnvVar() #FIXME: MAKE TYPE STABLE
      for xworld in x.worlds
        xworldid = xworld[1]
        ret.worlds[xworldid] = ($op)(xworld[2])
      end
      ret
    end
  end
end

function update_ret!(a::EnvVar,ret::EnvVar, constraints)
  for aworld in a.worlds
    ret.worlds[union(aworld[1],constraints)] = aworld[2]
  end
end

function update_ret!(a::ConcreteValue, ret::EnvVar, constraints)
  ret.worlds[constraints] = a
end

function add_path_constraints(e::EnvVar, constraints)
  ret = EnvVar()
  for world in e.worlds
    ret.worlds[union(world[1],constraints)] = world[2]
  end
  ret
end

add_path_constraints(e::ConcreteValue, constraints) = e

function update_ret!(a::Array, ret::EnvVar, constraints)
  amap = map(x->add_path_constraints(x,constraints),a)
  ret.worlds[constraints] = amap
end

macro Iff(condition, conseq, alt)
  local idtrue = singleton(gensym())
  local idfalse = singleton(gensym())
  q =
  quote
  c =  $(esc(condition));
  local ret
  if isa(c, EnvVar)
    ret = EnvVar()
    @show world[2]
    for world in c.worlds
      if world[2] === T || world[2] === true
        ret.worlds[world[1]] = $(esc(conseq))
      elseif world[2] === F || world[2] === false
        ret.worlds[world[1]] = $(esc(alt))
      elseif world[2] === TF
        a = $(esc(conseq))
        constraintstrue = union(world[1],$idtrue)
        update_ret!(a,ret, constraintstrue)

        b = $(esc(alt))
        constraintsfalse = union(world[1],$idfalse)
        update_ret!(b,ret, constraintsfalse)

      else
        println("error:", world[2])
        throw(DomainError())
      end
    end
  end
  ret
  end
  return q
end


text:
"""
    function add_cable(

       )
"""
function add_cable(
      network::Dict,
      cable_params::String,
      )

      params = JSON.parsefile(cable_params)

end


text:
depsfile = joinpath(dirname(@__FILE__),"path.jl")

if isfile(depsfile)
    rm(depsfile)
end

function write_depsfile(path)
    open(depsfile, "w") do io
        write(io, """const baron_exec = $(repr(path))\n""")
    end
end

function ci_installation()
    files = if Sys.islinux()
    [
        (ENV["SECRET_BARON_LIN64_JUMP_DEV"], "baron")
        (ENV["SECRET_BARON_LIC_JUMP_DEV"], "baronlice.txt")
    ]
    end
    for (url, file) in files
        local_filename = joinpath(@__DIR__, file)
        download(url, local_filename)
        chmod(local_filename, 0o777)
        if file == "baron"
            write_depsfile(local_filename)
        end
    end
end

if haskey(ENV, "BARON_JL_SKIP_LIB_CHECK")
    # Skip!
elseif get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", "false") == "true"
    write_depsfile("julia_registryci_automerge")
elseif get(ENV, "SECRET_BARON_LIN64_JUMP_DEV", "") != ""
    ci_installation()
elseif haskey(ENV, "BARON_EXEC")
    path = ENV["BARON_EXEC"]
    write_depsfile(path)
else
error(
"""
Unable to locate BARON executable.
Make sure the solver has been separately downloaded from https://minlp.com/baron-downloads
And that you properly set the `BARON_EXEC` environment variable.
"""
)
end


text:
# set paths for simulation -
const _PATH_TO_ROOT = @__DIR__
const _PATH_TO_SRC = joinpath(_PATH_TO_ROOT, "src")
const _PATH_TO_CONFIG = joinpath(_PATH_TO_ROOT, "config")

# load external packages -
using TOML
using Parameters
using UUIDs

# load my codes -
include(joinpath(_PATH_TO_SRC, "VLTypes.jl"))
include(joinpath(_PATH_TO_SRC, "VLExtensions.jl"))
include(joinpath(_PATH_TO_SRC, "VLDefaultExchangeLogic.jl"))
include(joinpath(_PATH_TO_SRC, "VLDefaultAgentLogic.jl"))
include(joinpath(_PATH_TO_SRC, "VLGameEngine.jl"))
include(joinpath(_PATH_TO_SRC, "VLFactory.jl"))