text:
"""
hmm_logic : logic associated with dealing with Hidden Markov Models and
associated files, as well as construction of internal versions of this model.
"""

using BioAlignments, BioSymbols, Unicode

islower(s) = all(c->islowercase(c) | isspace(c), s);
#isupper(s) = all(c->isuppercase(c) | isspace(c), s);

align_problem = GlobalAlignment()
align_blosum62 = BLOSUM62;
#Custom parameters that seem to handle "edges" of alignment well: we want a pretty heavy penalty on gaps.
align_model = AffineGapScoreModel(align_blosum62, match=2, mismatch=-1, gap_open=-20, gap_extend=-5)

global INF_NEG
INF_NEG = -1e12

global bases
bases = ["A", "C", "G", "T"]

function get_hmm_len(hmm_file)
	out_stats = read(pipeline(`hmmstat $hmm_file`), String)
	final_line = split(out_stats, "\n")[end-1]
	the_len = split(final_line)[6]
	return parse(Int64, the_len)
end

function state_probs(hmm::types.wholeHMM, node_num::Int)
	if node_num <= size(hmm.state_probs)[1]
		return hmm.state_probs[node_num, 1:end] #this should be 7 things long...
	else
		#Parameters appropriate for end of sequence.
		return [-0.02312, -3.7785, -1e5, -0.61958, -0.77255, 0.0000, -1e5]
	end
end

function emit_match_probs(hmm::types.wholeHMM, node_num::Int, aas::Array{Char, 1})
	return Dict{Char, Float64}(aas[i] => hmm.match_probs[node_num, i] for i in 1:20) #depends on aa, alphabetical. Think it's fine.
end

function emit_insert_probs(hmm::types.wholeHMM, node_num::Int, aas::Array{Char, 1})
	return Dict(aas[i] => hmm.insert_probs[node_num, i] for i in 1:20)
end

function hmm_get_consensus(hmm, hmm_length::Int)
	#This is a simple "return most likely residue at each position" consensus.
	aas = ['A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y']
	consensus = ""
	for pos in 1:(hmm_length-2)
		consensus *= string(aas[findmax(hmm.match_probs[pos, 1:end])[2]])
	end
	consensus *= "*"
	return consensus
end

function create_HMM(file_name)
	#Extracts parameters from hmm files. Returns them to construct HMM objects.
	in_file = open(file_name)
	in_read = readlines(in_file)
	close(in_file)

	line_count = 0
	block_mod = 0
	last_line = ""
	skip_block = false

	state_probs = zeros(Float64, 0, 0)
	insert_aa_probs = zeros(Float64, 0, 0)
	match_aa_probs = zeros(Float64, 0, 0)
	hmm_length = 0
	alphabet = ""
	alphabet_length = 20 #for now, will be set again downstream
	collecting = false
	cur_row = -1 #-1 because we start collecting at a distinctive "m->m" line which exists in every hmm file but is ahead of actual parameters.

	for line in in_read[1:end-1]

		if !(collecting)
			if line[1:4] == "LENG"
				hmm_length = parse(Int64, split(line)[2])
				state_probs = zeros(Float64, hmm_length, 7)
				insert_aa_probs = zeros(Float64, hmm_length, alphabet_length)
				match_aa_probs = zeros(Float64, hmm_length, alphabet_length)
			end
		end
		if (collecting)
			section = block_mod % 3
			elements = split(line)

			if section == 0
				cur_row += 1
				hmm_pos = elements[1]
				if hmm_pos == "COMPO"
					skip_block = true
				else
					skip_block = false
					probs = Float64[]
					for ele in elements[2:alphabet_length+1]
						if ele == "*"
							push!(probs, 0.0)
						else
							push!(probs, -parse(Float64, ele))
						end
					end
					match_aa_probs[cur_row, 1:end] = probs
				end

			elseif (section == 1) && !(skip_block)
				elements = split(line)
				probs = Float64[]
				for ele in elements[1:alphabet_length]
					if ele == "*"
						push!(probs, 0.0)
					else
						push!(probs, -parse(Float64, ele))
					end
				end
				insert_aa_probs[cur_row, 1:end] = probs
			elseif (section == 2) && !(skip_block)
				probs = Float64[]
				for ele in elements
					if ele == "*"
						push!(probs, 0.0)
					else
						push!(probs, -parse(Float64, ele))
					end
				end
				state_probs[cur_row, 1:end] = probs
			end

		elseif (occursin("m->m", line)) && (occursin("m->i", line)) && (occursin("i->i", line))
			collecting = true
			alphabet = split(last_line)[2:end]
			alphabet_length = length(alphabet)
		end

		block_mod += 1

		last_line = strip(line)
		line_count += 1
	end

	return hmm_length, state_probs, insert_aa_probs, match_aa_probs
end

function load_hmm_seq(prot_name, hmm, hmm_len)
	prot_dict = Dict{AbstractString, AbstractString}(prot_name => hmm_get_consensus(hmm, hmm_len))
	nucl_dict = Dict{AbstractString, AbstractString}(prot_name => optimize_codons(prot_dict[prot_name]))
	return prot_dict, nucl_dict
end

top_codon = Dict{Char, AbstractString}('A'=>"GCG", 'C'=>"TGC", 'D'=>"GAC", 'E'=>"GAA", 'F'=>"TTT",
																							'G'=>"GGT", 'H'=>"CAT", 'I'=>"ATT", 'K'=>"AAA", 'L'=>"CTG",
																							'M'=>"ATG", 'N'=>"AAC", 'P'=>"CCG", 'Q'=>"CAG", 'R'=>"CGT",
																							'S'=>"AGC", 'T'=>"ACC", 'V'=>"GTG", 'W'=>"TGG", 'Y'=>"TAT")

function get_hmm_consensus(hmm_file) #NOTE: extra dependency (hmmemit)
	hmm_con = read(`hmmemit -c $hmm_file`, String)
	return join(split(hmm_con, '\n')[2:end], "")
end

function get_tails(hmmalign) #HMM alignments sometimes struggle near end of sequence. This is a helper function to collect ends of sequence.
	#We correct the problems with the HMM alignment by just using a standard pairwise aligner.
	cur_pos = 1
	seq_to_date = Char[]
	while islower(hmmalign[cur_pos]) || hmmalign[cur_pos] == '-'
		push!(seq_to_date, hmmalign[cur_pos])
		cur_pos += 1
	end

	end_count = 7
	while end_count > 0 && cur_pos < length(hmmalign)
		if isuppercase(hmmalign[cur_pos])
			push!(seq_to_date, hmmalign[cur_pos])
			cur_pos += 1
			end_count -= 1
		else
			end_count = 0
		end
	end

	out_str = join(seq_to_date)
	return out_str
end

function get_consensus(hmm_file)
	get_hmm_consensus(hmm_file)
end

function get_seqs(seq_keys, gene_name)
	for k in seq_keys
		if occursin("gene:$gene_name ", k)
			return k
		end
	end
	return "N/A"
end

function align_grem(gremlin_prot, hmm_file)
	#This just aligns the HMM consensus with the MRF model. We do this to map between both sequence models.
	#Often the HMM has more insertion points than the MRF since HMMs are better at handling insertions. MRFs we just prune nodes.
	the_str = ">grem\n$gremlin_prot\n"
	hmm_consensus = get_hmm_consensus(hmm_file)
	hmm_out = read(pipeline(`printf "$the_str"`, `hmmalign --outformat A2M $hmm_file -`), String)
	out_hmm = join(split(hmm_out, '\n')[2:end], "")
	out_hmm = fix_hmmalign(out_hmm, hmm_consensus)
	return out_hmm
end

function align_seq(seq_prot, hmm_file, hmm_consensus) #this function is very similar to last, slightly different input parameters.
	the_str = ">seq\n$seq_prot\n"
	hmm_out = read(pipeline(`printf "$the_str"`, `hmmalign --outformat A2M $hmm_file -`), String)
	out_hmm = join(split(hmm_out, '\n')[2:end], "")
	out_hmm = fix_hmmalign(out_hmm, hmm_consensus)
	return out_hmm
end

function redo_alignment(grem_text, sample_nucs, is_mark_not_deg)
	#More logic around fixing issues with model alignments.
	#Here we're trying to map out insertions relative to MRF/HMM.
	k = 0
	skip = Dict{Int64, Array{Int64, 1}}()
	deletions = Dict{Int64, Array{Int64, 1}}()

	for ak in 1:length(sample_nucs) #making a skip/deletion vector for each member of sample_nucs.
		skip[ak] = Int64[]
		deletions[ak] = Int64[]
	end

	for sn in sample_nucs
		k += 1
		nat_path = sn.path[1:end]
		if is_mark_not_deg
			first_aa = sn.path[end][1]
		else
			first_aa = sn.path[end][2]
		end

		for_path = reverse(nat_path) #this is the forward path.
		if is_mark_not_deg
			short_path = [x[1] for x in for_path]
		else
			short_path = [x[2] for x in for_path]
		end

		hmm_pos = 0 #keeps track of HMM-node position. Doesn't increase on lower-case letters.
		grem_pos = 0 #keeps track of MRF-node position. What we are mapping to.
		path_pos = 0
		last_path_node = 0
		last_normal = 1
		while_flag = false

		for char in 1:length(grem_text)
			if !(islower(grem_text[char]))
				hmm_pos += 1
			end

			if !(grem_text[char] == '-')
				grem_pos += 1
			end

			if hmm_pos < minimum(short_path) || hmm_pos > maximum(short_path)
				#Then we are in unaligned zone, we just need to look at what grem/hmm are up to.
				if grem_text[char] == '-'
					push!(skip[k], grem_pos)
				elseif islower(grem_text[char])
					push!(deletions[k], last_normal)
				else
					last_normal = grem_pos
				end

			else
				if !(while_flag)
					path_pos += 1
				end
				while_flag = false

				if path_pos > length(short_path)
					println("Warning: when redoing alignment, HMM_pos is still in progress path is finished.")
				else
					cur_path_pos = short_path[path_pos]
				end

				#Here we're reacting to something specific to sampling trajectory, we deal with it. This is independent of alignment between HMM/GREM.
				while path_pos < length(short_path) && cur_path_pos == last_path_node
					while_flag = true
					#that means CAMEOS stayed at the same hmm node. An insertion.
					#if there is an insertion here from MRF, presumably it would consume the HMM "insertion"
					if islower(grem_text[char])
						#We don't consider this skippable because there is an insertion from the MRF offsetting the HMM insertion.
						grem_pos += 1
						char += 1
					else
						#If just business as usual, though, then we skip the hmm insertion.
						push!(skip[k], grem_pos)
					end
					#Move orward...
					last_path_node = short_path[path_pos]
					last_path_pos = path_pos
					path_pos += 1
					cur_path_pos = short_path[path_pos]
				end

				if !(while_flag)
				#So this is again the territory where not related to skips in path but between HMM/MRF.
					if grem_text[char] == '-'
						#Then an aa in HMM isn't seen in MRF. We skip this.
						#There's a chance it's deleted. In practice this almost never happens.
						push!(skip[k], grem_pos)
					elseif islower(grem_text[char])
						#Then there is an insertion from MRF vs the HMM.
						#Unless there is an insertion from HMM (which would be covered above) this means we have a deletion.
						push!(deletions[k], last_normal)
						path_pos -= 1 #this is not related to path, we need to go back one.
					elseif (short_path[path_pos] - last_path_node > 1) && last_path_node != 0 #because then it's just the first iteration and we don't have a last_path_node.
						for s in 1:(short_path[path_pos] - last_path_node)
							push!(deletions[k], last_normal)
						end
					else
						last_normal = grem_pos
					end
					cur_path_pos = short_path[path_pos]
					last_path_node = short_path[path_pos]
					last_path_pos = path_pos
				end
			end
		end
	end

	return skip, deletions
end

function to_aa(seq)
	return map(x -> convert(AminoAcid, x), collect(strip(seq)))
end

function fix_hmmalign(hmmalign, hmm_consensus)
	#So hmmalign can be very gappy and insertiony near the ends of sequences... so we want to do a global pairwise alignment on it.
	#We are now going to use Bio.jl's pairalign() function.
	pure_seq = uppercase(replace(hmmalign, "-"=>""))

	start_seq = get_tails(hmmalign)
	original_up_len = length(start_seq)
	edited_start = uppercase(replace(start_seq, "-"=>""))
	first_aligned = findfirst(z -> isuppercase(z), collect(start_seq)) #inserts before an aligned char are not inserts, they are unaligned.

	up_len = length(start_seq) - count(z -> islower(z), start_seq)
	rev_seq = reverse(get_tails(reverse(hmmalign)))
	original_down_len = length(rev_seq)
	edited_end = uppercase(replace(rev_seq, "-"=>""))
	down_len = length(rev_seq) - count(z -> islower(z), rev_seq)

	head_seq = uppercase(string(hmm_consensus[1:up_len]))
	tail_seq = uppercase(string(hmm_consensus[end - down_len + 1: end]))

	#align_problem is a global function that is just defining we want a global alignment here.

	head_align = alignment(pairalign(align_problem, to_aa(edited_start),
							to_aa(head_seq), align_model))
	tail_align = alignment(pairalign(align_problem, to_aa(edited_end),
							to_aa(tail_seq), align_model))

	head_ali_top = collect(join([x for (x, _) in head_align]))
	head_ali_bot = collect(join([x for (_, x) in head_align]))

	for head_count in 1:length(head_ali_top)
		if head_ali_bot[head_count] == '-'
			head_ali_top[head_count] = lowercase(head_ali_top[head_count])
		end
	end
	head_ali_top = join(head_ali_top)

	tail_ali_top = collect(join([x for (x, _) in tail_align]))
	tail_ali_bot = collect(join([x for (_, x) in tail_align]))

	for tail_count in 1:length(tail_ali_top)
		if tail_ali_bot[tail_count] == '-'
			tail_ali_top[tail_count] = lowercase(tail_ali_top[tail_count])
		end
	end
	tail_ali_top = join(tail_ali_top)

	return head_ali_top * hmmalign[original_up_len + 1 : (end - original_down_len)] * tail_ali_top
end

function hmm_name(gene_name)
	hmm_path_file = open("hmkt_gene_table.txt")
	hmm_read = readlines(hmm_path_file)
	close(hmm_path_file)
	hmm_read = filter(x -> !isempty(x), hmm_read)
	for line in hmm_read
		if '\t' in line
			spl = split(strip(line), '\t')
			if spl[1] == gene_name
				return strip(spl[4])
			end
		end
	end
	return ""
end

#Basically we also find it useful to trace the protein sequence to the hmm.
#We want wild-type sequence for non-overlapping regions, and we want to align it to hmm because hmm is sort of our "base positioning"
#of the models we're constructing.
function gen_hmm_trace(gene_name, gene_hmm)
	cds_seqs = bio_seq.load_fasta("cds.fasta")
	prot_seqs = bio_seq.load_fasta("proteins.fasta")

	local cds_seq, prot_seq
	cds_seq = cds_seqs[gene_name]
	prot_seq = prot_seqs[gene_name]

	hmm_consensus = get_consensus(gene_hmm)
	hmm_align = align_seq(prot_seq, gene_hmm, hmm_consensus)

	i_count = 1
	hmm_counter = 1
	to_delete = Int64[]
	to_add = Dict{Int64, Char}()

	for i in hmm_align
		if i == '-'
			to_add[i_count] = hmm_consensus[hmm_counter]
		end
		if islower(i)
			push!(to_delete, i_count)
		else
			hmm_counter += 1
		end
		i_count += 1
	end

	final_aas = Char[]
	final_codons = AbstractString[]

	aa_count = 1
	cds_count = 1
	for aa in hmm_align #hmm_align is a protein sequence, aligned.
		if aa_count in to_delete
			#Then we don't want to include this codon in subsequent step.
			cds_count += 1
		else
			if !(aa_count in keys(to_add)) #this is the normal case.
				push!(final_aas, aa)
				push!(final_codons, cds_seq[((cds_count - 1) * 3) + 1: cds_count * 3])
				cds_count += 1
			elseif aa_count in keys(to_add)
				push!(final_aas, to_add[aa_count])
				push!(final_codons, top_codon[to_add[aa_count]])
				#we don't iterate on cds_count because we haven't added natural codons.
			end
		end
		aa_count += 1 #aa_count has more to do with alignment.
	end

	final_cds = join(final_codons)
	final_prot_seq = join(final_aas)
	return final_cds, final_prot_seq
end

function optimize_codons(sequence::AbstractString)
	#most popular codons taken from first table at http://openwetware.org/wiki/Escherichia_coli/Codon_usage [super simple way to find codons, not claiming it's very smart]
	most_popular = Dict{Char, AbstractString}('G' => "GGC", 'E' => "GAA", 'D' => "GAT", 'V' => "GTG", 'A' => "GCG",
                  'R' => "CGC", 'K' => "AAA", 'N' => "AAC", 'M' => "ATG", 'I' => "ATT",
                  'T' => "ACC", 'W' => "TGG", 'C' => "TGC", '*' => "TAA", 'F' => "TTT",
                  'S' => "AGC", 'Q' => "CAG", 'H' => "CAT", 'L' => "CTG", 'P' => "CCG", 'Y'=>"TAT")
	return join([most_popular[aa] for aa in sequence])
end

function top_n(el_tuple, val_keys)
	a_max = findmax(el_tuple)
	max_key::Int = val_keys[a_max[2]]
	max_value = a_max[1]
	return (max_key, max_value)
end

### All of these give_log() functions look at the HMM probabilities of regions *around* the overlap.
#The idea is that this allows us to calculate full-length HMM probabilities, not just probabilities for overlap region.
function follow_prot_give_log_2d(hmm, prot, limit)
	sum_logs = Float32[0]
	aas = ['A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y']
	for position in 1:(limit-2)
		push!(sum_logs, sum_logs[position] + state_probs(hmm, position)[1] + emit_match_probs(hmm, position, aas)[prot[position]])
	end
	if limit - 1 > 0
		push!(sum_logs, sum_logs[limit-1] + state_probs(hmm, limit -1)[1])
	end
	return sum_logs * ones(Float32, 1, 64) #I think this makes a big matrix of this vector.
end

function follow_prot_give_log(hmm, prot, limit)
	sum_logs = Float64[0]
	aas = ['A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y']
	for position in 1:(limit-2)
		push!(sum_logs, sum_logs[position] + state_probs(hmm, position)[1] + emit_match_probs(hmm, position, aas)[prot[position]])
	end
	if limit - 1 > 0
		push!(sum_logs, sum_logs[limit-1] + state_probs(hmm, limit -1)[1])
	end
	return sum_logs
end

function tail_prot_give_log_2d(hmm, prot, limit)
	sum_logs = Float64[0]
	aas = ['A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y']
	for position in 1:(limit-2)
		push!(sum_logs, sum_logs[position] + state_probs(hmm, position)[1] + emit_match_probs(hmm, position, aas)[uppercase(prot[position])])
	end
	push!(sum_logs, sum_logs[limit-1] + state_probs(hmm, limit - 1)[1])
	total = sum_logs[end]
	return [total - sum_logs[y] for y in 1:length(sum_logs)]
end

function tail_prot_give_log(hmm, prot, limit)
	sum_logs = Float64[0]
	aas = ['A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y']
	for position in 1:(limit-2)
		push!(sum_logs, sum_logs[position] + state_probs(hmm, position)[1] + emit_match_probs(hmm, position, aas)[prot[position]])
	end
	push!(sum_logs, sum_logs[limit-1] + state_probs(hmm, limit - 1)[1])
	total = sum_logs[end]
	return [total - sum_logs[y] for y in 1:length(sum_logs)]
end


text:
using PSCOPF

using Test
using JuMP
using Dates
using DataStructures

@testset verbose=true "test_energy_market_slacks" begin

    #=
    If available power for a limitable exceeds consumption,
      excess power is capped.
    The capping variable is global (delocalised in space but has
      a per timestep and per scenario value)
    The power capped is distributed in the schedule (not visible here cause only one limitable unit).

    TS: [11h]
    S: [S1]
                      bus 1
                        |
    (limitable) wind_1_1|
    Pmin=0, Pmax=100    |
    Csta=0, Cprop=1     |     load_1
      S1: 20            | S1: 15
      S2: 25            | S2: 25
                        |
    =#
    @testset "energy_market_capping_limitable_power" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 20.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 25.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 15.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 25.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict{String, SortedDict{Dates.DateTime, PSCOPF.DecisionFirmness} }(),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE))
                    )
        # initial generators state : No need because all pmin=0 => ON by default
        generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)
        market = PSCOPF.EnergyMarket()
        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        # Limitable produces to the available level
        @test 15. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        @test 25. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S2")

        # Capping is not localised in market model
        # Limitable was capped when prod > load (ie. S1):
        @test 5. ≈ value(result.limitable_model.p_global_capping[TS[1], "S1"])
        # Limitable was not capped when prod <= load (ie. S2):
        @test value(result.limitable_model.p_global_capping[TS[1], "S2"]) < 1e-09

        # Capping is localised in the market schedule
        @test 5. ≈ PSCOPF.get_capping(context.market_schedule, "wind_1_1", TS[1], "S1")
        @test PSCOPF.get_capping(context.market_schedule, "wind_1_1", TS[1], "S2") < 1e-09

        # No penalty but we do pay for capped power
        #(otherwise we need a per unit variable for capping to localise and reduce the unpaid costs)
        @test value(result.objective_model.penalty) < 1e-09
        @test value(result.objective_model.start_cost) < 1e-09
        @test value(result.objective_model.prop_cost) < 1e-09 # no cost for limitables
    end

    #=
    It is possible to cut consumption.
    This variable should can be used in two cases :
    1- due to EOD constraints : we don't have enough production capacity (illustrated in S1)
    2- due to Pmin constraints : we don't have enough demand to start a unit (illustrated in S2)

    TS: [11h]
    S: [S1]
                        bus 1
                         |
    (pilotable) prod_1_1 |    load_1
    Pmin=20, Pmax=100    |  S1: 150
    Csta=100k, Cprop=1   |  S2: 15
                         |  S3: 25
    =#
    @testset "energy_market_cuts_consumption_due_to_pmin" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                20., 100.,
                                                100000., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 150.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 15.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S3", 25.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),),
                    SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE))
                    )
        # initial generators state :
        generators_init_state = SortedDict("prod_1_1" => PSCOPF.OFF)

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing, "DBGUG")
        market = PSCOPF.EnergyMarket()
        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # indicates using slacks for feasibility
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK
        # S1 : prod_capacity < load => cannot satisfy demand
        @test 100. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1")
        @test 50. ≈ value(result.lol_model.p_global_loss_of_load[TS[1], "S1"])
        # S2 : load < pmin => cannot start the unit for such load
        @test PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S2") < 1e-09
        @test value(result.limitable_model.p_global_capping[TS[1], "S2"]) < 1e-09
        @test 15. ≈ value(result.lol_model.p_global_loss_of_load[TS[1], "S2"])
        # S3 : works fine
        @test 25. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S3")
        @test value(result.lol_model.p_global_loss_of_load[TS[1], "S3"]) < 1e-09
        # penalize cutting consumption
        @test 1e7 ≈ market.configs.loss_of_load_penalty
        @test (50. * 1e7 + 15. * 1e7 + 0. ) ≈ value(result.objective_model.penalty)
        @test (1e5 + 0. + 1e5) ≈ value(result.objective_model.start_cost)
        @test (100. + 0. + 25. ) ≈ value(result.objective_model.prop_cost)

        #cut conso is localized in the ranscripted schedule:
        @test 50. ≈ PSCOPF.get_loss_of_load(context.market_schedule, "bus_1", TS[1], "S1")
        @test 15. ≈ PSCOPF.get_loss_of_load(context.market_schedule, "bus_1", TS[1], "S2")
        @test PSCOPF.get_loss_of_load(context.market_schedule, "bus_1", TS[1], "S3") < 1e-09
    end

    #=
    Cutting consumption is done through penalization.
    Careful for the chosen penalty cost with respect to production costs
    Here, the unit starting cost is higher than the penalty

    TS: [11h]
    S: [S1]
                        bus 1
                         |
    (pilotable) prod_1_1 |    load_1
    Pmin=20, Pmax=100    |  S1: 25
    Csta=100k, Cprop=1   |

    Desired solution :
        prod_1_1 at 25MW
    Retrieved solution with low penalization (1e3):
        prod_1_1 at 0
        loss_of_load : 25MW
    =#
    @testset "energy_market_careful_for_cuts_consumption_penalty" begin
        penalty = 1e3

        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                20., 100.,
                                                100000., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 25.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),),
                    SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE))
                    )
        # initial generators state :
        generators_init_state = SortedDict("prod_1_1" => PSCOPF.OFF)

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)
        market = PSCOPF.EnergyMarket()
        market.configs.loss_of_load_penalty = penalty
        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Desired Solution
        @test !( PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL )
        @test !( 25. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1") )
        @test !( value(result.lol_model.p_global_loss_of_load[TS[1], "S1"]) < 1e-09 )
        @test !( value(result.objective_model.penalty) < 1e-09 )

        # retrieved solution
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK
        @test PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1") < 1e-09
        @test 25. ≈ value(result.lol_model.p_global_loss_of_load[TS[1], "S1"])
        @test (25. * 1e3) ≈ value(result.objective_model.penalty)

    end

    #=
    In S1,
    demand = 15, wind_1_1 provides 10 => residual demand of 5MW but pmin(prod_1_1)=20
    => start prod_1_1 at 20MW exceeds the demand
    => use the 10MW from wind_1_1 and cut the remaining 5MW (ie LoL)
    In S2,
    demand = 25, wind_1_1 provides 10 => residual demand of 15MW but pmin(prod_1_1)=20
    => start prod_1_1 at 20MW
    => use 5MW from wind_1_1 and cap 5MW

    TS: [11h]
    S: [S1]
                        bus 1
                         |
    (limitable) wind_1_1 |    load_1
    Pmin=0, Pmax=100     |  S1: 15
    Csta=0, Cprop=0.     |  S2: 25
         S1 : 10         |
         S1 : 10         |
                         |
    (pilotable) prod_1_1 |
    Pmin=20, Pmax=100    |
    Csta=100k, Cprop=100 |
    =#
    @testset "energy_market_capping_limitables_due_to_pilotable_pmin" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                20., 100.,
                                                100000., 100.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 10.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 10.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 15.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 25.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                               "prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE))
                    )
        # initial generators state :
        generators_init_state = SortedDict("prod_1_1" => PSCOPF.OFF)

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)
        market = PSCOPF.EnergyMarket()
        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # indicates using slacks for feasibility because of scenario S1
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        # In S1 : Load=15, wind provides 10 => still missing 5 but pmin=20
        # => reduce consumption by 5
        @test 10. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1") < 1e-09
        @test value(result.limitable_model.p_global_capping[TS[1], "S1"]) < 1e-09
        @test 5. ≈ value(result.lol_model.p_global_loss_of_load[TS[1], "S1"])

        # In S2 : Load=25, wind provides 10 => still missing 15 but pmin=20
        # pilotable produces 20 => 5 extra prod (20+10 - 25) => reduce wind by 5.
        @test 5. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S2")
        @test 20. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S2")
        @test 5. ≈ value(result.limitable_model.p_global_capping[TS[1], "S2"])
        @test value(result.lol_model.p_global_loss_of_load[TS[1], "S2"]) < 1e-09
    end

    #=
    For now,
     capped power is distributed evenly (in terms of percentage of available power)
     among the limitable generators
    TODO? base ratio on prop_cost ? in that case, careful for null and negative prop_cost

    TS: [11h]
    S: [S1]
                      bus 1
                        |
    (limitable) wind_1_1|
    Pmin=0, Pmax=100    |
    Csta=0, Cprop=1     |     load_1
      S1: 50            | S1: 15
                        |
    (limitable) wind_1_2|
    Pmin=0, Pmax=100    |
    Csta=0, Cprop=2     |
      S1: 10            |
                        |
    (limitable) wind_1_3|
    Pmin=0, Pmax=100    |
    Csta=0, Cprop=1     |
      S1: 0             |
    =#
    @testset "energy_market_capping_distribution_in_schedule" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_2", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 2.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_3", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 50.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_2", DateTime("2015-01-01T11:00:00"), "S1", 10.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_3", DateTime("2015-01-01T11:00:00"), "S1", 0.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 15.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict{String, SortedDict{Dates.DateTime, PSCOPF.DecisionFirmness} }(),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                                "wind_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                                "wind_1_3" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE)),
                    )
        # initial generators state : No need because all pmin=0 => ON by default
        generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)
        market = PSCOPF.EnergyMarket()
        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        # Limitable produces to the available level
        @test (0.25 * 50.) ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        @test (0.25 * 10.) ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_2", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.market_schedule, "wind_1_3", TS[1], "S1") <= 1e-09

        # Total capped power
        @test 45. ≈ value(result.limitable_model.p_global_capping[TS[1], "S1"])

        # Capping distributed
        @test 0.75 ≈ 45 / 60 #75% of limitable power will be capped
        @test (0.75 * 50)  ≈ context.market_schedule.capping["wind_1_1", TS[1], "S1"]
        @test (0.75 * 10.) ≈ context.market_schedule.capping["wind_1_2", TS[1], "S1"]
        @test (0.75 * 0.)  ≈ context.market_schedule.capping["wind_1_3", TS[1], "S1"]
    end

    #=
    TS: [11h, 11h15]
    S: [S1,S2]
                        bus 1                   bus 2
                        |                      |
    (limitable) wind_1_1|       "1_2"          |
    Pmin=0, Pmax=200    |                      |
    Csta=0, Cprop=1     |                      |
      S1: 160           |----------------------|
                        |         5            |
                        |                      |
    load                |                      |load
      S1: 180           |                      |  S1: 20

    available prod : 160
    total load : 200
    => need to cut 40 => 20% of load
    cut 20% on each bus :
        bus1: 20% of 180 => 36
        bus2: 20% of 20  => 4
    =#
    @testset "energy_market_lol_distribution_by_bus_in_schedule" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        # Buses
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        PSCOPF.Networks.add_new_bus!(network, "bus_2")
        # Branches
        PSCOPF.Networks.add_new_branch!(network, "branch_1_2", 5.);
        # PTDF
        PSCOPF.Networks.add_ptdf_elt!(network, "branch_1_2", "bus_1", 0.5)
        PSCOPF.Networks.add_ptdf_elt!(network, "branch_1_2", "bus_2", -0.5)
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 200.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 160.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 180.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", DateTime("2015-01-01T11:00:00"), "S1", 20.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict{String, SortedDict{Dates.DateTime, PSCOPF.DecisionFirmness} }(),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE)),
                    )
        # initial generators state : No need because all pmin=0 => ON by default
        generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)
        market = PSCOPF.EnergyMarket()
        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK
        # Limitable produces to the available level
        @test 160. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")

        # Total cut conso
        @test 40. ≈ value(result.lol_model.p_global_loss_of_load[TS[1], "S1"])

        # cut conso distributed on buses
        @test 0.2 ≈ 40 / 200 # 40:loss_of_load, 200:total load
        @test (0.2 * 180) ≈ context.market_schedule.loss_of_load_by_bus["bus_1", TS[1], "S1"]
        @test (0.2 * 20.) ≈ context.market_schedule.loss_of_load_by_bus["bus_2", TS[1], "S1"]
    end

    #TODO : illustrate ptdf effect
    #TODO : test capping and load_shedding transcription in schedule for marketAtFO

end

ext: jl

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
#this file tests whether using a for loop to apply a function of two variables
# f(x, y) to all pairs of entries of two vectors u, v is faster or if
#if using syntactic sugar f.(X, Y) is faster

g = (x, y)->sin(sqrt(x*y+10000))

x = collect(1:1:400)
y = collect(1:1:400)
yadj = y'

println("for loop time")
@time begin
    z = zeros(400, 400)
    for i = 1:400
        for j = 1:400
            z[i, j] = g(x[i], y[j])
        end
    end
end
println("syntactic sugar time")
@time w = g.(x, yadj)
10
println(norm(w-z))

#why is broadcasting slower T_T

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
struct Differential <: Function
    x::Expression
    order::Int
end
Differential(x) = Differential(x,1)

Base.show(io::IO, D::Differential) = print(io,"($(D.x),$(D.order))")
Base.convert(::Type{Expr}, D::Differential) = D

(D::Differential)(x::Operation) = Operation(D, Expression[x])
function (D::Differential)(x::Variable)
    D.x === x             && return Constant(1)
    has_dependent(x, D.x) || return Constant(0)

    x′ = copy(x)
    x′.diff = D
    return x′
end
Base.:(==)(D1::Differential, D2::Differential) = D1.order == D2.order && D1.x == D2.x

function expand_derivatives(O::Operation)
    @. O.args = expand_derivatives(O.args)

    if O.op isa Differential
        D = O.op
        o = O.args[1]
        return simplify_constants(sum(i->Derivative(o,i)*expand_derivatives(D(o.args[i])),1:length(o.args)))
    end

    return O
end
expand_derivatives(x) = x

# Don't specialize on the function here
Derivative(O::Operation, idx) = Derivative(O.op, (O.args...,), Val(idx))

# Pre-defined derivatives
import DiffRules, SpecialFunctions, NaNMath
for (modu, fun, arity) ∈ DiffRules.diffrules()
    for i ∈ 1:arity
        @eval function Derivative(::typeof($modu.$fun), args::NTuple{$arity,Any}, ::Val{$i})
            M, f = $(modu, fun)
            partials = DiffRules.diffrule(M, f, args...)
            dx = @static $arity == 1 ? partials : partials[$i]
            convert(Expression, dx)
        end
    end
end

function count_order(x)
    @assert !(x isa Symbol) "The variable $x must have an order of differentiation that is greater or equal to 1!"
    n = 1
    while !(x.args[1] isa Symbol)
        n = n+1
        x = x.args[1]
    end
    n, x.args[1]
end

function _differential_macro(x)
    ex = Expr(:block)
    lhss = Symbol[]
    x = flatten_expr!(x)
    for di in x
        @assert di isa Expr && di.args[1] == :~ "@Deriv expects a form that looks like `@Deriv D''~t E'~t`"
        lhs = di.args[2]
        rhs = di.args[3]
        order, lhs = count_order(lhs)
        push!(lhss, lhs)
        expr = :($lhs = Differential($rhs, $order))
        push!(ex.args,  expr)
    end
    push!(ex.args, Expr(:tuple, lhss...))
    ex
end

macro Deriv(x...)
    esc(_differential_macro(x))
end

function calculate_jacobian(eqs,vars)
    Expression[Differential(vars[j])(eqs[i]) for i in 1:length(eqs), j in 1:length(vars)]
end

export Differential, Derivative, expand_derivatives, @Deriv, calculate_jacobian


text:
"""
Simplest form of image quantization by turning each pixel to the closest one
in the provided color palette `cs`.
Technically this not a dithering algorithm as the quatization error is not "randomized".
"""
struct ClosestColor <: AbstractDither end

function binarydither!(::ClosestColor, out::GenericGrayImage, img::GenericGrayImage)
    threshold = eltype(img)(0.5)
    return out .= img .> threshold
end

function colordither(
    ::ClosestColor, img::GenericImage, cs::AbstractVector{<:Pixel}, metric::DifferenceMetric
)
    cs = Lab.(cs)
    return map(px -> _closest_color_idx(px, cs, metric), img)
end


text:
using Genie
using Stipple, Stipple.ModelStorage.Sessions
using Stipple.Pages

isempty(Genie.SECRET_TOKEN[]) && Genie.Generator.write_secrets_file(@__DIR__)
Sessions.init()

@reactive mutable struct Person <: ReactiveModel
  name::String = "Bob"
end

function ui(model = init_from_storage(Person))
  on(model.isready) do _
    @show model.name[]
    @show getchannel(model)

    push!(model, channel = getchannel(model))
  end

  page(model, [
    input("", type="text", @bind(:name))
  ])
end

Page("/", view = ui())

up(async = false)


text:
@testset "Kronecker graphs" begin
    P1 = [0.1 0.0 0.0;
          0.3 0.4 0.5;
          1.0 0.7 0.2]

    P2 = [0.1 0.1 0.0;
          0.3 0.4 0.6;
          1.0 0.2 0.2]

    P = P1 ⊗ P2

    P4order = kronecker(P1, 4)

    @testset "Naive sample" begin
          G = naivesample(P)
          @test G isa AbstractSparseMatrix
          @test size(G) == size(P)
          @test G[9, 1]  # probability of 1
          @test !G[2, 5]  # probability of 0
    end

    @testset "Fast sample" begin
          G = fastsample(P)
          @test G isa AbstractSparseMatrix
          @test size(G) == size(P)
          # exact sample does not always generate this, so don't test
          #@test G[9, 1]  # probability of 1
          @test !G[2, 5]  # probability of 0
    end

    @testset "helpers" begin
         @test isprob(P1)
         @test !isprob([2 .01; 0.2 0])
         @test isprob(P)
         @test isprob(P4order)

         indices = sampleindices(P1, 100)
         @test all(maximum(indices) .<= size(P1))
         for ind in indices
               @test P1[ind...] > 0
         end

         indices = sampleindices(P, 100)
         @test all(maximum(indices) .<= size(P))
         for ind in indices
               @test P[ind...] > 0
         end

         indices = sampleindices(P4order, 1000)
         @test all(maximum(indices) .<= size(P4order))
         for ind in indices
               @test P4order[ind...] > 0
         end
   end

end


text:
function dubins_dynamics!(xdot,x,u)
    xdot[1] = u[1]*cos(x[3])
    xdot[2] = u[1]*sin(x[3])
    xdot[3] = u[2]
    xdot
end
n = 3
m = 2

model = Model(dubins_dynamics!,n,m)

x0 = [0.;0.;0.]
xf = [1.;1.;0.]

Q = Array((1e-2)*Diagonal(I,n))
Qf = Array(100.0*Diagonal(I,n))
R = Array((1e-2)*Diagonal(I,m))

tf = 5.0
dt = 0.01

obj = LQRObjective(Q, R, Qf, tf, x0, xf)

solver = TrajectoryOptimization.Solver(model,obj,dt=dt,integration=:rk3)
U = zeros(m,solver.N-1)

results, stats = solve(solver,U)

plot(to_array(results.X)')
plot(to_array(results.X)[1,:],to_array(results.X)[2,:],label="",title="Parallel Park",xlabel="x",ylabel="y")
results.X[end]


text:
module Models

import ACE
ACE.@aceimports
ACE.@baseimports
ACE.@extimports

import ACEatoms: environment, AtomicEnvironment

import JuLIP

import JuLIP: SitePotential,
              cutoff,
              z2i, i2z, numz,
              AbstractCalculator,
              Atoms,
              chemical_symbol, AtomicNumber,
              JVec, JMat,
              energy, forces, virial

import JuLIP.MLIPs: IPBasis

import JuLIP.Potentials: ZList, SZList, zlist

include("twolayer.jl")

end


text:
using Compat
@static if is_windows()
    using WinRPM
    WinRPM.install("gcc", yes=true)
    WinRPM.install("winpthreads-devel", yes=true)
end

@static if is_linux()
    DEPS_PATH = dirname(@__FILE__)
    DOWNLOADS_PATH = joinpath(DEPS_PATH, "downloads")
    INSTALL_PATH = joinpath(DEPS_PATH, "usr", "local")
    BIN_PATH = joinpath(INSTALL_PATH, "bin")
    TAR_FILE = "patchelf-0.9.tar.gz"
    EXTRACTED_DIR = "patchelf-0.9"

    mkpath(DOWNLOADS_PATH)
    mkpath(INSTALL_PATH)

    download("https://nixos.org/releases/patchelf/patchelf-0.9/$TAR_FILE", joinpath(DOWNLOADS_PATH, TAR_FILE))

    cd(DOWNLOADS_PATH) do
        run(`tar -xf $TAR_FILE`)
    end

    cd(joinpath(DOWNLOADS_PATH, EXTRACTED_DIR)) do
        run(`./configure --prefix="$INSTALL_PATH"`)
        run(`make install`)
        @assert ispath(joinpath(INSTALL_PATH, "bin", "patchelf"))
    end
end


text:
module Benchmark

import MathOptInterface
import SCS
import GLPK
import JSON
import Printf
import Random

using SparseArrays
using TimerOutputs

const MOI = MathOptInterface

struct PMedianData
    num_facilities::Int
    num_customers::Int
    num_locations::Int
    customer_locations::Vector{Float64}
end

# This is the LP relaxation.
function generate_moi_problem(model, data::PMedianData)
    NL = data.num_locations
    NC = data.num_customers

    ###
    ### 0 <= facility_variables <= 1
    ###

    facility_variables = MOI.add_variables(model, NL)

    for v in facility_variables
        MOI.add_constraint(model, MOI.SingleVariable(v), MOI.Interval(0.0, 1.0))
    end

    ###
    ### assignment_variables >= 0
    ###

    assignment_variables = reshape(MOI.add_variables(model, NC * NL), NC, NL)
    for v in assignment_variables
        MOI.add_constraint(model, MOI.SingleVariable(v), MOI.GreaterThan(0.0))
        # "Less than 1.0" constraint is redundant.
    end

    ###
    ### Objective function
    ###

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            [
                MOI.ScalarAffineTerm(
                    abs(data.customer_locations[i] - j),
                    assignment_variables[i, j]
                )
                for i in 1:NC for j in 1:NL
            ],
            0.0,
        ),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    ###
    ### assignment_variables[i, j] <= facility_variables[j]
    ###

    for i in 1:NC, j in 1:NL
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(
                [
                    MOI.ScalarAffineTerm(1.0, assignment_variables[i, j]),
                    MOI.ScalarAffineTerm(-1.0, facility_variables[j])
                ],
                0.0,
            ),
            MOI.LessThan(0.0),
        )
    end

    ###
    ### sum_j assignment_variables[i, j] = 1
    ###

    for i in 1:NC
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(
                [
                    MOI.ScalarAffineTerm(1.0, assignment_variables[i, j])
                    for j in 1:NL
                ],
                0.0,
            ),
            MOI.EqualTo(1.0),
        )
    end

    ###
    ### sum_j facility_variables[j] == num_facilities
    ###

    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.(1.0, facility_variables),
            0.0,
        ),
        MOI.EqualTo{Float64}(data.num_facilities),
    )

    return assignment_variables, facility_variables
end

function generate_moi_problem_vector(model, data::PMedianData)
    NL = data.num_locations
    NC = data.num_customers

    ###
    ### facility_variables    ∈ R₊
    ### facility_variables -1 ∈ R₋
    ###

    facility_variables = MOI.add_variables(model, NL)
    MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [
                MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(1.0, x))
                for (i, x) in enumerate(facility_variables)
            ],
            zeros(NL),
        ),
        MOI.Nonnegatives(NL),
    )

    MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [
                MOI.VectorAffineTerm(i, MOI.ScalarAffineTerm(1.0, x))
                for (i, x) in enumerate(facility_variables)
            ],
            -ones(NL),
        ),
        MOI.Nonpositives(NL),
    )

    ###
    ### assignment_variables ∈ R₊
    ###

    assignment_variables = reshape(MOI.add_variables(model, NC * NL), NC, NL)
    MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            [
                MOI.VectorAffineTerm(
                    i, MOI.ScalarAffineTerm(1.0, assignment_variables[i])
                )
                for i in eachindex(assignment_variables)
            ],
            zeros(NC * NL),
        ),
        MOI.Nonnegatives(NC * NL),
    )
    # "Less than 1.0" constraint is redundant.

    ###
    ### Objective function
    ###

    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            [
                MOI.ScalarAffineTerm(
                    abs(data.customer_locations[i] - j),
                    assignment_variables[i, j],
                )
                for i in 1:NC for j in 1:NL
            ],
            0.0,
        )
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    ###
    ### assignment_variables[i,j] <= facility_variables[j]
    ###     ⟺
    ### assignment_variables[i,j] - facility_variables[j] ∈ R₋
    ###

    terms = Vector{MOI.VectorAffineTerm{Float64}}(undef, 2 * NC * NL)
    row = 1
    for i in 1:NC, j in 1:NL
        terms[2 * row - 1] = MOI.VectorAffineTerm(
            row, MOI.ScalarAffineTerm(1.0, assignment_variables[i, j])
        )
        terms[2 * row] = MOI.VectorAffineTerm(
            row, MOI.ScalarAffineTerm(-1.0, facility_variables[j])
        )
        row += 1
    end

    MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(terms, zeros(NC * NL)),
        MOI.Nonpositives(NC * NL),
    )

    ###
    ### sum_j assignment_variables[i,j] == 1
    ###     ⟺
    ### sum_j assignment_variables[i,j] - 1 ∈ 0
    ###

    terms_2 = Vector{MOI.VectorAffineTerm{Float64}}(undef, NC * NL)
    k = 1
    for i in 1:NC, j in 1:NL
        terms_2[k] = MOI.VectorAffineTerm(
            i, MOI.ScalarAffineTerm(1.0, assignment_variables[i, j])
        )
        k += 1
    end
    MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(terms_2, -ones(NC)),
        MOI.Zeros(NC),
    )

    ###
    ### sum_j facility_variables[j] - num_facilities ∈ 0
    ###

    MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                1, MOI.ScalarAffineTerm.(1.0, facility_variables)
            ),
            Float64[-data.num_facilities],
        ),
        MOI.Zeros(1),
    )

    return assignment_variables, facility_variables
end

function generate_glpk_problem(prob, data::PMedianData)
    NL = data.num_locations
    NC = data.num_customers

    facility_variables = [GLPK.glp_add_cols(prob, 1) for i in 1:NL]

    for j in 1:NL
        GLPK.glp_set_col_bnds(prob, facility_variables[j], GLPK.GLP_DB, 0.0, 1.0)
    end

    assignment_variables = [
        GLPK.glp_add_cols(prob, 1)
        for i in 1:data.num_customers, j in 1:data.num_locations
    ]

    for i in 1:data.num_customers, j in 1:data.num_locations
        GLPK.glp_set_col_bnds(prob, assignment_variables[i, j], GLPK.GLP_DB, 0.0, 1.0)
        GLPK.glp_set_obj_coef(prob, assignment_variables[i, j], abs(data.customer_locations[i] - j))
    end
    GLPK.glp_set_obj_dir(prob, GLPK.GLP_MIN)

    I = Cint[0] # Extra padding because GLPK starts reading from the second element
    J = Cint[0]
    V = Float64[0.0]
    for i in 1:data.num_customers, j in 1:data.num_locations
        # assignment_variables[i,j] <= facility_variables[j]
        index = GLPK.glp_add_rows(prob, 1)
        GLPK.glp_set_row_bnds(prob, index, GLPK.GLP_UP, 0.0, 0.0)
        push!(I, index)
        push!(J, assignment_variables[i, j])
        push!(V, 1.0)
        push!(I, index)
        push!(J, facility_variables[j])
        push!(V, -1.0)
    end
    for i in 1:data.num_customers
        # sum_j assignment_variables[i,j] = 1
        index = GLPK.glp_add_rows(prob, 1)
        GLPK.glp_set_row_bnds(prob, index, GLPK.GLP_FX, 1.0, 1.0)
        for j in 1:data.num_locations
            push!(I, index)
            push!(J, assignment_variables[i, j])
            push!(V, 1.0)
        end
    end
    # sum_j facility_variables[j] = num_facilities
    index = GLPK.glp_add_rows(prob, 1)
    GLPK.glp_set_row_bnds(prob, index, GLPK.GLP_FX, data.num_facilities, data.num_facilities)
    for j in 1:data.num_locations
        push!(I, index)
        push!(J, facility_variables[j])
        push!(V, 1.0)
    end
    GLPK.glp_load_matrix(prob, length(I) - 1, I, J, V)
    return
end

function solve_glpk_direct(data::PMedianData; time_limit_sec=Inf)
    @timeit "GLPK direct" begin
        prob = GLPK.glp_create_prob()
        param = GLPK.glp_smcp()
        GLPK.glp_init_smcp(param)
        param.msg_lev = GLPK.GLP_MSG_ERR
        if isfinite(time_limit_sec)
            param.tm_lim = 1000 * time_limit_sec
        end
        @timeit "generate" generate_glpk_problem(prob, data)
        @timeit "solve" GLPK.glp_simplex(prob, param)
        objval = GLPK.glp_get_obj_val(prob)
        GLPK.glp_delete_prob(prob)
    end
    return objval
end

function generate_scs_problem(data::PMedianData)
    # Equalities must be ordered first.
    num_variables = 0
    facility_variables = [num_variables += 1 for i in 1:data.num_locations]
    assignment_variables = [num_variables += 1 for i in 1:data.num_customers, j in 1:data.num_locations]

    b = Float64[]
    I = Int[]
    J = Int[]
    V = Float64[]

    num_constraints = 0
    for i in 1:data.num_customers
        # sum_j assignment_variables[i,j] = 1
        num_constraints += 1
        push!(b, 1.0)
        for j in 1:data.num_locations
            push!(I, num_constraints)
            push!(J, assignment_variables[i, j])
            push!(V, 1.0)
        end
    end
    # sum_j facility_variables[j] = num_facilities
    num_constraints += 1
    push!(b, data.num_facilities)
    for j in 1:data.num_locations
        push!(I, num_constraints)
        push!(J, facility_variables[j])
        push!(V, 1.0)
    end

    # Now inequality constraints in the form b - a'x >= 0.

    for i in 1:data.num_customers, j in 1:data.num_locations
        # assignment_variables[i,j] <= facility_variables[j]
        num_constraints += 1
        push!(b, 0.0)
        push!(I, num_constraints)
        push!(J, assignment_variables[i, j])
        push!(V, 1.0)
        push!(I, num_constraints)
        push!(J, facility_variables[j])
        push!(V, -1.0)
    end

    # All variables are nonnegative.
    for v in 1:num_variables
        num_constraints += 1
        push!(b, 0.0)
        push!(I, num_constraints)
        push!(J, v)
        push!(V, -1.0)
    end

    # facility_variables <= 1
    for j in 1:data.num_locations
        num_constraints += 1
        push!(b, 1.0)
        push!(I, num_constraints)
        push!(J, facility_variables[j])
        push!(V, 1.0)
    end

    c = zeros(num_variables)
    for i in 1:data.num_customers, j in 1:data.num_locations
        c[assignment_variables[i, j]] = abs(data.customer_locations[i] - j)
    end
    num_zero_cones = data.num_customers + 1
    num_positive_orthants = num_constraints - num_zero_cones

    A = sparse(I, J, V, num_constraints, num_variables)

    return num_constraints, num_variables, A, b, c, num_zero_cones, num_positive_orthants,
        Int[], Int[], 0, 0, Float64[]
end

function solve_scs_direct(data::PMedianData; max_iters)
    @timeit "SCS direct" begin
        @timeit "generate" scs_prob = generate_scs_problem(data)
        @timeit "solve" solution = SCS.SCS_solve(
            SCS.IndirectSolver,
            scs_prob...;
            max_iters=max_iters,
            acceleration_lookback=0,
            verbose=0,
        )
    end
    objval = scs_prob[5]' * solution.x
    return objval
end

function solve_moi(data::PMedianData, optimizer; vector_version, params)
    cache = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        optimizer()
    )
    model = MOI.Bridges.full_bridge_optimizer(cache, Float64)
    # resetting optimizer is necessary to use copy_to to transfer all the data
    # from caching optimizer to GLPK in a single batch. If that is not called
    # constraints are passed one by one during model generation phase.
    MOI.Utilities.reset_optimizer(cache)
    for (param, value) in params
        MOI.set(model, param, value)
    end
    @timeit "generate" x, y = if vector_version
        generate_moi_problem_vector(model, data)
    else
        generate_moi_problem(model, data)
    end
    @timeit "solve" MOI.optimize!(model)
    return MOI.get(model, MOI.ObjectiveValue())
end

function solve_glpk_moi(data::PMedianData; vector_version, time_limit_sec=Inf)
    params = []
    if isfinite(time_limit_sec)
        push!(params, (MOI.TimeLimitSec(), time_limit_sec))
    end
    s_type = vector_version ? "vector" : "scalar"
    @timeit(
        "GLPK MOI $(s_type)",
        solve_moi(
            data, GLPK.Optimizer; vector_version=vector_version, params=params
        )
    )
end

function solve_scs_moi(data::PMedianData; vector_version, max_iters::Int)
    params = Tuple{MOI.RawParameter, Int}[
        (MOI.RawParameter("max_iters"), max_iters),
        (MOI.RawParameter("verbose"), 0),
        (MOI.RawParameter("acceleration_lookback"), 0)
    ]
    s_type = vector_version ? "vector" : "scalar"
    @timeit(
        "SCS MOI $(s_type)",
        solve_moi(
            data, SCS.Optimizer; vector_version=vector_version, params=params
        )
    )
end

function run_benchmark(;
    num_facilities, num_customers, num_locations, time_limit_sec, max_iters
)
    Random.seed!(10)
    reset_timer!()
    data = PMedianData(num_facilities, num_customers, num_locations, rand(num_customers) .* num_locations)
    GC.gc()
    glpk_moi_vector_obj = solve_glpk_moi(data, vector_version=true, time_limit_sec=time_limit_sec)
    @show glpk_moi_vector_obj
    GC.gc()
    glpk_moi_scalar_obj = solve_glpk_moi(data, vector_version=false, time_limit_sec=time_limit_sec)
    @show glpk_moi_scalar_obj
    GC.gc()
    glpk_direct_obj = solve_glpk_direct(data, time_limit_sec=time_limit_sec)
    @show glpk_direct_obj
    GC.gc()
    scs_moi_vector_obj = solve_scs_moi(data, vector_version=true, max_iters=max_iters)
    @show scs_moi_vector_obj
    GC.gc()
    scs_moi_scalar_obj = solve_scs_moi(data, vector_version=false, max_iters=max_iters)
    @show scs_moi_scalar_obj
    GC.gc()
    scs_direct_obj = solve_scs_direct(data, max_iters=max_iters)
    @show scs_direct_obj

    print_timer()
    println()
end

function run_paper_benchmark(L)
    to = TimerOutputs.get_defaulttimer()
    s = ""
    if isfile("cvxpy.json")
        cvxpy = JSON.parsefile("cvxpy.json")
    else
        cvxpy = nothing
        println("WARNING: cvxpy.json not found.")
    end
    for l in L
        row_starts = Dict(
            "generate" => "\t\t\\multirow{3}{*}{$(l)} & generate",
            "solve" => "\t\t& load",
            "total" => "\t\t& total",
        )
        row_computation = Dict(
            "generate" => (to, key) -> round(
                TimerOutputs.time(to[key]["generate"]) / 1e9; digits = 2
            ),
            "solve" => (to, key) -> round(
                TimerOutputs.time(to[key]["solve"]) / 1e9; digits = 2
            ),
            "total" => (to, key) -> round(
                (
                    TimerOutputs.time(to[key]["generate"]) +
                    TimerOutputs.time(to[key]["solve"])
                ) / 1e9;
                digits = 2
            ),
        )
        run_benchmark(
            num_facilities = 100,
            num_customers = 100,
            num_locations = l,
            time_limit_sec = 0.001,
            max_iters = 1,
        )
        for row in ["generate", "solve", "total"]
            s *= row_starts[row]
            for key in [
                "GLPK MOI scalar",
                "GLPK MOI vector",
                "GLPK direct",
                "GLPK CVXPY",
                "SCS MOI scalar",
                "SCS MOI vector",
                "SCS direct",
                "SCS CVXPY",
            ]
                t = if key == "GLPK CVXPY"
                    cvxpy !== nothing ? cvxpy["glpk"]["$(l)"][row] : NaN
                elseif key == "SCS CVXPY"
                    cvxpy !== nothing ? cvxpy["scs"]["$(l)"][row] : NaN
                else
                    row_computation[row](to, key)
                end
                s *= Printf.@sprintf " & %.2f" t
            end
            s *= " \\\\\n"
        end
    end
    println(s)
    return s
end

end # module

max_stars_repo_name: jump-dev/MOIPaperBenchmarks

max_issues_repo_name: jump-dev/MOIPaperBenchmarks

max_forks_repo_name: jump-dev/MOIPaperBenchmarks


text:
GrowthPowerDict13 = Dict([("δ3", 0), ("η3", 1) , ("2tr[KK2]", 0), ("δη2", 1),
                           ("2ηη2", 2), ("2KK2∥", 0), ("Otd", 0), ("δΠ2∥", 0),
                           ("ηΠ2∥", 1), ("Π2K∥", 0), ("s^k∂kΠ2∥", 0),
                           ("u2∥∂∥δ", 1), ("u2∥∂∥η", 2), ("u∥∂∥η2", 2),
                           ("u∥∂∥Π2∥", 1), ("Π3∥", 0)])

@memoize function P_13ln(Pk,ℓ,n)
    r, ξℓn = ξ(Pk,ℓ,n)

    function get_int_func(f)
        Spline1D(r, f./r, bc="extrapolate")
    end
    if n == 0
        q = .5
    elseif n == 1
        q = .1
    elseif n == -1
        q = .5
    elseif n == -2
        q = 1.1
    elseif (ℓ == 2) && (n == 2)
        q = -.8
    elseif n == 2
        q = -1
    end
    k, mcP = xicalc(get_int_func(ξℓn), ℓ, 0, kmin=rmini, kmax=rmaxi, r0=r0i, N=N, q=q)
    return k, 2π^2 .* mcP
end

@memoize function I1(Pk)
    k, mcP1m1 = P_13ln(Pk,1,-1)
    _, mcP3m1 = P_13ln(Pk,3,-1)
    mcI1 = @. (2k^3/5)*(mcP1m1 - mcP3m1)
    return k, mcI1
end

@memoize function I2(Pk)
    k, mcP00 = P_13ln(Pk,0,0)
    _, mcP20 = P_13ln(Pk,2,0)
    mcI2 = @. (2k^2/3)*(mcP00 - mcP20)
    return k, mcI2
end

@memoize function I3(Pk)
    k, mcP02 = P_13ln(Pk,0,2)
    _, mcP22 = P_13ln(Pk,2,2)
    mcI3 = @. (2/3)*(mcP02 - mcP22)
    return k, mcI3
end

@memoize function I4(Pk)
    k, mcP02 = P_13ln(Pk,0,2)
    _, mcP22 = P_13ln(Pk,2,2)
    _, mcP42 = P_13ln(Pk,4,2)
    mcI4 = @. (2/15)*mcP02 + (2/21)*mcP22 - (8/35)*mcP42
    return k, mcI4
end

@memoize function I5(Pk)
    k, mcP02 = P_13ln(Pk,0,2)
    _, mcP22 = P_13ln(Pk,2,2)
    _, mcP42 = P_13ln(Pk,4,2)
    _, mcP62 = P_13ln(Pk,6,2)
    mcI5 = @. (2/35)*mcP02 + (2/21)*mcP22 - (32/385)*mcP42 - (16/231)*mcP62
    return k, mcI5
end

function fO(Pk,O::String,μs)
    m = M13(O)
    k, i1 = I1(Pk)
    _, i2 = I2(Pk)
    _, i3 = I3(Pk)
    _, i4 = I4(Pk)
    _, i5 = I5(Pk)
    fout = Array{Float64}(undef,length(k),length(μs))
    for j = 1:length(μs)
        for i = 1:length(k)
            fout[i,j] = ([1 μs[j]^2 μs[j]^4] * m * [i1[i], i2[i], i3[i], i4[i], i5[i]])[1]
        end
    end
    return k, fout
end

function P13ggs(Pk, μs; BiasOperatorDict13, f)
    p13 = zeros(Float64, N, length(μs))
    kout = Array{Float64}(undef, N)
    i = 1
    b1 = BiasOperatorDict13["δ3"]
    bη = BiasOperatorDict13["η3"]
    local k, fo
    for op in keys(BiasOperatorDict13)
        try
            k, fo = fO(Pk,op,μs)
        catch
            continue
        end

        if i == 1
            kout .= k
            i += 1
        end
        bo = BiasOperatorDict13[op]
        nf = GrowthPowerDict13[op]
        for j = 1:length(μs)
            p13[:,j] .+= (b1-bη*f*μs[j]^2) * bo * f^nf .* fo[:,j] .* Pk.(k)
        end
    end
    return kout, p13
end

function P13lggs(Pk, l; BiasOperatorDict13, f, Cln)
    if l == 8
        return [0.0], [0.0]
    end
    k, i1 = I1(Pk)
    _, i2 = I2(Pk)
    _, i3 = I3(Pk)
    _, i4 = I4(Pk)
    _, i5 = I5(Pk)
    In = [i1,i2,i3,i4,i5]
    Pkvals = Pk.(k)
    p13l = zeros(Float64, length(k))
    for n = 1:5
        p13l .+= Cln[l+1,n] .* In[n] .* Pkvals
    end
    return k, p13l
end

function P13ggsmultipole(Pk, μs; BiasOperatorDict13, f)
    L(n,x) = Legendre(n,x)
    Cln = generateCln(BiasOperatorDict13, f)
    Pl = l -> P13lggs(Pk,l, BiasOperatorDict13=BiasOperatorDict13, f=f, Cln=Cln)
    k, P0 = Pl(0)
    _, P2 = Pl(2)
    _, P4 = Pl(4)
    _, P6 = Pl(6)
    p13mp = Array{Float64}(undef, length(k), length(μs))
    for i = 1:length(μs)
        @. p13mp[:,i] = L(0,μs[i])*P0 + L(2,μs[i])*P2 + L(4,μs[i])*P4 + L(6,μs[i])*P6
    end
    return k, p13mp
end

max_stars_repo_name: JosephTomlinson/GeneralBiasPk

max_issues_repo_name: JosephTomlinson/GeneralBiasPk

max_forks_repo_name: JosephTomlinson/GeneralBiasPk


text:
if false
using Test
using GraphPlot
using Neo4j
using DistributedFactorGraphs
using Pkg
using Dates
using UUIDs

include("testBlocks.jl")

testDFGAPI = CloudGraphsDFG
end
# Build a basic graph.

testDFGAPI = LightDFG
dfg = testDFGAPI{NoSolverParams}()
#add types for softtypes
# struct TestInferenceVariable1 <: InferenceVariable end
# v1 = DFGVariable(:a, TestInferenceVariable1())
# v2 = DFGVariable(:b, TestInferenceVariable1())
# f1 = DFGFactor{TestFunctorInferenceType1}(:f1, [:a,:b]);
# #add tags for filters
# union!(v1.tags, [:VARIABLE, :POSE])
# union!(v2.tags, [:VARIABLE, :LANDMARK])
# union!(f1.tags, [:FACTOR])
# # @testset "Creating Graphs" begin

var1, var2, var3, vorphan, v1_tags = DFGVariableSCA()
fac0, fac1, fac2 = DFGFactorSCA()

addVariable!(dfg, var1)
addVariable!(dfg, var2)
addFactor!(dfg, fac1)

# Stores to test
testStores = [InMemoryDataStore(), FileDataStore("/tmp/dfgFilestore")]

if false
testStore = testStores[1]
testStore = testStores[2]
end

for testStore in testStores
    # Create a data store and a dataset
ds = testStore
# ds = FileDataStore("/tmp/filestore")
dataset = rand(UInt8, 1000)
dataset2 = rand(UInt8, 1000)

entry1 = GeneralDataEntry(dfg, var1, :test1)
# Set it in the store
@test addDataBlob!(ds, entry1, dataset) == dataset
@test getDataBlob(ds, entry1) == dataset
# Now add it to the variable
@test addDataEntry!(var1, entry1) == entry1
@test entry1 in getDataEntries(var1)
# Update test
copyEntry = deepcopy(entry1)
sleep(0.1)
@test updateDataBlob!(ds, entry1, dataset2) == dataset2
# Data updated?
@test getDataBlob(ds, entry1) == dataset2
# Timestamp updated?
@test entry1.lastUpdatedTimestamp > copyEntry.lastUpdatedTimestamp
# Delete data
@test deleteDataBlob!(ds, entry1) == dataset2
# Delete entry
@test deleteDataEntry!(var1, entry1) == var1

# @test addData!()
# gde_td1 = DFG.addData!(dfg, :a, ds, :td1, "text/plain", rand(UInt8,10))
# DFG.getDataEntryBlob(dfg, :a, ds, :td1)
#
#
#
# @test getDataBlob(ds, entry1) == dataset2

#TODO
    #listStoreEntries(ds)
end


text:
# This has the same job as the LogBinner Compressor
mutable struct EPCompressor{T}
    # One value per input argument
    values::Vector{T}
    switch::Bool
end

struct ErrorPropagator{T, N}
    compressors::NTuple{N, EPCompressor{T}}
    # ∑x
    sums1D::NTuple{N, Vector{T}}
    #∑xy
    sums2D::NTuple{N, Matrix{T}}

    count::Vector{Int64}
end

# Overload some basic Base functions
Base.eltype(ep::ErrorPropagator{T,N}) where {T,N} = T
Base.length(ep::ErrorPropagator) = ep.count[1]
Base.ndims(ep::ErrorPropagator{T,N}) where {T,N} = ndims(eltype(ep))
Base.isempty(ep::ErrorPropagator) = length(ep) == 0

function _print_header(io::IO, ep::ErrorPropagator{T,N}) where {T,N}
    print(io, "ErrorPropagator{$(T),$(N)}")
    nothing
end

function _println_body(io::IO, ep::ErrorPropagator{T,N}) where {T,N}
    n = length(ep)
    println(io)
    print(io, "| Count: ", n)
    if n > 0 && ndims(ep) == 0
        print(io, "\n| Means: ", round.(means(ep), digits=5))
        print(io, "\n| StdErrors: ", round.(std_errors(ep), digits=5))
    end
    nothing
end

# short version (shows up in arrays etc.)
function Base.show(io::IO, ep::ErrorPropagator{T,N}) where {T,N}
    print(io, "ErrorPropagator{$(T),$(N)}()")
end
# verbose version (shows up in the REPL)
function Base.show(io::IO, m::MIME"text/plain", ep::ErrorPropagator)
    (_print_header(io, ep); _println_body(io, ep))
end

"""
    empty!(ep::ErrorPropagator)

Clear the error propagator, i.e. reset it to its inital state.
"""
function Base.empty!(ep::ErrorPropagator)
    !isempty(ep) || return

    # get zero
    z = zero(eltype(eltype(ep)))

    # reset sums to all zeros
    @inbounds for i in eachindex(ep.sums1D)
        if ndims(ep) == 0 # compile time
            # numbers
            ep.sums1D[i] .= z
            ep.sums2D[i] .= z
        else
            # arrays
            for j in eachindex(ep.sums1D[i])
                ep.sums1D[i][j] .= z
            end
            for j in eachindex(ep.sums2D[i])
                ep.sums2D[i][j] .= z
            end
        end
    end

    # reset counts
    fill!(ep.count, 0)

    # reset compressors
    for c in ep.compressors
        if ndims(ep) == 0 # compile time
            c.values .= z
        else
            for i in eachindex(c.values)
                c.values[i] .= z
            end
        end
        c.switch = false
    end

    nothing
end

"""
    capacity(ep)

Capacity of the error propagator, i.e. how many values can be handled before
overflowing.
"""
capacity(ep::ErrorPropagator{T,N}) where {T,N} = _nlvls2capacity(N)
nlevels(ep::ErrorPropagator{T,N}) where {T,N} = N

"""
    ErrorPropagator([::Type{T}; N_args, capacity::Int])

Creates an `ErrorPropagator` which can handle (at least) `capacity` many values
for each of `N_args` inputs of type `T`.

The default is `T = Float64`, `N_args = 2` and `capacity = 2^32-1 ≈ 4e9`.
"""
function ErrorPropagator(::Type{T} = Float64; N_args::Int64 = 2, kw...) where T
    ErrorPropagator([zero(T) for _ in 1:N_args]...; kw...)
end

"""
    ErrorPropagator(zero_element::T...[; capacity::Int])

Creates a new `ErrorPropagator` which can take (at least) `capacity` many
values for each input of type `T`. The type and the size are inherited from the
given `zero_element`s, which must exclusively contain zeros.

Values can be added using `push!` and `append!`.

---

    ErrorPropagator(timeseries::AbstractVector{T}...)

Creates a new `ErrorPropagator` and adds all elements from each given timeseries.
"""
function ErrorPropagator(
            xs::T...;
            capacity::Int64 = _nlvls2capacity(32)
        ) where {T <: Union{Number, AbstractArray}}

    # check keyword args
    capacity <= 0 && throw(ArgumentError("`capacity` must be finite and positive."))

    # got_timeseries = didn't receive a zero && is a vector
    # got_timeseries = count(!iszero, x[1]) > 0 && ndims(T) == 1
    got_timeseries = any(count(!iszero, x) > 0 && ndims(T) == 1 for x in xs)

    if got_timeseries
        # xs[i] = time series
        N = _capacity2nlvls(length(xs[1]))
        S = _sum_type_heuristic(eltype(T), ndims(xs[1][1]))
        el = zero(xs[1][1])
    else
        # xs[i] = zero_element
        N = _capacity2nlvls(capacity)
        S = _sum_type_heuristic(T, ndims(xs[1]))
        el = xs[1]
    end

    ep = ErrorPropagator{S, N}(
        tuple([
            EPCompressor{S}([copy(el) for _ in 1:length(xs)], false)
            for __ in 1:N
        ]...),
        tuple([
            [copy(el) for _ in 1:length(xs)]
            for __ in 1:N
        ]...),
        tuple([
            [copy(el) for _ in 1:length(xs), __ in 1:length(xs)]
            for __ in 1:N
        ]...),
        zeros(Int64, N)
    )

    got_timeseries && append!(ep, xs...)

    return ep
end

# TODO typing?
"""
    append!(ErrorPropagator, values...)

Adds multiple arrays of values to the Error Propagator by `push!`ing each element.
"""
function Base.append!(ep::ErrorPropagator, values::AbstractArray...)
    # TODO does this perform well?
    @inbounds for args in zip(values...)
        _push!(ep, 1, args)
    end
    nothing
end

"""
    push!(ep::ErrorPropagator, args...)

Pushes a set of argumentes `args` into a given error propagator. Note that the
number of arguments must match the initially defined `N_arguments` of the
error propagator.
"""
@inline function Base.push!(ep::ErrorPropagator{T, N}, args::T...) where {T, N}
    # @boundscheck allows you to define a boundscheck that can be omitted with
    # @inbounds
    # but that requires an @inline
    @boundscheck begin
        (length(args) == length(ep.sums1D[1])) || throw(DimensionMismatch(
            "Number of arguments $(length(args)) does not match the number " *
            "of arguments accepted by the error propagator (N_arguments = " *
            "$(length(ep.sums1D[1])))"
        ))
    end
    @inbounds _push!(ep, 1, args)
end

@inline _mult(x, y) = x * y
@inline _mult(x::Complex, y::Complex) = x * conj(y)
@inline _mult(x::AbstractArray, y::AbstractArray) = _mult.(x, y)

# recursion, back-end function
@inline function _push!(ep::ErrorPropagator{T, N}, lvl::Int64, args) where {T, N}
    @boundscheck begin
        (length(args) == length(ep.sums1D[1])) || throw(DimensionMismatch(
            "Number of arguments $(length(args)) does not match the number " *
            "of arguments accepted by the error propagator (N_arguments = " *
            "$(length(ep.sums1D[1])))"
        ))
    end
    C = ep.compressors[lvl]

    # any value propagating through this function is new to lvl. Therefore we
    # add it to the sums. Note that values pushed to the output arrays are not
    # added here until the array drops to the next level. (New compressors are
    # added)
    @simd for i in eachindex(ep.sums1D[lvl])
        ep.sums1D[lvl][i] += args[i]
        for j in eachindex(ep.sums1D[lvl])
            ep.sums2D[lvl][i, j] += _mult(args[i], args[j])
        end
    end
    ep.count[lvl] += 1

    if !C.switch
        # Compressor has space -> save value
        C.values .= args
        C.switch = true
        return nothing
    else
        # Do averaging
        if lvl == N
            # No more propagation possible -> throw error
            throw(OverflowError("The Error Propagator has exceeded its maximum capacity."))
        else
            # propagate to next lvl
            C.switch = false
            @inbounds _push!(ep, lvl+1, 0.5 * (C.values .+ args))
            return nothing
        end
    end

    return nothing
end


text:
struct Loop{T} <: FiniteLengthAnimation{T}
    animation::FiniteLengthAnimation{T}
    start::Float64
    gap::Float64
    repetitions::Int

    function Loop(a::FiniteLengthAnimation{T}, tstart::Real, gap::Real, repetitions::Int) where T
        if gap < 0
            error("Gap cannot be smaller than zero, but is $gap")
        end

        if repetitions < 1
            error("Repetitions must be a positive integer, not $repetitions")
        end

        # animation should start at zero
        if start(a) != 0
            error("First timestamp of an animation in a loop must be 0, not $(start(a)).")
        end

        new{T}(a, tstart, gap, repetitions)
    end
end

function at(l::Loop, t::Real)
    s = stop(l.animation)
    duration_gapped = s + l.gap

    n_repetitions_done, t_within_animation = divrem(t - l.start, duration_gapped)
    if n_repetitions_done >= l.repetitions
        at(l.animation, t - l.start)
    else
        at(l.animation, t_within_animation)
    end
end

(l::Loop)(t::Real) = at(l, t)

duration(l::Loop) = (duration(l.animation) + l.gap) * l.repetitions


text:
#=
File with utility functions that don't really fit anywhere else. May change
locations/name in the future.

Author: Ivan Alejandro Moreno Soto
Last updated: 2019/July/15
=#

module Misc

export InvalidCodewordError,
    string2bits,
    bits2string

struct InvalidCodewordError <: Exception end

"""
    string2bits(plaintext, k)

Takes a plaintext and returns its binary representation with as much padding as
needed to make its length divisible by k.
"""
function string2bits(plaintext, k)
    plaintextBytes = Vector{UInt8}(plaintext)
    bits = Array{UInt32, 1}([digit for number in digits.(plaintextBytes, base=2, pad=8) for digit in number])

    while length(bits) % k != 0
        push!(bits, false)
    end

    return bits
end

"""
    bits2string(msg)

Returns the string that msg represents.
"""
function bits2string(msg)
    decomposed_msg = []
    # Taking 8 bits segments.
    for i in 1:div(length(msg), 8)
        current_character = msg[(8*(i-1) + 1):(8*i)]

        # Stopping if padding is found.
        if !(1 in current_character)
            break
        end

        current_string = ""
        for c in current_character
            current_string = current_string * string(c)
        end
        current_string = reverse(current_string)
        push!(decomposed_msg, Char(parse(Int, current_string, base = 2)))
    end

    # Concatenating every character in a single string.
    text = ""
    for c in decomposed_msg
        text = text * c
    end
    return text
end

end

max_stars_repo_name: rexemin/McEliece

max_issues_repo_name: rexemin/McEliece

max_forks_repo_name: rexemin/McEliece


text:
using LinearAlgebra
using Distributions, Random

#softmax(x) = exp.(x) / sum(exp.(x))
function softmax(x)
    x = clamp.(x, -32., 32.)
    return exp.(x) / sum(exp.(x))
end

function sample_discrete(p, rng)
    n = length(p)
    i = 1
    c = p[1]
    u = rand(rng)
    while c < u && i < n
        c += p[i += 1]
    end

    return i
end

mutable struct LinearSoftmaxPolicy{TP, TB} <: AbstractPolicy where {TP<:Real,TB<:AbstractFuncApprox}
    f::LinearFunction{TP, TB}
    action::Int
    probs::Array{TP,1}
    function LinearSoftmaxPolicy(state_dim::Int, num_actions::Int)
        f = LinearFunction(state_dim, num_actions)
        new{eltype(f.θ),typeof(f.ϕ)}(f, 0, zeros(eltype(f.θ), num_actions))
    end
    function LinearSoftmaxPolicy(::Type{T}, state_dim::Int, num_actions::Int) where {T<:Real}
        f = LinearFunction(T, state_dim, num_actions)
        new{T,typeof(f.ϕ)}(f, 0, zeros(T, num_actions))
    end
    function LinearSoftmaxPolicy(f::LinearFunction{T}) where {T<:Real}
        new{T,typeof(f.ϕ)}(clone(f), 0, zeros(T, get_num_outputs(f)))
    end
end

function get_num_params(π::LinearSoftmaxPolicy)::Int
    return get_num_params(π.f)
end

function get_action_probabilities!(π::LinearSoftmaxPolicy, x)
    call!(π.probs,π.f,x)
    π.probs .= softmax(π.probs)
end

function get_action_probabilities(π::LinearSoftmaxPolicy{T}, x)::Array{T,1} where {T<:Real}
    get_action_probabilities!(π, x)
    return π.probs
end

function get_action!(π::LinearSoftmaxPolicy{T}, x, rng::AbstractRNG)::T where {T<:Real}
    get_action_probabilities!(π, x)
    π.action = sample_discrete(π.probs, rng)
    logp = log(π.probs[π.action])
    return logp
end

function gradient_logp!(grad::Array{T}, π::LinearSoftmaxPolicy{T, IdentityBasis}, x::Int, action::Int)::T where {T<:Real}
    grad = reshape(grad, size(get_params(π.f)))
    fill!(grad, 0.)

    get_action_probabilities!(π, x)
    # probs = -π.probs
    logp = log(π.probs[action])
    # probs[action] += 1

    grad[x, :] .-= π.probs
    grad[x, action] += 1.

    return logp
end

function gradient_logp!(grad::Array{T}, π::LinearSoftmaxPolicy{T, IdentityBasis}, x::Array{T, 1}, action::Int)::T where {T<:Real}
    grad = reshape(grad, size(get_params(π.f)))
    fill!(grad, 0.)

    get_action_probabilities!(π, x)
    probs = -π.probs
    logp = log(π.probs[action])
    probs[action] += 1.

    for i in 1:length(probs)
        @. grad[:, i] = x * probs[i]
    end
    # grad .= x*probs'

    return logp
end

function gradient_logp!(grad::Array{T}, π::LinearSoftmaxPolicy{T}, x, action::Int)::T where {T<:Real}
    grad = reshape(grad, size(get_params(π.f)))
    fill!(grad, 0.)

    get_action_probabilities!(π, x)
    probs = -π.probs
    logp = log(π.probs[action])
    probs[action] += 1.
    for i in 1:length(probs)
        @. grad[:, i] = π.f.feats * probs[i]
    end
    # grad .= π.f.feats*probs'

    return logp
end

function gradient_entropy!(grad, π::LinearSoftmaxPolicy{T, IdentityBasis}, x::Int; reuse::Bool=false) where {T}
    grad = reshape(grad, size(get_params(π.f)))
    fill!(grad, 0.)
    if !reuse
        get_action_probabilities!(π, x)
    end

    p1 = zeros(length(π.probs))
    @. p1 = -π.probs * log(π.probs)
    H = sum(p1)  # entropy (p1 has negative in it)
    @. p1 += π.probs * H
    grad[x, :] .= p1
    return H
end

function gradient_entropy!(grad, π::LinearSoftmaxPolicy{T, IdentityBasis}, x::Array{T, 1}; reuse::Bool=false) where {T}
    grad = reshape(grad, size(get_params(π.f)))
    fill!(grad, 0.)
    if !reuse
        get_action_probabilities!(π, x)
    end

    p1 = zeros(length(π.probs))
    @. p1 = -π.probs * log(π.probs)
    H = sum(p1)  # entropy (p1 has negative in it)
    @. p1 += π.probs * H

    for i in 1:length(p1)
        @. grad[:, i] = x * p1[i]
    end
    # grad .= x*p1'
    return H
end

function gradient_entropy!(grad, π::LinearSoftmaxPolicy{T}, x; reuse::Bool=false) where {T}
    grad = reshape(grad, size(get_params(π.f)))
    fill!(grad, 0.)
    if !reuse
        get_action_probabilities!(π, x)
    end

    p1 = zeros(length(π.probs))
    @. p1 = -π.probs * log(π.probs)
    H = sum(p1)  # entropy (p1 has negative in it)
    @. p1 += π.probs * H

    for i in 1:length(p1)
        @. grad[:, i] = π.f.feats * p1[i]
    end

    # grad .= π.f.feats*p1'
    return H
end

function get_action_gradient_logp!(grad::Array{T}, π::LinearSoftmaxPolicy{T, IdentityBasis}, x::Int, rng::AbstractRNG)::T where {T<:Real}
    get_action_probabilities!(π, x)
    π.action = sample_discrete(π.probs, rng)
    logp = log(π.probs[π.action])
    # probs = -π.probs
    # probs[π.action] += 1.
    fill!(grad, 0.)
    grad[x, :] .-= π.probs
    grad[x, π.action] += 1.
    # grad[x, :] .= probs
    return logp
end

function get_action_gradient_logp!(grad::Array{T}, π::LinearSoftmaxPolicy{T, IdentityBasis}, x::Array{T, 1}, rng::AbstractRNG)::T where {T<:Real}
    get_action_probabilities!(π, x)
    π.action = sample_discrete(π.probs, rng)
    logp = log(π.probs[π.action])
    probs = -π.probs
    probs[π.action] += 1.
    grad .= x*probs'
    return logp
end

function get_action_gradient_logp!(grad::Array{T}, π::LinearSoftmaxPolicy{T}, x, rng::AbstractRNG)::T where {T<:Real}
    get_action_probabilities!(π, x)
    π.action = sample_discrete(π.probs, rng)
    logp = log(π.probs[π.action])
    probs = deepcopy(-π.probs)
    probs[π.action] += 1.
    grad .= π.f.feats*probs'
    return logp
end

function get_action_gradient_logp(grad::Array{T}, π::LinearSoftmaxPolicy{T}, x, rng::AbstractRNG) where {T<:Real}
    grad = zeros(T, size(get_params(π)))
    logp = get_action_gradient_logp!(grad, π, x)
    action = π.action
    return action, grad, logp
end

function set_params!(π::LinearSoftmaxPolicy{T}, θ) where {T<:Real}
    set_params!(π.f, θ)
end

function get_params(π::LinearSoftmaxPolicy)
    get_params(π.f)
end

function copy_params!(params::Array{T}, π::LinearSoftmaxPolicy{T}) where {T<:Real}
    copy_params!(params, π.f)
end

function copy_params(π::LinearSoftmaxPolicy{T})::Array{T} where {T<:Real}
    return copy_params(π.f)
end

function add_to_params!(π::LinearSoftmaxPolicy{T}, Δθ) where {T<:Real}
    add_to_params!(π.f, Δθ)
end

function clone(π::LinearSoftmaxPolicy)::LinearSoftmaxPolicy
    π₂ = LinearSoftmaxPolicy(π.f)
    π₂.action = deepcopy(π.action)
    π₂.probs = deepcopy(π.probs)
    return π₂
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
# include("constructor.jl")

# Tables
Profs1 = DataFrame(course = ["Artificial Intelligence", "Robotics"], professor = ["Prof. A", "Prof. B"] )

P3 = Person(age=24, field="Mechanical Engineering", courses = ["Artificial Intelligence", "Robotics"], professors= Profs1, id = 1)

pprint(P3)
#=
<person id="1">
<age>24</age>
<study-field>Mechanical Engineering</study-field>
<GPA>4.5</GPA>
<taken-courses>Artificial Intelligence</taken-courses>
<taken-courses>Robotics</taken-courses>

course
professor

String
String

Artificial Intelligence
Prof. A

Robotics
Prof. B

</person>
=#


text:
# start from the REPL of with:
# julia http_websockets.jl
# --> Client says hello!
using HTTP

@async HTTP.WebSockets.listen("127.0.0.1", UInt16(8081)) do ws  # 1
  while !eof(ws)
    data = readavailable(ws)    # 2
    write(ws, data)             # 3   - echo
  end
end
# Task (runnable) @0x000000000bdd2100

HTTP.WebSockets.open("ws://127.0.0.1:8081") do ws   # 4
  write(ws, "Client says hello!")  # 5
  x = readavailable(ws)            # 6
  println(String(x))               # 7
end

# Response:
# Client says hello!
# HTTP.Messages.Response:
# """
# HTTP/1.1 101 Switching Protocols
# Upgrade: websocket
# Connection: Upgrade
# Sec-WebSocket-Accept: yZuObElXGNUo/FJrtSlJpYuup7g=

# """


text:
module CategoricalTest

using Test
using ReactiveMP
using Distributions
using Random

@testset "Categorical" begin

    # Categorical comes from Distributions.jl and most of the things should be covered there
    # Here we test some extra ReactiveMP.jl specific functionality

    @testset "vague" begin
        @test_throws MethodError vague(Categorical)

        d1 = vague(Categorical, 2)

        @test typeof(d1) <: Categorical
        @test probvec(d1) ≈ [ 0.5, 0.5 ]

        d2 = vague(Categorical, 4)

        @test typeof(d2) <: Categorical
        @test probvec(d2) ≈ [ 0.25, 0.25, 0.25, 0.25 ]
    end

    @testset "prod" begin
        @test prod(ProdAnalytical(), Categorical([ 0.1, 0.4, 0.5 ]), Categorical([ 1/3, 1/3, 1/3 ])) == Categorical([ 0.1, 0.4, 0.5 ])
        @test prod(ProdAnalytical(), Categorical([ 0.1, 0.4, 0.5 ]), Categorical([ 0.8, 0.1, 0.1 ])) == Categorical([0.47058823529411764, 0.23529411764705882, 0.2941176470588235])
        @test prod(ProdAnalytical(), Categorical([ 0.2, 0.6, 0.2 ]), Categorical([ 0.8, 0.1, 0.1 ])) == Categorical([0.6666666666666666, 0.24999999999999994, 0.08333333333333333])
    end

    @testset "probvec" begin
        @test probvec(Categorical([ 0.1, 0.4, 0.5 ])) == [ 0.1, 0.4, 0.5 ]
        @test probvec(Categorical([ 1/3, 1/3, 1/3 ])) == [ 1/3, 1/3, 1/3 ]
        @test probvec(Categorical([ 0.8, 0.1, 0.1 ])) == [ 0.8, 0.1, 0.1 ]
    end

end

end


text:
using KhepriGrasshopper
using Test

@testset "KhepriGrasshopper.jl" begin
    # Write your tests here.
end


text:
## Standard stabilizers (non-exhaustive)
#==========================================================================================#

HorizontalTail(; root_chord = 1.0, taper = 1.0, span = 1.0, sweep = 0., w_sweep = 0., root_twist = 0., tip_twist = 0., root_foil = naca4(0,0,1,2), tip_foil = root_foil, angle = 0.) = WingSection(root_chord = root_chord, taper = taper, span = span, sweep = sweep, w_sweep = w_sweep, root_twist = root_twist, tip_twist = tip_twist, root_foil = naca4(0,0,1,2), tip_foil = tip_foil, angle = angle, axis = [0., 1., 0.])

VerticalTail(; root_chord = 1.0, taper = 1.0, span = 1.0, sweep = 0., w_sweep = 0., root_foil = naca4(0,0,1,2), tip_foil = root_foil) = HalfWingSection(root_chord = root_chord, taper = taper, span = span, sweep = sweep, w_sweep = w_sweep, tip_foil = tip_foil, angle = 90., axis = [1., 0., 0.])

## Prototype for generic control deflection on panel mesh
#==========================================================================================#

abstract type AbstractControlSurface end

struct BoundingBox{T <: Real} <: AbstractControlSurface
    location    :: SVector{3,T}
    span_ratio  :: T
    chord_ratio :: T
    function BoundingBox{T}(x_o, s_r, c_r) where T <: Real
        @assert all(0. < x < 1., x_o) "Origin of bounding box must lie within [0,1] ⊂ ℝ!"
        @assert 0. <= s_r <= 1. "The normalized span length must lie within [0,1] ⊂ ℝ!"
        @assert 0. <= c_r <= 1. "The normalized chord length must lie within [0,1] ⊂ ℝ!"

        new{T}(x_o, s_r, c_r)
    end
end

function Aileron(bb :: BoundingBox, δ)

end