text:
using Pkg
if isfile("Project.toml") && isfile("Manifest.toml")
    Pkg.activate(".")
end

###################
# Terminal ########
###################

# Read from terminal user input
case =  ARGS[1]		# First argument sets case
timeperiod = ARGS[2]	# Second argument sets time period
trm = parse(Float64, ARGS[3])
η_E = parse(Float64, ARGS[4])
η_M = parse(Float64, ARGS[5])

timeperiod_str1 = replace(string(timeperiod), ":" => "-")
timeperiod_str = timeperiod_str1*"__"*ARGS[4]*"__"*ARGS[5]

println()
println("Case >" * case * "< selected.")
println("Calculating for time period >" * timeperiod * "<.")
println()

###################
# Run CM ##########
###################

println("Loading required packages: ProgressMeter, CSV, DataFrames, Dates, JuMP, Gurobi.")
using ProgressMeter
using CSV
using DataFrames
using Dates
using JuMP
using Gurobi
println("Loading of packages complete.")

include("model/src/Sesam_ptdf.jl")
using .Sesam

###################
# Data ############
###################

println("Reading data from csv...")
prog = Progress(14)

# Read nodes and load from .csv
	nodes_df = CSV.read(joinpath("data", case, "nodes.csv"))
		next!(prog)
	load_df = CSV.read(joinpath("data", case, "load.csv"))
		next!(prog)
	exchange_df = CSV.read(joinpath("data", case, "exchange.csv"))
		next!(prog)
# Read lines from .csv
	lines_df = CSV.read(joinpath("data", case, "lines.csv"))
		next!(prog)
# Read power plants from .csv
	powerplants_df = CSV.read(joinpath("data", case, "powerplants.csv"))
		next!(prog)
	avail_powerplants_df = CSV.read(joinpath("data", case, "avail_powerplants.csv"), header = 2, skipto = 3)
		next!(prog)
	mustrun_powerplants_df = CSV.read(joinpath("data", case, "mustrun_powerplants.csv"), header = 2, skipto = 3)
		next!(prog)
# Read fuel costs from .csv
	fuelcost_df = CSV.read(joinpath("data", case, "fuelcost.csv"))
		next!(prog)
# Read renewables and availabilities from .csv
	renewables_df = CSV.read(joinpath("data", case, "renewables.csv"))
		next!(prog)
	avail_solar_pv_df = CSV.read(joinpath("data", case, "avail_solar_pv.csv"))
		next!(prog)
	avail_wind_offshore_df = CSV.read(joinpath("data", case, "avail_wind_offshore.csv"))
		next!(prog)
	avail_wind_onshore_df = CSV.read(joinpath("data", case, "avail_wind_onshore.csv"))
		next!(prog)
	avail_other_res_df = CSV.read(joinpath("data", case, "avail_other_res.csv"))
		next!(prog)
	renewables_avail_dict = Dict(:SolarPV => avail_solar_pv_df,
		:WindOffshore => avail_wind_offshore_df,
		:WindOnshore => avail_wind_onshore_df,
		:other => avail_other_res_df)
		next!(prog)

println("Reading of data complete.")
println()

###################
# Model setup #####
###################

println("Preparing model technologies, nodes, and grid...")
prog = Progress(6)

# Create set of nodes and lines
	nodes = Nodes(nodes_df, load_df, exchange_df, 1000, 1)
		next!(prog)
	lines = Lines(lines_df, 1000, trm)
		next!(prog)
# Create set of power plants	ED_L_S_opt, ED_p
	pp = PowerPlants(powerplants_df, avail_powerplants_df, mustrun_powerplants_df, fuelcost_df)
		next!(prog)
# Create set of renewables
	res = Renewables(renewables_df, renewables_avail_dict)
		next!(prog)
# Create a set of storages
	storages_df = CSV.read(joinpath("data", case, "storages.csv"))
		next!(prog)
	storages = Storages(storages_df)
		next!(prog)

println("Preparation complete.")
println()

session_time = Dates.now()
session = Dates.format(session_time, "yyyy-mm-dd__HH-MM-SS")*"__ptdf"

# Create session folder for output
if isdir("output")
	if !isdir(joinpath("output", session*"__"*case*"__"*timeperiod_str))
		mkdir(joinpath("output", session*"__"*case*"__"*timeperiod_str))
	end
else
	mkdir("output")
	mkdir(joinpath("output", session*"__"*case*"__"*timeperiod_str))
end

prog = Progress(length(dayslicer(timeperiod)))

### Parameter outputs
println("Exporting parameters.")
println("Exporting merit order...")
# Merit order
output_ED_merit_order_pp = DataFrame(unit=Symbol[],fuel=String[],capacity=Float64[],mc=Float64[])

try
	for g in pp.unit
			mc = pp.mc[g]
			fuel = pp.fuel[g]
			capacity = pp.capacity[g]
			push!(output_ED_merit_order_pp, [g fuel capacity mc])
	end

	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_merit_order_pp.csv"), output_ED_merit_order_pp)
catch e
	println("The marginal costs for power plants are not available.")
end

output_ED_merit_order_res = DataFrame(unit=Symbol[],fuel=Symbol[],capacity=Float64[],mc=Float64[])

try
	for r in res.unit
			mc = 0
			fuel = res.fuel[r]
			capacity = res.capacity[r]
			push!(output_ED_merit_order_res, [r fuel capacity mc])
	end

	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_merit_order_res.csv"), output_ED_merit_order_res)
catch e
	println("The marginal costs for renewables are not available.")
end

output_ED_merit_order_storages = DataFrame(unit=Symbol[],fuel=String[],capacity=Float64[],mc=Float64[])

try
	for s in storages.unit
			mc = 0
			fuel = storages.fuel[s]
			capacity = storages.power[s]
			push!(output_ED_merit_order_storages, [s fuel capacity mc])
	end

	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_merit_order_storages.csv"), output_ED_merit_order_storages)
catch e
	println("The marginal costs for renewables are not available.")
end

println("Export complete.")
println()

### Initial run of first dayslice
println("Initial run. Calculating first dayslice: " * string(dayslicer(timeperiod)[1]))
println()
println()

ED_mod, ED_P_opt, ED_P_R_opt, ED_P_R_opt_dist, ED_P_S_opt, ED_D_S_opt,
	ED_L_S_opt, ED_price = EconomicDispatch80(dayslicer(timeperiod)[1], nodes, pp, res, storages)

CM_PtG_mod, CM_PtG_ΔP_up_opt, CM_PtG_ΔP_dn_opt, CM_PtG_ΔP_R_up_opt, CM_PtG_ΔP_R_dn_opt,
	CM_PtG_Θ_opt, CM_PtG_P_flow_opt, CM_PtG_price, CM_PtG_P_gen_lost_opt, CM_PtG_P_load_lost_opt,
	CM_PtG_P_S_up_opt, CM_PtG_P_syn_opt, CM_PtG_D_PtG_opt, CM_PtG_L_syn_opt = CongestionManagementPHS80PtG(dayslicer(timeperiod)[1],
	nodes, lines, pp, res, storages, ED_price, ED_P_opt, ED_P_R_opt_dist,
	ED_P_S_opt, ED_D_S_opt, ED_L_S_opt, η_E, η_M)

next!(prog)

	# Dispatchable power plants
	output_ED_P_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])

	try
		for t in dayslicer(timeperiod)[1]
			for g in pp.unit
				val = ED_P_opt[t, g]
				fuel =  pp.fuel[g]
				push!(output_ED_P_opt, [t g fuel val pp.node[g]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_P.csv"), output_ED_P_opt)
	catch e
		println("ED_P_opt does not exist.")
	end

	# RES generation
	output_ED_P_R_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for r in res.unit
				val = ED_P_R_opt[t, r]
				fuel =  res.fuel[r]
				push!(output_ED_P_R_opt, [t r fuel val res.node[r]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_P_R.csv"), output_ED_P_R_opt)
	catch e
		println("ED_P_R_opt does not exist.")
	end

	# Maximum available RES infeed from input data
	output_MaxRES = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for r in res.unit
				output = res.infeed[r][t]
				fuel =  res.fuel[r]
				push!(output_MaxRES, [t r fuel output res.node[r]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "MaxRES.csv"), output_MaxRES)
	catch e
		println("res.infeed does not exist.")
	end

	# Generation from storages
	output_ED_P_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for s in storages.unit
				val = ED_P_S_opt[t, s]
				fuel =  storages.fuel[s]
				push!(output_ED_P_S_opt, [t s fuel val storages.node[s]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_P_S.csv"), output_ED_P_S_opt)
	catch e
		println("ED_P_S_opt does not exist.")
	end

	# Electricity demand from storages
	output_ED_D_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for s in storages.unit
				val = -1*ED_D_S_opt[t, s]
				fuel =  storages.fuel[s]
				push!(output_ED_D_S_opt, [t s fuel val storages.node[s]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_D_S.csv"), output_ED_D_S_opt)
	catch e
		println("ED_D_S_opt does not exist.")
	end

	# Storage level
	output_ED_L_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],storage=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for s in storages.unit
				val = ED_L_S_opt[t, s]
				fuel =  storages.fuel[s]
				push!(output_ED_L_S_opt, [t s fuel val storages.node[s]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_L_S.csv"), output_ED_L_S_opt)
	catch e
		println("ED_L_S_opt does not exist.")
	end

	# Load
	output_load = DataFrame(time=Int64[],load=Float64[])
	try
		for t in dayslicer(timeperiod)[1]
			load = nodes.systemload[t]
			push!(output_load, [t load])
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "load.csv"), output_load)
	catch e
		println("nodes.systemload does not exist.")
	end

	# Prices
	output_ED_price = DataFrame(time=Int64[],price=Float64[])
	try
		for t in dayslicer(timeperiod)[1]
			price =  ED_price[t]
			push!(output_ED_price, [t price])
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_price.csv"), output_ED_price)
	catch e
		println("ED_price does not exist.")
	end

	# Must-run
	output_pmin = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for g in pp.unit
				val = pp.pmin[g][t]
				fuel =  pp.fuel[g]
				push!(output_pmin, [t g fuel val pp.node[g]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "pmin.csv"), output_pmin)
	catch e
		println("pp.pmin does not exist.")
	end

	### Model parameters
	output_ED_info = DataFrame(timeperiod=String[],
									   termination_status=String[],
									   primal_status=String[],
									   dual_status=String[],
									   objective_value=Float64[],
									   solve_time=Float64[])

	push!(output_ED_info, [string(dayslicer(timeperiod)[1]) string(termination_status(ED_mod)) string(primal_status(ED_mod)) string(dual_status(ED_mod)) objective_value(ED_mod) solve_time(ED_mod)])
	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_info.csv"), output_ED_info)

	# Nodal exchange
	output_exchange = DataFrame(time=Int64[],node=Int64[],output=Float64[])
	try
		for t in dayslicer(timeperiod)[1]
			for j in nodes.id
				push!(output_exchange, [t j nodes.exchange[j][t]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "exchange.csv"), output_exchange)
	catch e
		println("output_exchange does not exist.")
	end

	###################
	# CM ##############
	###################

	# Dispatchable power plants
	output_CM_PtG_ΔP_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for g in pp.unit
				val = CM_PtG_ΔP_up_opt[t, g]
				fuel =  pp.fuel[g]
				push!(output_CM_PtG_ΔP_up_opt, [t g fuel val pp.node[g]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_up.csv"), output_CM_PtG_ΔP_up_opt)
	catch e
		println("CM_PtG_ΔP_up_opt does not exist.")
	end

	output_CM_PtG_ΔP_dn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for g in pp.unit
				val = -1*CM_PtG_ΔP_dn_opt[t, g]
				fuel =  pp.fuel[g]
				push!(output_CM_PtG_ΔP_dn_opt, [t g fuel val pp.node[g]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_dn.csv"), output_CM_PtG_ΔP_dn_opt)
	catch e
		println("CM_PtG_ΔP_dn_opt does not exist.")
	end

	# RES generation
	output_CM_PtG_ΔP_R_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for r in res.unit
				val = CM_PtG_ΔP_R_up_opt[t, r]
				fuel =  res.fuel[r]
				push!(output_CM_PtG_ΔP_R_up_opt, [t r fuel val res.node[r]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_R_up.csv"), output_CM_PtG_ΔP_R_up_opt)
	catch e
		println("CM_PtG_ΔP_R_up_opt does not exist.")
	end

	output_CM_PtG_ΔP_R_dn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for r in res.unit
				val = -1*CM_PtG_ΔP_R_dn_opt[t, r]
				fuel =  res.fuel[r]
				push!(output_CM_PtG_ΔP_R_dn_opt, [t r fuel val res.node[r]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_R_dn.csv"), output_CM_PtG_ΔP_R_dn_opt)
	catch e
		println("CM_PtG_ΔP_R_dn_opt does not exist.")
	end

	### Model parameters
	output_CM_PtG_info = DataFrame(timeperiod=String[],
									   termination_status=String[],
									   primal_status=String[],
									   dual_status=String[],
									   objective_value=Float64[],
									   solve_time=Float64[])

	push!(output_CM_PtG_info, [string(dayslicer(timeperiod)[1]) string(termination_status(CM_PtG_mod)) string(primal_status(CM_PtG_mod)) string(dual_status(CM_PtG_mod)) objective_value(CM_PtG_mod) solve_time(CM_PtG_mod)])
	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_info.csv"), output_CM_PtG_info)

	### Power flow results
	output_CM_PtG_P_flow_opt = DataFrame(time=Int64[],
									 line=Int64[],
									 from=Int64[],
									 to=Int64[],
									 pmax=Float64[],
									 pflow=Float64[])

	 try
 		for t in dayslicer(timeperiod)[1]
 			for l in lines.id
				from = parse(Int, string(lines.from[l]))
				to = parse(Int, string(lines.to[l]))
				pmax = lines.pmax[l]
				pflow = CM_PtG_P_flow_opt[t, l]
 				push!(output_CM_PtG_P_flow_opt, [t l from to pmax pflow])
 			end
 		end
 		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_flow.csv"), output_CM_PtG_P_flow_opt)
 	catch e
 		println("CM_PtG_P_flow_opt does not exist.")
 	end

	# Lost generation
	output_CM_PtG_P_gen_lost = DataFrame(time=Int64[],node=Int64[],output=Float64[])
	try
		for t in dayslicer(timeperiod)[1]
			for j in nodes.id
				output = CM_PtG_P_gen_lost_opt[t, j]
				push!(output_CM_PtG_P_gen_lost, [t, j, output])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_gen_lost.csv"), output_CM_PtG_P_gen_lost)
	catch e
		println("CM_PtG_P_gen_lost does not exist.")
	end

	# Lost load
	output_CM_PtG_P_load_lost = DataFrame(time=Int64[],node=Int64[],output=Float64[])
	try
		for t in dayslicer(timeperiod)[1]
			for j in nodes.id
				output = CM_PtG_P_load_lost_opt[t, j]
				push!(output_CM_PtG_P_load_lost, [t, j, output])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_load_lost.csv"), output_CM_PtG_P_load_lost)
	catch e
		println("CM_PtG_P_load_lost does not exist.")
	end

	output_CM_PtG_price = DataFrame(time=Int64[],node=Int64[],price=Float64[])
	 try
		 for t in dayslicer(timeperiod)[1]
			 for j in nodes.id
				 price =  CM_PtG_price[j, t]
				 push!(output_CM_PtG_price, [t j price])
			 end
		 end
		 CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_price.csv"), output_CM_PtG_price)
	 catch e
		 println("CM_PtG_price does not exist.")
	 end

	# PHS output increase
	output_CM_PtG_P_S_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for s in storages.unit
				val = CM_PtG_P_S_up_opt[t, s]
				fuel =  storages.fuel[s]
				push!(output_CM_PtG_P_S_up_opt, [t s fuel val storages.node[s]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_S_up.csv"), output_CM_PtG_P_S_up_opt)
	catch e
		println("CM_PtG_P_S_up_opt does not exist.")
	end

	# Electricity generation from synthetic methane
	output_CM_PtG_P_syn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for e in getKeyVector(pp.fuel, "NaturalGas")
				val = CM_PtG_P_syn_opt[t, e]
				fuel = "SNG"
				push!(output_CM_PtG_P_syn_opt, [t e fuel val pp.node[e]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_syn.csv"), output_CM_PtG_P_syn_opt)
	catch e
		println("CM_PtG_P_syn_opt does not exist.")
	end

	# Electricity demand from PtG units
	output_CM_PtG_D_PtG_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	try
		for t in dayslicer(timeperiod)[1]
			for r in CM_PtG_D_PtG_opt.axes[2]
				val = CM_PtG_D_PtG_opt[t, r]
				fuel = "PtG"
				push!(output_CM_PtG_D_PtG_opt, [t r fuel val res.node[r]])
			end
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_D_PtG.csv"), output_CM_PtG_D_PtG_opt)
	catch e
		println("CM_PtG_D_PtG_opt does not exist.")
	end

	# Synthetic methane storage level
	output_CM_PtG_L_syn_opt = DataFrame(time=Int64[],storage=Float64[])
	try
		for t in dayslicer(timeperiod)[1]
			storage = CM_PtG_L_syn_opt[t]
			push!(output_CM_PtG_L_syn_opt, [t storage])
		end
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_L_syn.csv"), output_CM_PtG_L_syn_opt)
	catch e
		println("CM_PtG_L_syn_opt does not exist.")
	end

### Remaining periods: Initialising first hour with last hour of previous day
# 	and append data

for hours in dayslicer(timeperiod)[2:end]
	hours_extended = [hours[1]-1; hours]

	println()
	println()
	println("Calculating economic dispatch: " * string(hours))

	global ED_mod, ED_P_opt, ED_P_R_opt, ED_P_R_opt_dist, ED_P_S_opt, ED_D_S_opt,
		ED_L_S_opt, ED_price = EconomicDispatch80(hours, nodes, pp, res, storages)

	    ###################
	    # ED ##############
	    ###################

	    # Dispatchable power plants
	    output_ED_P_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])

	    try
	    	for t in hours
	    		for g in pp.unit
	    			val = ED_P_opt[t, g]
	    			fuel =  pp.fuel[g]
	    			push!(output_ED_P_opt, [t g fuel val pp.node[g]])
	    		end
	    	end
	    	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_P.csv"), output_ED_P_opt, append = true)
	    catch e
	    	println("ED_P_opt does not exist.")
	    end

	    # RES generation
	    output_ED_P_R_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for r in res.unit
	    			val = ED_P_R_opt[t, r]
	    			fuel =  res.fuel[r]
	    			push!(output_ED_P_R_opt, [t r fuel val res.node[r]])
	    		end
	    	end
	    	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_P_R.csv"), output_ED_P_R_opt, append = true)
	    catch e
	    	println("ED_P_R_opt does not exist.")
	    end

	    # Maximum available RES infeed from input data
	    output_MaxRES = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for r in res.unit
	    			output = res.infeed[r][t]
	    			fuel =  res.fuel[r]
	    			push!(output_MaxRES, [t r fuel output res.node[r]])
	    		end
	    	end
	    	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "MaxRES.csv"), output_MaxRES, append = true)
	    catch e
	    	println("res.infeed does not exist.")
	    end

	    # Generation from storages
	    output_ED_P_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for s in storages.unit
	    			val = ED_P_S_opt[t, s]
	    			fuel =  storages.fuel[s]
	    			push!(output_ED_P_S_opt, [t s fuel val storages.node[s]])
	    		end
	    	end
	    	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_P_S.csv"), output_ED_P_S_opt, append = true)
	    catch e
	    	println("ED_P_S_opt does not exist.")
	    end

	    # Electricity demand from storages
	    output_ED_D_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for s in storages.unit
	    			val = -1*ED_D_S_opt[t, s]
	    			fuel =  storages.fuel[s]
	    			push!(output_ED_D_S_opt, [t s fuel val storages.node[s]])
	    		end
	    	end
	    	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_D_S.csv"), output_ED_D_S_opt, append = true)
	    catch e
	    	println("ED_D_S_opt does not exist.")
	    end

	    # Storage level
	    output_ED_L_S_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],storage=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for s in storages.unit
	    			val = ED_L_S_opt[t, s]
	    			fuel =  storages.fuel[s]
	    			push!(output_ED_L_S_opt, [t s fuel val storages.node[s]])
	    		end
	    	end
	    	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_L_S.csv"), output_ED_L_S_opt, append = true)
	    catch e
	    	println("ED_L_S_opt does not exist.")
	    end

	    # Load
	    output_load = DataFrame(time=Int64[],load=Float64[])
	    try
	    	for t in hours
	    		load = nodes.systemload[t]
	    		push!(output_load, [t load])
	    	end
	    	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "load.csv"), output_load, append = true)
	    catch e
	    	println("nodes.systemload does not exist.")
	    end

	    # Prices
	    output_ED_price = DataFrame(time=Int64[],price=Float64[])
	    try
	    	for t in hours
	    		price =  ED_price[t]
	    		push!(output_ED_price, [t price])
	    	end
	    	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_price.csv"), output_ED_price, append = true)
	    catch e
	    	println("ED_price does not exist.")
	    end

	    # Must-run
	    output_pmin = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
	    try
	    	for t in hours
	    		for g in pp.unit
	    			val = pp.pmin[g][t]
	    			fuel =  pp.fuel[g]
	    			push!(output_pmin, [t g fuel val pp.node[g]])
	    		end
	    	end
	    	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "pmin.csv"), output_pmin, append = true)
	    catch e
	    	println("pp.pmin does not exist.")
	    end

		### Model parameters
		output_ED_info = DataFrame(timeperiod=String[],
										   termination_status=String[],
										   primal_status=String[],
										   dual_status=String[],
										   objective_value=Float64[],
										   solve_time=Float64[])

		push!(output_ED_info, [string(hours) string(termination_status(ED_mod)) string(primal_status(ED_mod)) string(dual_status(ED_mod)) objective_value(ED_mod) solve_time(ED_mod)])
		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "ED_info.csv"), output_ED_info, append = true)

		###################
		# CM ##############
		###################

		println()
		println()
		println("Calculating congestion management: " * string(hours))

			global CM_PtG_mod, CM_PtG_ΔP_up_opt, CM_PtG_ΔP_dn_opt, CM_PtG_ΔP_R_up_opt, CM_PtG_ΔP_R_dn_opt,
				CM_PtG_Θ_opt, CM_PtG_P_flow_opt, CM_PtG_price, CM_PtG_P_gen_lost_opt, CM_PtG_P_load_lost_opt,
				CM_PtG_P_S_up_opt, CM_PtG_P_syn_opt, CM_PtG_D_PtG_opt, CM_PtG_L_syn_opt = CongestionManagementPHS80PtG(hours, nodes, lines, pp, res, storages, ED_price,
				ED_P_opt,  ED_P_R_opt_dist, ED_P_S_opt, ED_D_S_opt, ED_L_S_opt, η_E, η_M)

			# Dispatchable power plants
			output_CM_PtG_ΔP_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
			try
				for t in hours
					for g in pp.unit
						val = CM_PtG_ΔP_up_opt[t, g]
						fuel =  pp.fuel[g]
						push!(output_CM_PtG_ΔP_up_opt, [t g fuel val pp.node[g]])
					end
				end
				CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_up.csv"), output_CM_PtG_ΔP_up_opt, append = true)
			catch e
				println("CM_PtG_ΔP_up_opt does not exist.")
			end

			output_CM_PtG_ΔP_dn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
			try
				for t in hours
					for g in pp.unit
						val = -1*CM_PtG_ΔP_dn_opt[t, g]
						fuel =  pp.fuel[g]
						push!(output_CM_PtG_ΔP_dn_opt, [t g fuel val pp.node[g]])
					end
				end
				CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_dn.csv"), output_CM_PtG_ΔP_dn_opt, append = true)
			catch e
				println("CM_PtG_ΔP_dn_opt does not exist.")
			end

			# RES generation
			output_CM_PtG_ΔP_R_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
			try
				for t in hours
					for r in res.unit
						val = CM_PtG_ΔP_R_up_opt[t, r]
						fuel =  res.fuel[r]
						push!(output_CM_PtG_ΔP_R_up_opt, [t r fuel val res.node[r]])
					end
				end
				CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_R_up.csv"), output_CM_PtG_ΔP_R_up_opt, append = true)
			catch e
				println("CM_PtG_ΔP_R_up_opt does not exist.")
			end

			output_CM_PtG_ΔP_R_dn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=Symbol[],output=Float64[],node=Int64[])
			try
				for t in hours
					for r in res.unit
						val = -1*CM_PtG_ΔP_R_dn_opt[t, r]
						fuel =  res.fuel[r]
						push!(output_CM_PtG_ΔP_R_dn_opt, [t r fuel val res.node[r]])
					end
				end
				CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_R_dn.csv"), output_CM_PtG_ΔP_R_dn_opt, append = true)
			catch e
				println("CM_PtG_ΔP_R_dn_opt does not exist.")
			end

			### Model parameters
			output_CM_PtG_info = DataFrame(timeperiod=String[],
											   termination_status=String[],
											   primal_status=String[],
											   dual_status=String[],
											   objective_value=Float64[],
											   solve_time=Float64[])

			push!(output_CM_PtG_info, [string(hours) string(termination_status(CM_PtG_mod)) string(primal_status(CM_PtG_mod)) string(dual_status(CM_PtG_mod)) objective_value(CM_PtG_mod) solve_time(CM_PtG_mod)])
			CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_info.csv"), output_CM_PtG_info, append = true)

			### Power flow results
			output_CM_PtG_P_flow_opt = DataFrame(time=Int64[],
											 line=Int64[],
											 from=Int64[],
											 to=Int64[],
											 pmax=Float64[],
											 pflow=Float64[])

			try
			   for t in hours
				   for l in lines.id
					   from = parse(Int, string(lines.from[l]))
					   to = parse(Int, string(lines.to[l]))
					   pmax = lines.pmax[l]
					   pflow = CM_PtG_P_flow_opt[t, l]
					   push!(output_CM_PtG_P_flow_opt, [t l from to pmax pflow])
				   end
			   end
			   CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_flow.csv"), output_CM_PtG_P_flow_opt, append = true)
		   catch e
			   println("CM_PtG_P_flow_opt does not exist.")
		   end

		   # Lost generation
		   output_CM_PtG_P_gen_lost = DataFrame(time=Int64[],node=Int64[],output=Float64[])
		   try
			   for t in hours
				   for j in nodes.id
					   output = CM_PtG_P_gen_lost_opt[t, j]
					   push!(output_CM_PtG_P_gen_lost, [t, j, output])
				   end
			   end
			   CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_gen_lost.csv"), output_CM_PtG_P_gen_lost, append = true)
		   catch e
			   println("CM_PtG_P_gen_lost does not exist.")
		   end

		   # Lost load
		   output_CM_PtG_P_load_lost = DataFrame(time=Int64[],node=Int64[],output=Float64[])
		   try
			   for t in hours
				   for j in nodes.id
					   output = CM_PtG_P_load_lost_opt[t, j]
					   push!(output_CM_PtG_P_load_lost, [t, j, output])
				   end
			   end
			   CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_load_lost.csv"), output_CM_PtG_P_load_lost, append = true)
		   catch e
			   println("CM_PtG_P_load_lost does not exist.")
		   end

		   # Noral Price
		   output_CM_PtG_price = DataFrame(time=Int64[],node=Int64[],price=Float64[])
	   	    try
	   	    	for t in hours
					for j in nodes.id
		   	    		price =  CM_PtG_price[j, t]
		   	    		push!(output_CM_PtG_price, [t j price])
					end
	   	    	end
	   	    	CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_price.csv"), output_CM_PtG_price, append = true)
	   	    catch e
	   	    	println("CM_PtG_price does not exist.")
	   	    end

			output_CM_PtG_P_S_up_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
			try
				for t in hours
					for s in storages.unit
						val = CM_PtG_P_S_up_opt[t, s]
						fuel =  storages.fuel[s]
						push!(output_CM_PtG_P_S_up_opt, [t s fuel val storages.node[s]])
					end
				end
				CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_S_up.csv"), output_CM_PtG_P_S_up_opt, append = true)
			catch e
				println("CM_PtG_P_S_up_opt does not exist.")
			end

			# Electricity generation from synthetic methane
		   	output_CM_PtG_P_syn_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
		   	try
		   		for t in hours
		   			for e in getKeyVector(pp.fuel, "NaturalGas")
		   				val = CM_PtG_P_syn_opt[t, e]
		   				fuel = "SNG"
		   				push!(output_CM_PtG_P_syn_opt, [t e fuel val pp.node[e]])
		   			end
		   		end
		   		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_P_syn.csv"), output_CM_PtG_P_syn_opt, append = true)
		   	catch e
		   		println("CM_PtG_P_syn_opt does not exist.")
		   	end

			# Electricity demand from PtG units
		   	output_CM_PtG_D_PtG_opt = DataFrame(time=Int64[],unit=Symbol[],fuel=String[],output=Float64[],node=Int64[])
		   	try
		   		for t in hours
		   			for r in CM_PtG_D_PtG_opt.axes[2]
		   				val = CM_PtG_D_PtG_opt[t, r]
		   				fuel = "PtG"
		   				push!(output_CM_PtG_D_PtG_opt, [t r fuel val res.node[r]])
		   			end
		   		end
		   		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_D_PtG.csv"), output_CM_PtG_D_PtG_opt, append = true)
		   	catch e
		   		println("CM_PtG_D_PtG_opt does not exist.")
		   	end

			# Synthetic methane storage level
		   	output_CM_PtG_L_syn_opt = DataFrame(time=Int64[],storage=Float64[])
		   	try
		   		for t in hours
		   			storage = CM_PtG_L_syn_opt[t]
		   			push!(output_CM_PtG_L_syn_opt, [t storage])
		   		end
		   		CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "CM_PtG_L_syn.csv"), output_CM_PtG_L_syn_opt, append = true)
		   	catch e
		   		println("CM_PtG_L_syn_opt does not exist.")
		   	end

			# Nodal exchange
			output_exchange = DataFrame(time=Int64[],node=Int64[],output=Float64[])
			try
				for t in hours
					for j in nodes.id
						push!(output_exchange, [t j nodes.exchange[j][t]])
					end
				end
				CSV.write(joinpath("output", session*"__"*case*"__"*timeperiod_str, "exchange.csv"), output_exchange, append = true)
			catch e
				println("output_exchange does not exist.")
			end

	next!(prog)

end

println()
println("Calculation complete.")

max_stars_repo_name: bobbyxiong/redispatch-ptg

max_issues_repo_name: bobbyxiong/redispatch-ptg

max_forks_repo_name: bobbyxiong/redispatch-ptg


text:
__precompile__()

module FirstOrderOptimization

# shared among all algorithms

include("stepsize.jl")
include("optparams.jl")
include("utilities.jl")
include("convergence.jl")

# functions

abstract AbstractData
include("MNL.jl")

# algorithms

include("algorithms/gradient_descent.jl")
include("algorithms/prox_grad.jl")
include("algorithms/PRISMA.jl")
include("algorithms/frank_wolfe.jl")

# fancy data structures

include("sketch/sketch.jl")
include("sketch/operator.jl")

# algorithms with fancy data structures

include("algorithms/frank_wolfe_sketched.jl")
include("algorithms/frank_wolfe_thin.jl")
# etc

end # module

ext: jl

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
module Proto

using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

struct __enum_Version <: ProtoEnum
    _START_VERSION::Int32
    IR_VERSION_2017_10_10::Int32
    IR_VERSION_2017_10_30::Int32
    IR_VERSION::Int32
    __enum_Version() = new(0,1,2,3)
end #struct __enum_Version
const Version = __enum_Version()

mutable struct StringStringEntryProto <: ProtoType
    key::AbstractString
    value::AbstractString
    StringStringEntryProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct StringStringEntryProto
hash(v::StringStringEntryProto) = ProtoBuf.protohash(v)
isequal(v1::StringStringEntryProto, v2::StringStringEntryProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::StringStringEntryProto, v2::StringStringEntryProto) = ProtoBuf.protoeq(v1, v2)

struct __enum_TensorProto_DataType <: ProtoEnum
    UNDEFINED::Int32
    FLOAT::Int32
    UINT8::Int32
    INT8::Int32
    UINT16::Int32
    INT16::Int32
    INT32::Int32
    INT64::Int32
    STRING::Int32
    BOOL::Int32
    FLOAT16::Int32
    DOUBLE::Int32
    UINT32::Int32
    UINT64::Int32
    COMPLEX64::Int32
    COMPLEX128::Int32
    __enum_TensorProto_DataType() = new(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)
end #struct __enum_TensorProto_DataType
const TensorProto_DataType = __enum_TensorProto_DataType()

mutable struct TensorProto_Segment <: ProtoType
    _begin::Int64
    _end::Int64
    TensorProto_Segment(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TensorProto_Segment
hash(v::TensorProto_Segment) = ProtoBuf.protohash(v)
isequal(v1::TensorProto_Segment, v2::TensorProto_Segment) = ProtoBuf.protoisequal(v1, v2)
==(v1::TensorProto_Segment, v2::TensorProto_Segment) = ProtoBuf.protoeq(v1, v2)

mutable struct TensorProto <: ProtoType
    dims::Vector{Int64}
    data_type::Int32
    segment::TensorProto_Segment
    float_data::Vector{Float32}
    int32_data::Vector{Int32}
    string_data::Vector{Array{UInt8,1}}
    int64_data::Vector{Int64}
    name::AbstractString
    doc_string::AbstractString
    raw_data::Array{UInt8,1}
    double_data::Vector{Float64}
    uint64_data::Vector{UInt64}
    TensorProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TensorProto
const __fnum_TensorProto = Int[1,2,3,4,5,6,7,8,12,9,10,11]
const __pack_TensorProto = Symbol[:float_data,:int32_data,:int64_data,:double_data,:uint64_data]
meta(t::Type{TensorProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_TensorProto, ProtoBuf.DEF_VAL, true, __pack_TensorProto, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)
hash(v::TensorProto) = ProtoBuf.protohash(v)
isequal(v1::TensorProto, v2::TensorProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::TensorProto, v2::TensorProto) = ProtoBuf.protoeq(v1, v2)

mutable struct TensorShapeProto_Dimension <: ProtoType
    dim_value::Int64
    dim_param::AbstractString
    TensorShapeProto_Dimension(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TensorShapeProto_Dimension
const __oneofs_TensorShapeProto_Dimension = Int[1,1]
const __oneof_names_TensorShapeProto_Dimension = [Symbol("value")]
meta(t::Type{TensorShapeProto_Dimension}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, __oneofs_TensorShapeProto_Dimension, __oneof_names_TensorShapeProto_Dimension, ProtoBuf.DEF_FIELD_TYPES)
hash(v::TensorShapeProto_Dimension) = ProtoBuf.protohash(v)
isequal(v1::TensorShapeProto_Dimension, v2::TensorShapeProto_Dimension) = ProtoBuf.protoisequal(v1, v2)
==(v1::TensorShapeProto_Dimension, v2::TensorShapeProto_Dimension) = ProtoBuf.protoeq(v1, v2)

mutable struct TensorShapeProto <: ProtoType
    dim::Vector{TensorShapeProto_Dimension}
    TensorShapeProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TensorShapeProto
hash(v::TensorShapeProto) = ProtoBuf.protohash(v)
isequal(v1::TensorShapeProto, v2::TensorShapeProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::TensorShapeProto, v2::TensorShapeProto) = ProtoBuf.protoeq(v1, v2)

mutable struct TypeProto_Tensor <: ProtoType
    elem_type::Int32
    shape::TensorShapeProto
    TypeProto_Tensor(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TypeProto_Tensor
hash(v::TypeProto_Tensor) = ProtoBuf.protohash(v)
isequal(v1::TypeProto_Tensor, v2::TypeProto_Tensor) = ProtoBuf.protoisequal(v1, v2)
==(v1::TypeProto_Tensor, v2::TypeProto_Tensor) = ProtoBuf.protoeq(v1, v2)

mutable struct TypeProto <: ProtoType
    tensor_type::TypeProto_Tensor
    TypeProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct TypeProto
const __oneofs_TypeProto = Int[1]
const __oneof_names_TypeProto = [Symbol("value")]
meta(t::Type{TypeProto}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, __oneofs_TypeProto, __oneof_names_TypeProto, ProtoBuf.DEF_FIELD_TYPES)
hash(v::TypeProto) = ProtoBuf.protohash(v)
isequal(v1::TypeProto, v2::TypeProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::TypeProto, v2::TypeProto) = ProtoBuf.protoeq(v1, v2)

mutable struct ValueInfoProto <: ProtoType
    name::AbstractString
    _type::TypeProto
    doc_string::AbstractString
    ValueInfoProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ValueInfoProto
hash(v::ValueInfoProto) = ProtoBuf.protohash(v)
isequal(v1::ValueInfoProto, v2::ValueInfoProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::ValueInfoProto, v2::ValueInfoProto) = ProtoBuf.protoeq(v1, v2)

mutable struct OperatorSetIdProto <: ProtoType
    domain::AbstractString
    version::Int64
    OperatorSetIdProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct OperatorSetIdProto
hash(v::OperatorSetIdProto) = ProtoBuf.protohash(v)
isequal(v1::OperatorSetIdProto, v2::OperatorSetIdProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::OperatorSetIdProto, v2::OperatorSetIdProto) = ProtoBuf.protoeq(v1, v2)

struct __enum_AttributeProto_AttributeType <: ProtoEnum
    UNDEFINED::Int32
    FLOAT::Int32
    INT::Int32
    STRING::Int32
    TENSOR::Int32
    GRAPH::Int32
    FLOATS::Int32
    INTS::Int32
    STRINGS::Int32
    TENSORS::Int32
    GRAPHS::Int32
    __enum_AttributeProto_AttributeType() = new(0,1,2,3,4,5,6,7,8,9,10)
end #struct __enum_AttributeProto_AttributeType
const AttributeProto_AttributeType = __enum_AttributeProto_AttributeType()

mutable struct AttributeProto <: ProtoType
    name::AbstractString
    doc_string::AbstractString
    _type::Int32
    f::Float32
    i::Int64
    s::Array{UInt8,1}
    t::TensorProto
    g::Any
    floats::Vector{Float32}
    ints::Vector{Int64}
    strings::Vector{Array{UInt8,1}}
    tensors::Vector{TensorProto}
    graphs::Any
    AttributeProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AttributeProto
const __fnum_AttributeProto = Int[1,13,20,2,3,4,5,6,7,8,9,10,11]
const __ftype_AttributeProto = Dict(:g => "GraphProto", :graphs => "Vector{GraphProto}")
meta(t::Type{AttributeProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_AttributeProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, __ftype_AttributeProto)
hash(v::AttributeProto) = ProtoBuf.protohash(v)
isequal(v1::AttributeProto, v2::AttributeProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::AttributeProto, v2::AttributeProto) = ProtoBuf.protoeq(v1, v2)

mutable struct NodeProto <: ProtoType
    input::Vector{AbstractString}
    output::Vector{AbstractString}
    name::AbstractString
    op_type::AbstractString
    domain::AbstractString
    attribute::Vector{AttributeProto}
    doc_string::AbstractString
    NodeProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct NodeProto
const __fnum_NodeProto = Int[1,2,3,4,7,5,6]
meta(t::Type{NodeProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_NodeProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)
hash(v::NodeProto) = ProtoBuf.protohash(v)
isequal(v1::NodeProto, v2::NodeProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::NodeProto, v2::NodeProto) = ProtoBuf.protoeq(v1, v2)

mutable struct GraphProto <: ProtoType
    node::Vector{NodeProto}
    name::AbstractString
    initializer::Vector{TensorProto}
    doc_string::AbstractString
    input::Vector{ValueInfoProto}
    output::Vector{ValueInfoProto}
    value_info::Vector{ValueInfoProto}
    GraphProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct GraphProto
const __fnum_GraphProto = Int[1,2,5,10,11,12,13]
meta(t::Type{GraphProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_GraphProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)
hash(v::GraphProto) = ProtoBuf.protohash(v)
isequal(v1::GraphProto, v2::GraphProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::GraphProto, v2::GraphProto) = ProtoBuf.protoeq(v1, v2)

mutable struct ModelProto <: ProtoType
    ir_version::Int64
    opset_import::Vector{OperatorSetIdProto}
    producer_name::AbstractString
    producer_version::AbstractString
    domain::AbstractString
    model_version::Int64
    doc_string::AbstractString
    graph::GraphProto
    metadata_props::Vector{StringStringEntryProto}
    ModelProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct ModelProto
const __fnum_ModelProto = Int[1,8,2,3,4,5,6,7,14]
meta(t::Type{ModelProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_ModelProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)
hash(v::ModelProto) = ProtoBuf.protohash(v)
isequal(v1::ModelProto, v2::ModelProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::ModelProto, v2::ModelProto) = ProtoBuf.protoeq(v1, v2)

export Version, AttributeProto_AttributeType, AttributeProto, ValueInfoProto, NodeProto, ModelProto, StringStringEntryProto, GraphProto, TensorProto_DataType, TensorProto_Segment, TensorProto, TensorShapeProto_Dimension, TensorShapeProto, TypeProto_Tensor, TypeProto, OperatorSetIdProto, AttributeProto_AttributeType, AttributeProto

end


text:
export printpoly

## Poly{T} is basically T[x], with T a Ring.
## T[x] may not have an order so abs, comparing to 0 may not be defined.

## to handle this case we create some functions
## which can be modified by users for other Ts

_iszero(x::T) where {T} = (x == zero(T)) === true
_iszero(x::AbstractArray{T}) where {T} = all(isequal.(x, zero(T)))

_isone(x::T) where {T} = (x == one(T)) === true
_isone(x::AbstractArray{T}) where {T} = all(isequal.(x, one(T)))

"`hasneg(::T)` attribute is true if: `pj < zero(T)` is defined."
hasneg(::Type{T}) where {T} = false

"Could value possibly be negative and if so, is it?"
isneg(pj::T) where {T} = hasneg(T) && _isone(sign(-pj))

"Make `pj` positive if it is negative. (Don't call `abs` as that may not be defined, or appropriate.)"
aspos(pj::T) where {T} = (hasneg(T) && isneg(pj)) ? -pj : pj

"Should a value of `one(T)` be shown as a coefficient of monomial `x^i`, `i >= 1`? (`1.0x^2` is shown, `1 x^2` is not)"
showone(::Type{T}) where {T} = true

#####

## Numbers
hasneg(::Type{T}) where {T<:Real} = true

### Integer
showone(::Type{T}) where {T<:Integer} = false
showone(::Type{Rational{T}}) where {T<:Integer} = false

### Complex coefficients
hasneg(::Type{Complex{T}}) where {T} = true      ## we say neg if real(z) < 0 || real(z) == 0 and imag(g) < 0

function isneg(pj::Complex{T}) where {T}
    real(pj) < 0 && return true
    (real(pj) == 0 && imag(pj) < 0) && return(true)
    return false
end

showone(pj::Type{Complex{T}}) where {T} = showone(T)

### Polynomials as coefficients
hasneg(::Type{<:AbstractPolynomial{S}}) where {S} = false
showone(::Type{<:AbstractPolynomial{S}}) where {S} = false

#=
Common Printing
=#

Base.show(io::IO, p::AbstractPolynomial) = show(io, MIME("text/plain"), p)

function Base.show(io::IO, mimetype::MIME"text/plain", p::P) where {P<:AbstractPolynomial}
    print(io,"$(P.name.wrapper)(")
    printpoly(io, p, mimetype)
    print(io,")")
end

function Base.show(io::IO, mimetype::MIME"text/latex", p::AbstractPolynomial)
    print(io, "\$")
    printpoly(io, p, mimetype)
    print(io, "\$")
end

function Base.show(io::IO, mimetype::MIME"text/html", p::AbstractPolynomial)
    printpoly(io, p, mimetype)
end

# print uses compact representation
Base.print(io::IO, p::AbstractPolynomial) = printpoly(io, p, MIME("text/plain"), compact=true)

#####

"Show different operations depending on mimetype. `l-` is leading minus sign."
function showop(::MIME"text/plain", op)
    d = Dict("*" => "*", "+" => " + ", "-" => " - ", "l-" => "-")
    get(d, op, "")
end

function showop(::MIME"text/latex", op)
    d = Dict("*" => "\\cdot ", "+" => " + ", "-" => " - ", "l-" => "-")
    get(d, op, "")
end

function showop(::MIME"text/html", op)
    d = Dict("*" => "∙", "+" => " + ", "-" => " - ", "l-" => "-")
    get(d, op, "")
end

###

"""
    printpoly(io::IO, p::AbstractPolynomial, mimetype = MIME"text/plain"(); descending_powers=false, offset::Int=0, var=indeterminate(p), compact=false, mulsymbol="*")

Print a human-readable representation of the polynomial `p` to `io`. The MIME
types "text/plain" (default), "text/latex", and "text/html" are supported. By
default, the terms are in order of ascending powers, matching the order in
`coeffs(p)`; specifying `descending_powers=true` reverses the order.
`offset` allows for an integer number to be added to the exponent, just for printing.
`var` allows for overriding the variable used for printing. Setting `mulsymbol=""`
will avoid an operator being printed. Setting `compact=true` will use a compact style for floating point numbers.

# Examples
```jldoctest show
julia> using Polynomials

julia> printpoly(stdout, Polynomial([1,2,3], :y))
1 + 2*y + 3*y^2

julia> printpoly(stdout, Polynomial([1,2,3], :y), descending_powers=true)
3*y^2 + 2*y + 1

julia> printpoly(stdout, Polynomial([2, 3, 1], :z), descending_powers=true, offset=-2)
1 + 3*z^-1 + 2*z^-2

julia> printpoly(stdout, Polynomial([-1, 0, 1], :z), offset=-1, descending_powers=true)
z - z^-1

julia> printpoly(stdout, Polynomial([-1, 0, 1], :z), offset=-1, descending_powers=true, var=:x)
x - x^-1

julia> p = Polynomial([sqrt(i) for i in 1:4])
Polynomial(1.0 + 1.4142135623730951*x + 1.7320508075688772*x^2 + 2.0*x^3)

julia> printpoly(stdout, p, compact=true)
1.0 + 1.41421*x + 1.73205*x^2 + 2.0*x^3

julia> printpoly(stdout, map(x -> round(x, digits=12), p))  # more control on rounding
1.0 + 1.414213562373*x + 1.732050807569*x^2 + 2.0*x^3
```
"""
function printpoly(io::IO, p::P, mimetype=MIME"text/plain"();
                   descending_powers=false, offset::Int=0, var=indeterminate(p),
                   compact=false, mulsymbol="*") where {T,P<:AbstractPolynomial{T}}
    first = true
    printed_anything = false
    for i in (descending_powers ? reverse(eachindex(p)) : eachindex(p))
        ioc = IOContext(io,
                        :compact=>get(io, :compact, compact),
                        :multiplication_symbol => get(io, :multiplication_symbol, mulsymbol)
                        )
        printed = showterm(ioc, P, p[i], var, i+offset, first, mimetype)
        first &= !printed
        printed_anything |= printed
    end
    printed_anything || print(io, zero(T))
    return nothing
end

"""
    showterm(io::IO, ::Type{<:AbstractPolynomial} pj, var, j, first, mimetype)

Shows the j'th term of the given polynomial. Returns `true` after successfully printing.

For example. for a `Polynomial` this would show the term `pj * var^j`.
"""
function showterm(io::IO, ::Type{AbstractPolynomial}, pj::T, var, j, first::Bool, mimetype) where {T} end

## print the sign
## returns aspos(pj)
function printsign(io::IO, pj::T, first, mimetype) where {T}
    neg = hasneg(T) && isneg(pj) === true
    if first
        neg && print(io, showop(mimetype, "l-"))    #Prepend - if first and negative
    else
        neg ? print(io, showop(mimetype, "-")) : print(io,showop(mimetype, "+"))
    end
    aspos(pj)
end

## print * or cdot, ...
## pass `:multiplication_symbol => "" to IOContext have no sign
function printproductsign(io::IO, pj::T, j, mimetype) where {T}
    j == 0 && return
    multiplication_symbol = showop(mimetype, get(io, :multiplication_symbol,"*"))
    (showone(T) || !_isone(pj)) &&  print(io, multiplication_symbol)
end

function printproductsign(io::IO, pj::T, j, mimetype) where {T<:Complex}
    j == 0 && return
    (a,b) = reim(pj)
    !iszero(a) && !iszero(b) && return # parentheses inserted, no * needed
    !iszero(a) && return printproductsign(io, a, j, mimetype)
    print(io, showop(mimetype, "*"))
end

# show a single term
# Other types can overload Polynomials.printcofficient with a mimetype
# or Base.show_unquoted(io, pj, indent, prec)
# For example

"""
    printcoefficient(io::IO, pj, j, mimetype)

Print coefficient pj of monomial pj * x^j with the given mimetype.

For pretty printing different number types, or for adding parentheses,
methods can be added to this function. If no mimetype is desired,
adding a method to `Base.show_unquoted` is suggested, as this will
also be useful for the default `show` methods. The following example
shows how `Dual` objects of `DualNumbers` may be printed with
parentheses.

```jldoctest
julia> using Polynomials, DualNumbers

julia> Polynomial([Dual(1,2), Dual(3,4)])
Polynomial(1 + 2ɛ + 3 + 4ɛ*x)
```

```jldoctest
julia> using DualNumbers, Polynomials

julia> function Base.show_unquoted(io::IO, pj::Dual, indent::Int, prec::Int)
            if Base.operator_precedence(:+) <= prec
                print(io, "(")
                show(io, pj)
                print(io, ")")
            else
                show(io, pj)
            end
        end

julia> Polynomial([Dual(1,2), Dual(3,4)])
Polynomial((1 + 2ɛ) + (3 + 4ɛ)*x)
```
"""
printcoefficient(io::IO, pj::Any, j, mimetype) = Base.show_unquoted(io, pj, 0, Base.operator_precedence(:*))

# pretty print rational numbers in latex
function printcoefficient(io::IO, a::Rational{T}, j, mimetype::MIME"text/latex") where {T}
    abs(a.den) == one(T) ? print(io, a.num) : print(io, "\\frac{$(a.num)}{$(a.den)}")
end

# print complex numbers with parentheses as needed
function printcoefficient(io::IO, pj::S, j, mimetype) where {T,S <: Complex{T}}

    (a,b) = reim(pj)
    hasreal = !iszero(a) || isnan(a) || isinf(a)
    hasimag = !iszero(b) || isnan(b) || isinf(b)

    if hasreal && hasimag
        iszero(j) || print(io, "(")
        print(io, a)

        # print b
        if isone(b) || isone(-b)
            if hasneg(S) && b < 0
                print(io, showop(mimetype, "-"))
            else
                print(io, showop(mimetype, "+"))
            end
        else
            if hasneg(S) && b < 0
                print(io, showop(mimetype, "-"))
                (showone(S) || !isone(-b)) && print(io, -b)
            else
                print(io, showop(mimetype,"+"))
                print(io, b)
            end
            (isnan(b) || isinf(b)) && print(io, showop(mimetype, "*"))
        end

        print(io, imagsymbol(mimetype))
        iszero(j) || print(io, ")")

    elseif hasreal

        (iszero(j) || showone(T) || !isone(a)) && printcoefficient(io, a, j, mimetype)

    elseif hasimag
        (showone(T) || !isone(b)) && printcoefficient(io, b, j, mimetype)
        (isnan(b) || isinf(b)) && print(io, showop(mimetype, "*"))
        print(io, imagsymbol(mimetype))

    end

    return nothing

end

## show exponent

exponent_text(i, ::MIME) = "^$(i)"
exponent_text(i, ::MIME"text/html") = "<sup>$(i)</sup>"
exponent_text(i, ::MIME"text/latex") = "^{$(i)}"

function printexponent(io, var, i, mimetype::MIME)
    if i == 0
        return
    elseif i == 1
        print(io,var)
    else
        print(io, var, exponent_text(i, mimetype))
    end
end

function unicode_exponent(io, j)
    a = ("⁻","","","⁰","¹","²","³","⁴","⁵","⁶","⁷","⁸","⁹")
    for i in string(j)
        print(io, a[Int(i)-44])
    end
end

function unicode_subscript(io, j)
    a = ("₋","","","₀","₁","₂","₃","₄","₅","₆","₇","₈","₉")
    for i in string(j)
        print(io, a[Int(i)-44])
    end
end

imagsymbol(::Any) = "im"
imagsymbol(::MIME"text/latex") = "i"


text:
@testset "OpenSpielEnv" begin

    for name in [
        "tic_tac_toe",
        "kuhn_poker",
        "goofspiel(imp_info=True,num_cards=4,points_order=descending)",
    ]
        @info "testing OpenSpiel: $name"
        env = OpenSpielEnv(name)
        RLBase.test_runnable!(env)
    end
end


text:
@syms 𝑥
@syms u[20]

mutable struct Transform
    k::Int
    sub::Dict
end

function next_variable!(f, eq)
    μ = u[f.k]
    f.k += 1
    f.sub[μ] = eq
    return μ
end

transformer(eq::SymbolicUtils.Add, f) = sum(transformer(t,f) for t in arguments(eq); init=0)
transformer(eq::SymbolicUtils.Mul, f) = prod(transformer(t,f) for t in arguments(eq); init=1)
transformer(eq::SymbolicUtils.Div, f) = transformer(arguments(eq)[1],f) * transformer(inv(arguments(eq)[2]),f)

function transformer(eq::SymbolicUtils.Pow, f)
    y, k = arguments(eq)

    if is_pos_int(k)
        μ = next_variable!(f, y)
        return μ ^ k
    elseif is_neg_int(k)
        μ = next_variable!(f, inv(y))
        return μ ^ -k
    else
        return next_variable!(f, y^k)
    end
end

function transformer(eq, f)
    if isdependent(eq, 𝑥)
        return next_variable!(f, eq)
    else
        return 1
    end
end

function transform(eq, x)
    eq = substitute(eq, Dict(x => 𝑥))
    f = Transform(1, Dict())
    q = transformer(eq, f)
    if !any(is_poly, values(f.sub))
        q *= next_variable!(f, 1)
    end
    return q, f.sub
end

##############################################################################

function substitute_x(eq, x, sub)
    eq = substitute(eq, sub)
    substitute(eq, Dict(𝑥 => x))
end

function generate_homotopy(eq, x)
    q, sub = transform(eq, x)
    S = 0

    for i = 1:length(sub)
        μ = u[i]
        h₁, ∂h₁ = apply_partial_int_rules(sub[μ])
        h₂ = expand_derivatives(Differential(μ)(q))

        h₁ = substitute_x(h₁, x, sub)
        h₂ = substitute_x(h₂ * ∂h₁^-1, x, sub)

        S += expand((1 + h₁) * (1 + h₂))
    end

    unique([one(x); [equivalent(t,x) for t in terms(S)]])
end

##############################################################################

function ∂(x)
    d = expand_derivatives(Differential(𝑥)(x))
    return isequal(d, 0) ? 1 : d
end

partial_int_rules = [
    @rule 𝛷(^(sin(~x), ~k::is_neg)) => 𝛷(^(csc(~x), -~k))
    @rule 𝛷(^(cos(~x), ~k::is_neg)) => 𝛷(^(sec(~x), -~k))
    @rule 𝛷(^(tan(~x), ~k::is_neg)) => 𝛷(^(cot(~x), -~k))
    @rule 𝛷(^(csc(~x), ~k::is_neg)) => 𝛷(^(sin(~x), -~k))
    @rule 𝛷(^(sec(~x), ~k::is_neg)) => 𝛷(^(cos(~x), -~k))
    @rule 𝛷(^(cot(~x), ~k::is_neg)) => 𝛷(^(tan(~x), -~k))

    @rule 𝛷(^(sinh(~x), ~k::is_neg)) => 𝛷(^(csch(~x), -~k))
    @rule 𝛷(^(cosh(~x), ~k::is_neg)) => 𝛷(^(sech(~x), -~k))
    @rule 𝛷(^(tanh(~x), ~k::is_neg)) => 𝛷(^(coth(~x), -~k))
    @rule 𝛷(^(csch(~x), ~k::is_neg)) => 𝛷(^(sinh(~x), -~k))
    @rule 𝛷(^(sech(~x), ~k::is_neg)) => 𝛷(^(cosh(~x), -~k))
    @rule 𝛷(^(coth(~x), ~k::is_neg)) => 𝛷(^(tanh(~x), -~k))

    @rule 𝛷(sin(~x)) => (cos(~x), ∂(~x))
    @rule 𝛷(cos(~x)) => (sin(~x), ∂(~x))
    @rule 𝛷(tan(~x)) => (log(cos(~x)), ∂(~x))
    @rule 𝛷(csc(~x)) => (log(sin(~x)^-1 - cos(~x)*sin(~x)^-1), ∂(~x))
    @rule 𝛷(sec(~x)) => (log(cos(~x)^-1 + sin(~x)*cos(~x)^-1), ∂(~x))
    @rule 𝛷(cot(~x)) => (log(sin(~x)), ∂(~x))

    @rule 𝛷(sinh(~x)) => (cosh(~x), ∂(~x))
    @rule 𝛷(cosh(~x)) => (sinh(~x), ∂(~x))
    @rule 𝛷(tanh(~x)) => (log(cosh(~x)), ∂(~x))
    @rule 𝛷(csch(~x)) => (log(sinh(~x)^-1 + cosh(~x)*sinh(~x)^-1), ∂(~x))
    @rule 𝛷(sech(~x)) => (log(cosh(~x)^-1 + sinh(~x)*cosh(~x)^-1), ∂(~x))
    @rule 𝛷(coth(~x)) => (log(sinh(~x)), ∂(~x))

    @rule 𝛷(asin(~x)) => (~x*asin(~x) + sqrt(1 - ~x*~x), ∂(~x))
    @rule 𝛷(acos(~x)) => (~x*acos(~x) + sqrt(1 - ~x*~x), ∂(~x))
    @rule 𝛷(atan(~x)) => (~x*atan(~x) + log(~x*~x + 1), ∂(~x))
    @rule 𝛷(acsc(~x)) => (~x*acsc(~x) + acosh(~x), ∂(~x))     # needs an abs inside acosh
    @rule 𝛷(asec(~x)) => (~x*asec(~x) + acosh(~x), ∂(~x))     # needs an abs inside acosh
    @rule 𝛷(acot(~x)) => (~x*acot(~x) + log(~x*~x + 1), ∂(~x))

    @rule 𝛷(asinh(~x)) => (~x*asinh(~x) + sqrt(~x*~x + 1), ∂(~x))
    @rule 𝛷(acosh(~x)) => (~x*acosh(~x) + sqrt(~x*~x - 1), ∂(~x))
    @rule 𝛷(atanh(~x)) => (~x*atanh(~x) + log(~x + 1), ∂(~x))
    @rule 𝛷(acsch(~x)) => (acsch(~x), ∂(~x))
    @rule 𝛷(asech(~x)) => (asech(~x), ∂(~x))
    @rule 𝛷(acoth(~x)) => (~x*acot(~x) + log(~x + 1), ∂(~x))

    @rule 𝛷(log(~x)) => (~x + ~x * log(~x), ∂(~x))

    @rule 𝛷(^(~x, ~k::is_abs_half)) => (sum(candidate_sqrt(~x,~k); init=one(~x)), 1)
    @rule 𝛷(^(~x::is_poly, ~k::is_neg)) => (sum(candidate_pow_minus(~x, ~k); init=one(~x)), 1)
    @rule 𝛷(sqrt(~x)) => (sum(candidate_sqrt(~x,0.5); init=one(~x)), 1)
    @rule 𝛷(^(sqrt(~x),-1)) => 𝛷(^(~x,-0.5))

    @rule 𝛷(^(~x, -1)) => (log(~x), ∂(~x))
    @rule 𝛷(1 / ~x) => 𝛷(^(~x, -1))
    @rule 𝛷(^(~x, ~k)) => (^(~x, ~k+1), ∂(~x))

    @rule 𝛷(exp(~x)) => (exp(~x), ∂(~x))
    @rule 𝛷(1) => (𝑥, 1)
    @rule 𝛷(~x) => ((~x + ^(~x,2)), ∂(~x))
]

apply_partial_int_rules(eq) = expand(Fixpoint(Prewalk(Chain(partial_int_rules))))(𝛷(value(eq)))


text:
struct WMSlayer
	name::String
	title::String
	srs::String
	crs::String
	bbox::NTuple{4,Float64}
	imgformat::String
	tilesize::Int
	overviewcount::Int
	resolution::Float64
	tiled::Bool
	transparent::Bool
end
struct WMS
	serverURL::String
	OnlineResource::String
	version::String
	request::String
	layernames::Vector{String}
	layer::Vector{WMSlayer}
end

# ---------------------------------------------------------------------------------------------------
"""
    wmsinfo(server::String)

Read the xml information from the WebMapServer service and create a WMS data type that holds the information
necessary to download the data. The ``show`` method will display the contents of the WMS data type.

- `server`: The server URL address.

## Example
    wms = wmsinfo("http://tiles.maps.eox.at/wms?")

As an option, use the form

    wmsinfo(wms::WMS; layer)

to get further information, in particular the number of bands and sizes, of the layer number or layer name `layer`. `wms` is
returned by the first form.
"""
function wmsinfo(server::String)::WMS
	w = gdalinfo("WMS:" * server)
	(w === nothing) && error("unable to open '$(server)'")
	sds = split(w, "SUBDATASET_")
	lastguy = collect(2:2:length(sds))[end]

	names, layers = String[], WMSlayer[]

	onlineRes::String, ver::SubString{String}, req::SubString{String}, name::SubString{String}, srs::SubString{String}, crs::SubString{String}, bbox, fmt::SubString{String}, ts, cnt, res, til, trans = "", "", "", "", "", "", zeros(4), "", 0, 0, 0., false, false
	for k = 2:2:length(sds)
		f = split(sds[k], '&')
		if (k == 2)					# "f[1] = 1_NAME=WMS:http://tiles.maps.eox.at/?SERVICE=WMS"
			ind = findfirst("=", f[1]); t = f[1][ind[1]+5:end]
			ind = findfirst("?", t);	onlineRes = t[1:ind[1]]
		end
		for n = 2:length(f)			# Loop over the fields od the SUBDATASET_XX_NAME
			if     (startswith(f[n], "VERSION"))  ver = f[n][9:end];  continue
			elseif (startswith(f[n], "REQUEST"))  req = f[n][9:end];  continue
			elseif (startswith(f[n], "LAYERS"))   name = f[n][8:end]; continue
			elseif (startswith(f[n], "SRS"))      srs = f[n][5:end];  continue
			elseif (startswith(f[n], "CRS"))      crs = f[n][5:end];  continue
			elseif (startswith(f[n], "BBOX"))     bbox = parse.(Float64, split(f[n][6:end], ',')); continue
			elseif (startswith(f[n], "FORMAT"))         fmt = f[n][8:end];	continue
			elseif (startswith(f[n], "TILESIZE"))       ts  = parse(Int, f[n][10:end]);	continue
			elseif (startswith(f[n], "OVERVIEWCOUNT"))  cnt = parse(Int, f[n][15:end]);	continue
			elseif (startswith(f[n], "MINRESOLUTION"))  res = parse(Float64, f[n][15:end]);	continue
			elseif (startswith(f[n], "TILED"))          til = contains(f[n], "true");	continue
			elseif (startswith(f[n], "TRANSPARENT"))    trans = contains(f[n], "true");	continue
			end
		end
		ind = findfirst("=", sds[k+1])
		title = (ind !== nothing) ? chomp(rstrip(sds[k+1][ind[1]+1:end])) : ""	# It comes with a damn "...\n  "
		if (k == lastguy)			# Last one contains trailing trash with fake "Corner Coordinates: ..."
			ind = findfirst("Corner", title)
			title = title[1:ind[1]-2]
		end
		append!(names, [name])
		append!(layers, [WMSlayer(name, title, srs, crs, (bbox[1],bbox[3],bbox[2],bbox[4]), fmt, ts, cnt, res, til, trans)])
		name, srs, crs, bbox, fmt, ts, cnt, res, til, trans = "", "", "", zeros(4), "", 0, 0, 0., false, false
	end
	WMS(server, onlineRes, ver, req, names, layers)
end

# ---------------------------------------------------------------------------------------------------
function wmsinfo(wms::WMS; layer=0, stronly::Bool=false)
	# This method inquires a specific layer
	layer_n = get_layer_number(wms, layer)
	nt = (region=wms.layer[layer_n].bbox, size=(10,10))
	str, = wms_helper(wms; layer=layer_n, nt...)
	println(str)
	(!stronly) && println(gdalinfo(str))
	return nothing
end

# ---------------------------------------------------------------------------------------------------
"""
    wmsread(wms::WMS; layer, kwargs...)

Read the `layer` number provided by the service from which the `wms` type was created.

### Parameters
- `wms`: A WMS type obtained from the ``wmsinfo`` function.
- `layer`: The layer number or layer name of interest from those provided by the WMS service. That is,
   both of these forms are allowed: `layer=3` or `layer="Invented layer name"`

### `kwargs` is the keywords/values pairs used to set
- `region | limits`: The region limits. This can be a Tuple or Array with 4 elements defining the ``(xmin, xmax, ymin, ymax)``
   or a string defining the limits in all ways that GMT can recognize. When the layer has the data projected, we can
   a Tuple or Array with 3 elements ``(lon0, lat0, width)``, where first two set the center of a square in geographical
   coordinates and the third (``width``) is the width of that box in meters.
- `cellsize | pixelsize | resolution | res`: Sets the requested cell size in meters [default]. Use a string appended with a 'd'
   (e.g. `resolution="0.001d"`) if the resolution is given in degrees. This way works only when the layer is in geogs.
- `size`: Alternatively to the the `cellsize` use this option, a tuple or array with two elements, to specify
   the image dimensions. Example, `size=(1200, 100)` to get an image with 1200 rows and 100 columns.

### Returns
A GMTimage

## Example
    img = wmsread(wms, layer=3, region=(-10,-5,37,44), pixelsize=500);
"""
function wmsread(wms::WMS; layer=0, kw...)
	layer_n = get_layer_number(wms, layer)
	str, dim_x, dim_y = wms_helper(wms; layer=layer_n, kw...)
	opts = ["-outsize", "$(dim_x)", "$(dim_y)"]
	gdaltranslate(str, opts)
end

# ---------------------------------------------------------------------------------------------------
"""
    wmstest(wms::WMS; layer, size::Bool=false, kwargs...)

Test function that generates the GetMap request string or the size of the resulting image
given the requested resolution. It is meant to generate first the command that gets the image/grid
but not running it. Specially usefull to check that the resulting image size is not huge.

### Parameters
- `wms`: A WMS type obtained from the ``wmsinfo`` function.
- `layer`: The layer number or layer name of interest from those provided by the WMS service. That is,
   both of these forms are allowed: `layer=3` or `layer="Invented layer name"`
- `size`: If ``false``, returns the GetMap request string, otherwise the image size given the requested resolution.

### `kwargs` is the keywords/values pairs used to set
- `region | limits`: The region limits. This can be a Tuple or Array with 4 elements defining the ``(xmin, xmax, ymin, ymax)``
   or a string defining the limits in all ways that GMT can recognize. When the layer has the data projected, we can
   a Tuple or Array with 3 elements ``(lon0, lat0, width)``, where first two set the center of a square in geographical
   coordinates and the third (``width``) is the width of that box in meters.
- `cellsize | pixelsize | resolution | res`: Sets the requested cell size in meters [default]. Use a string appended with a 'd'
   (e.g. `resolution="0.001d"`) if the resolution is given in degrees. This way works only when the layer is in geogs.

### Returns
Either a the GetMap request string (when `size=false`) or the resulting image/grid size `dim_x, dim_y`

## Example
    wmstest(wms, layer=34, region=(-8,39, 100000), pixelsize=100)

	"WMS:http://tiles.maps.eox.at/?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=s2cloudless-2020_3857&SRS=EPSG:900913&BBOX=-940555.9,4671671.6,-840555.9,4771671.68&FORMAT=image/jpeg&TILESIZE=256&OVERVIEWCOUNT=18&MINRESOLUTION=0.59716428347794&TILED=true"
"""
function wmstest(wms::WMS; layer=0, size::Bool=false, kw...)
	layer_n = get_layer_number(wms, layer)
	str, dim_x, dim_y = wms_helper(wms; layer=layer_n, kw...)
	return (size) ? (dim_y, dim_x) : str
end

# ---------------------------------------------------------------------------------------------------
function get_layer_number(wms::WMS, layer::Union{Int, AbstractString})::Int
	# Return the layer number in case 'layer' is a string
	if (isa(layer, AbstractString))
		layer_n = findfirst(wms.layernames .== layer)
		(layer_n === nothing) && error("Layer name not found in this WMS service")
	else
		layer_n = layer
	end
	(layer_n <= 0) && error("Must set the layer number or layer name that you want information about.")
	(layer_n > length(wms.layer)) && error("The requested layer is greater then the number of layers in this dataset.")
	layer_n
end

# ---------------------------------------------------------------------------------------------------
function wms_helper(wms::WMS; layer=0, kw...)
	# Helper function that creates the Getmap string and also optionaly computes the image size from resolution.
	layer_n = get_layer_number(wms, layer)
	d = KW(kw)
	((reg = find_in_dict(d, [:R :region :limits], false)[1]) === nothing) && error("Must provide the `region` option.")

	if (isa(reg, Tuple) || isa(reg, VMr))
		len::Int = length(reg)				# reg is a damn Any
		(len > 4 || len < 3) && error("The region array must have THREE or FOUR elements.")
		if (len == 4)
			lims::Vector{<:Real} = [reg[1], reg[2], reg[3], reg[4]]		# Make a copy to have the same name as in the other branch.
		else
			(contains(wms.layer[layer_n].srs, "4326") || contains(wms.layer[layer_n].crs, ":84")) &&
				error("This is a Geographical layer so you cannot set the region with a center and a width.")
			t_srs = epsg2proj(parse(Int, wms.layer[layer_n].srs[6:end]))	# Convert to proj4 because lonlat2xy() does not know EPSG
			center::Matrix{Float64} = lonlat2xy([reg[1] reg[2]], t_srs)		# Box center in projected coordinates
			half_w::Float64 = reg[3] / 2
			lims = [center[1]-half_w, center[1]+half_w, center[2]-half_w, center[2]+half_w]
		end
		BB = @sprintf("%.12g,%.12g,%.12g,%.12g", lims[1], lims[3], lims[2], lims[4])
	else
		_, opt_R = parse_R(d, "")
		lims = opt_R2num(opt_R)
		BB = @sprintf("%.12g,%.12g,%.12g,%.12g", lims[1], lims[3], lims[2], lims[4])
	end
	(wms.layer[layer_n].bbox[1] > lims[1] || wms.layer[layer_n].bbox[2] < lims[2] || wms.layer[layer_n].bbox[3] > lims[3] || wms.layer[layer_n].bbox[4] < lims[4]) && @warn("Requested region overflows this layer BoundingBox")

	function getsize(res::Float64, lims::Vector{<:Real})::Tuple{Int, Int}
		_dim_x = round(Int, (lims[2] - lims[1]) / res)		# Assume data has pixel registration
		_dim_y = round(Int, (lims[4] - lims[3]) / res)
		return _dim_x, _dim_y
	end

	if ((dim = find_in_dict(d, [:size])[1]) !== nothing)	# Got the image size directly
		(!isa(dim, Tuple) && !isa(dim, VMr) && length(dim) != 2) &&
			error("The 'size' option must be a Tuple or Array with TWO elements")
		dim_x::Int, dim_y::Int = dim[2], dim[1]
	elseif ((res = find_in_dict(d, [:res :cellsize :pixelsize :resolution])[1]) !== nothing)
		if (isa(res, String) && res[end] == 'd')	# Got resolution in degrees (a string ending with a 'd')
			(!contains(wms.layer[layer_n].srs, "4326") && !contains(wms.layer[layer_n].crs, ":84")) &&
				error("This is not a Geographical dataset so you cannot set the resolution in degrees")
			_res::Float64 = parse(Float64, res[1:end-1])
			dim_x, dim_y = getsize(_res, lims)
		else
			if (!contains(wms.layer[layer_n].srs, "4326") && !contains(wms.layer[layer_n].crs, ":84"))	# Easy case
				_res = res
			else				# Here BB is in degrees and resolution in meters
				deg_per_m = 360 / (6371000 * 2pi)		# Use spherical authalic radius
				_res = res * deg_per_m
			end
			dim_x, dim_y = getsize(_res, lims)
		end
	else
		error("Must provide either the 'size' or the 'resolution' option")
	end

	#WMS:http://tiles.maps.eox.at/?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=overlay_base_bright&SRS=EPSG:4326&BBOX=-180.000000,-90.000000,180.000000,90.000000&FORMAT=image/png&TILESIZE=256&OVERVIEWCOUNT=17&MINRESOLUTION=0.0000053644180298&TILED=true

	TILED = (wms.layer[layer_n].tiled) ? @sprintf("&TILED=%s", string(wms.layer[layer_n].tiled)) : ""
	TRANS = (wms.layer[layer_n].transparent) ? @sprintf("&TRANSPARENT=%s", string(wms.layer[layer_n].transparent)) : ""
	MINRES= (wms.layer[layer_n].resolution > 0.) ? @sprintf("&MINRESOLUTION=%.14f", wms.layer[layer_n].resolution) : ""
	OVCNT = (wms.layer[layer_n].overviewcount > 0) ? @sprintf("&OVERVIEWCOUNT=%d", wms.layer[layer_n].overviewcount) : ""
	TILSZ = (wms.layer[layer_n].tilesize > 0) ? @sprintf("&TILESIZE=%d", wms.layer[layer_n].tilesize) : ""
	FMT   = (wms.layer[layer_n].imgformat != "") ? @sprintf("&FORMAT=%s", wms.layer[layer_n].imgformat) : "&FORMAT=image/png"

	str = @sprintf("WMS:%sSERVICE=WMS&VERSION=%s&REQUEST=GetMap&LAYERS=%s&SRS=%s&BBOX=%s%s%s%s%s%s%s", wms.OnlineResource, wms.version, wms.layer[layer_n].name, wms.layer[layer_n].srs, BB, FMT, TILSZ, OVCNT, MINRES, TILED, TRANS)

	return str, dim_x, dim_y
end

# ---------------------------------------------------------------------------------------------------
function Base.:show(io::IO, W::WMS)
	println("serverURL:\t", W.serverURL);   println("OnlineResource:\t", W.OnlineResource)
	println("version:\t", W.version);       println("request:\t", W.request);
	println("layernames:\t", W.layernames)
	println("\nlayer:\t", length(W.layer), " Layers. Use layer[k] to see the contents of layer k")
end
function Base.:show(io::IO, W::WMSlayer)
	println("name:\t", W.name);		println("title:\t", W.title);		println("srs:\t", W.srs);
	println("crs:\t", W.crs);		println("bbox:\t", W.bbox)
	println("imgformat:\t", W.imgformat);          println("tilesize:\t", W.tilesize)
	println("overviewcount:\t", W.overviewcount);  println("resolution:\t", W.resolution)
	println("tiled:\t", W.tiled); println("transparent:\t", W.transparent)
end


text:
module MyDSP

export resample, bandpass, savitsky_golay, plot_spectrogram

using LinearAlgebra
using Plots
using Plots.PlotMeasures
using DSP

using ..MyUtils
#using ERec

# functions

function plot_spectrogram(s, fs, plt, sp)
    s = flatten(s)
    S = spectrogram(s, Int(25e-3*fs), Int(10e-3*fs); window=hanning)
    t = time(S)
    f = freq(S)
    #

    xtick_id = collect(1:100:size(S.power, 2))
    xtick_label = round.(Int.(collect(t))[xtick_id] ./ fs .*1000, digits=1)

    ytick_id = collect(1:50:size(S.power, 1))
    ytick_label = round.(collect(f)[ytick_id].*fs, digits=1)

    pl = heatmap!(log.(power(S)),
        xticks=( xtick_id, xtick_label ),
        yticks=( ytick_id, ytick_label ),
        grid=false,
        left_margin=50px,
        xlabel="ms",
        ylabel="Hz",
        subplot=sp
    );

    return pl
end

function savitsky_golay(x::Vector, windowSize::Integer, polyOrder::Integer; dv::Integer=0)
    #Polynomial smoothing with the Savitsky Golay filters
    #
    # Sources
    # ---------
    # Theory: http://www.ece.rutgers.edu/~orfanidi/intro2sp/orfanidis-i2sp.pdf
    # Python Example: http://wiki.scipy.org/Cookbook/SavitzkyGolay

    #Some error checking
    @assert isodd(windowSize) "Window size must be an odd integer."
    @assert polyOrder < windowSize "Polynomial order must me less than window size."
    halfWindow = Int((windowSize-1)/2)
    #Setup the S matrix of basis vectors.
    S = zeros(windowSize, polyOrder+1)
    for ct = 0:polyOrder
    	S[:,ct+1] = collect(-halfWindow:halfWindow).^(ct)
    end

    #Compute the filter coefficients for all orders
    #From the scipy code it seems pinv(S) and taking rows should be enough
    G = S*pinv(S'*S)

    #Slice out the derivative order we want
    filterCoeffs = G[:,dv+1] * factorial(dv);

    #Pad the signal with the endpoints and convolve with filter
    paddedX = [x[1]*ones(halfWindow), x, x[end]*ones(halfWindow)]
    y = conv(filterCoeffs[end:-1:1],  collect(Iterators.flatten(paddedX)))

    #Return the valid midsection
    return y[2*halfWindow+1:end-2*halfWindow]
end

end

max_stars_repo_name: maalaria/Newbe

max_issues_repo_name: maalaria/Newbe

max_forks_repo_name: maalaria/Newbe


text:
abstract type GaussianItem end

mutable struct Gaussian <: GaussianItem
    pi::Float64
    tau::Float64
end

function Gaussian(; mu::Float64=0.0, sigma::Float64=0.0)
    if sigma == 0
        return Gaussian(0.0, 0.0)
    end
    pi = sigma^(-2)
    tau = pi * mu
    Gaussian(pi, tau)
end

@inline function mu(g::Gaussian)
    if g.pi != 0
        return g.tau / g.pi
    else
        return 0
    end
end

@inline function sigma(g::Gaussian)
    if g.pi != 0
        return(sqrt(1/g.pi))
    else
        return Inf
    end
end

@inline function isless(a::Gaussian, b::Gaussian)
    isless(mu(a), mu(b))
end

@inline (*)(a::Gaussian, b::Gaussian) = Gaussian(a.pi + b.pi, a.tau + b.tau)
@inline (/)(a::Gaussian, b::Gaussian) = Gaussian(a.pi - b.pi, a.tau - b.tau)

function show(io::IO, g::Gaussian)
    print(io, "Gaussian(mu=$(mu(g)), sigma=$(sigma(g)))")
end


text:
using Dates, Random, UUIDs
using Logging

import Base: convert

export
    Model,
    Network,
    InstantNetwork,
    DynamicNetwork,
    Placeholder,
    Node,
    SFunc,
    Dist,
    Score,
    ValueTyped,
    Variable,
    VariableGraph,
    get_initial_graph,
    get_transition_graph,
    get_children,
    get_initial_children,
    get_transition_children,
    get_parents,
    get_initial_parents,
    get_transition_parents,
    get_sfunc,
    get_node,
    get_variables,
    get_placeholders,
    get_nodes,
    input_type,
    make_initial,
    make_transition,
    output_type

"""
    abstract type ValueTyped{O}

Supertype for typing all Scruff Variables; `O` is the
actual type of the variable
"""
abstract type ValueTyped{O} end

"""
value_type(v::ValueTyped{O}) where {O}

return the actual type (i.e. `O`) of the `ValueTyped`
"""
value_type(v::ValueTyped{O}) where {O} = O

"""
    abstract type SFunc{I<:Tuple, O}

A Stochastic Function type with input variables defined by `I` and output type `O`.
This is an abstract representation for a collection of operators with the same
input and output types.

All sfuncs have the following operators defined:
- `compute_lambda`
- `compute_bel`
- `send_pi`
- `outgoing_pis`
- `outgoing_lambdas`

`SFunc` _also_ has both the operators `cpdf` and `logcpdf` defined in terms of the other.
All sfuncs should implement one or the other of these operators.
"""
abstract type SFunc{I<:Tuple, O} end

Base.print_without_params(x::Type{<:SFunc}) = false

"""
    Dist{T} = SFunc{Tuple{}, T}

# Additional supported operators
- `make_factors`
- `send_lambda`
"""
Dist{T} = SFunc{Tuple{}, T}

"""
    Score{I} = SFunc{Tuple{I}, Nothing}

`Score` supports two (2) operators:  `get_score` and `get_log_score`.  `get_log_score`
is defined, by default, using `get_score.`  Every subtype of `Score` _must_ implement
`get_score`.
"""
Score{I} = SFunc{Tuple{I}, Nothing}

"""
    input_type(::Type{<:SFunc{I,O}}) where {I,O}
    input_type(::SFunc{I,O}) where {I,O}

Return the input type (i.e. `I`) of the `SFunc`
"""
input_type(::Type{<:SFunc{I,O}}) where {I,O} = I
input_type(::SFunc{I,O}) where {I,O} = I
"""
    output_type(::Type{<:SFunc{I,O}}) where {I,O}
    output_type(::SFunc{I,O}) where {I,O}

Return the output type (i.e. `O`) of the `SFunc`
"""
output_type(::Type{<:SFunc{I,O}}) where {I,O} = O
output_type(::SFunc{I,O}) where {I,O} = O

"""
    Placeholder{O} <: ValueTyped{O}

A type for typing Scruff variables that do not reference models
"""
struct Placeholder{O} <: ValueTyped{O}
    name::Symbol
end

"""
    Model{I, J, O} <: ValueTyped{O}

Supertype for all Scruff models.
The model represents a variable that varies over time and has output type `O`.
The model may return an initial sfunc with input type `I` using `make_initial`,
which takes the current time as argument,
and a transition sfunc with input type `J` using `make_transition`,
which takes both the parent times (a tuple of times of the same length as `J`)
and the current time as arguments.
These two functions need to be defined for every sfunc.

# Type parameters
- `I`: the input type to the `SFunc` returned by the model's `make_initial` function
- `J`: the input type to the `SFunc` used during the `make_trasition` function call
- `O`: the actual type of the variables represented by this model
"""
abstract type Model{I, J, O} <: ValueTyped{O} end
# is_fixed(m) = false

function make_initial(m,t) end

function make_transition(m, parenttimes, time) end

(model::Model)(symbol::Symbol) = instantiate(model, symbol)

"""
    instantiate(model::Model, name::Symbol)

Create a new `Variable` with the given name for the given model.  Every
`Model` instance is also a function that, given a `Symbol`, will call
this method.
"""
function instantiate(model::Model, name::Symbol)
    return Variable(name, model)
end

"""
    mutable struct Variable{I,J,O} <: ValueTyped{O}

A Variable describes the (time-series) of some set of random values.
It must be named, described by a model, and references to the model inputs must be defined.

For the type variables, see [`Model`](@ref)
"""
mutable struct Variable{I,J,O} <: ValueTyped{O}
    name::Symbol
    model::Model{I,J,O}
end

"""`Node{O} = Union{Placeholder{O}, Variable{I,J,O} where {I,J}}`"""
Node{O} = Union{Placeholder{O}, Variable{I,J,O} where {I,J}}

"""
    output_type(::Variable{I,J,O})
    output_type(::Placeholder{O})

returns the output type `O` of the given parameter
"""
function output_type(::Variable{I,J,O}) where {I,J,O}
    return O
end

function output_type(::Placeholder{O}) where O
    return O
end

get_name(node::Node) = node.name

"""VariableGraph = Dict{Node, Vector{Node}}"""
VariableGraph = Dict{Node, Vector{Node}}

"""
    abstract type Network{I,J,O}

Collects variables, with placeholders for inputs, and defined outputs,
along with up to two directed graphs for initial and transition.
Must implement `get_nodes`, `get_outputs`,
`get_initial_graph`, `get_transition_graph` (which returns a `VariableGraph`)

For the type parameters, see the underlying [`Model`](@ref) class for the mapped
[`Variable`](@ref)s.
"""
abstract type Network{I,J,O} end

"""
    get_initial_parents(n::Network, node::Node)::Vector{Node}

Returns a list of the initial parent nodes in the network from the network's
initial graph (i.e. the return graph from get_initial_graph(network)).
"""
function get_initial_parents(n::Network, node::Node)::Vector{Node}
    g = get_initial_graph(n)
    get(g, node, Node[])
end

"""
    get_transition_parents(n::Network, node::Node)::Vector{Node}

Returns a list of transitioned parent nodes in the network from the network's
transition graph (i.e. the return graph from get_transition_graph(network))
"""
function get_transition_parents(n::Network, node::Node)::Vector{Node}
    g = get_transition_graph(n)
    get(g, node, Node[])
end

"""
    get_node(n::Network, name::Symbol)::Union{Node, Nothing}

Returns the node with the given name, or `nothing`.
"""
function get_node(n::Network, name::Symbol)::Union{Node, Nothing}
    vs = get_nodes(n)
    for v in vs
        if v.name == name
            return v
        end
    end
    return nothing
end

"""
    get_initial_children(n::Network, var::Node)::Vector{Node}

Returns a list of the initial child nodes in the network from the network's
initial graph (i.e. the return graph from get_initial_graph(network)).
"""
function get_initial_children(n::Network, var::Node)::Vector{Node}
    [c for c in get_nodes(n) if var in get_initial_parents(n, c)]
end

"""
    get_transition_children(n::Network, var::Node)::Vector{Node}

Returns a list of transitioned child nodes in the network from the network's
transition graph (i.e. the return graph from get_transition_graph(network))
"""
function get_transition_children(n::Network, var::Node)::Vector{Node}
    [c for c in get_nodes(n) if var in get_transition_parents(n, c)]
end

"""
    complete_graph!(variables::Vector{<:Variable}, graph::VariableGraph)

Populate the `graph` with `v=>Node[]` for all `variables` not in the graph.
"""
function complete_graph!(variables::Vector{<:Variable}, graph::VariableGraph)
    ks = keys(graph)
    for v in variables
        if !(v in ks)
            graph[v] = Node[]
        end
    end
    return graph
end

"""
    struct InstantNetwork{I,J,O} <: Network{I,Nothing,O}

A network that only supports initial models.  The `Nothing` in the supertype
shows that there is no transition to another time.

See [`Network`](@ref) for the type parameters
"""
struct InstantNetwork{I,O} <: Network{I,Nothing,O}
    variables::Vector{<:Variable}
    placeholders::Vector{<:Placeholder}
    outputs::Vector{<:Variable} # outputs is a subset of variables
    parents::VariableGraph
    function InstantNetwork(variables::Vector{<:Variable}, parents::VariableGraph,
                placeholders::Vector{<:Placeholder} = Placeholder[], outputs::Vector{<:Variable} = Variable[])
        I = Tuple{[value_type(p) for p in placeholders]...}
        O = Tuple{[output_type(v) for v in outputs]...}
        return new{I,O}(variables, placeholders, outputs, complete_graph!(variables, parents))
    end
end

get_variables(n::InstantNetwork) = n.variables
get_placeholders(n::InstantNetwork) = n.placeholders
get_nodes(n::InstantNetwork) = union(Set(n.variables), Set(n.placeholders))
get_initial_placeholders(n::InstantNetwork) = n.placeholders
get_transition_placeholders(::InstantNetwork) = error("InstantNetwork does not have a transition graph")
get_outputs(n::InstantNetwork) = n.outputs
get_graph(n::InstantNetwork) = n.parents
get_initial_graph(n::InstantNetwork) = n.parents
get_transition_graph(::InstantNetwork) = error("InstantNetwork does not have a transition graph")
get_parents(n::InstantNetwork, v) = get_initial_parents(n,v)
get_children(n::InstantNetwork, v) = get_initial_children(n,v)

"""
    DynamicNetwork{I,J,O} <: Network{I,J,O}

A network that can transition over time.
See [`Network`](@ref) for the type parameters.
"""
struct DynamicNetwork{I,J,O} <: Network{I,J,O}
    variables::Vector{<:Variable}
    initial_placeholders::Vector{<:Placeholder}
    transition_placeholders::Vector{<:Placeholder}
    outputs::Vector{<:Variable} # outputs is a subset of variables
    initial_parents::VariableGraph
    transition_parents::VariableGraph
    function DynamicNetwork(variables::Vector{<:Variable},
                initial_parents::VariableGraph, transition_parents::VariableGraph,
                init_placeholders::Vector{<:Placeholder} = Placeholder[],
                trans_placeholders::Vector{<:Placeholder} = Placeholder[],
                outputs::Vector{<:Variable} = Variable[])
        I = Tuple{[value_type(p) for p in init_placeholders]...}
        J = Tuple{[value_type(p) for p in trans_placeholders]...}
        O = Tuple{[output_type(v) for v in outputs]...}
        return new{I,J,O}(variables, init_placeholders, trans_placeholders, outputs,
            complete_graph!(variables, initial_parents), complete_graph!(variables, transition_parents))
    end
end

get_variables(n::DynamicNetwork) = n.variables
get_initial_placeholders(n::DynamicNetwork) = n.initial_placeholders
get_transition_placeholders(n::DynamicNetwork) = n.transition_placeholders
get_nodes(n::DynamicNetwork) =
    union(Set(n.variables), Set(n.initial_placeholders), Set(n.transition_placeholders))
get_outputs(n::DynamicNetwork) = n.outputs
get_initial_graph(n::DynamicNetwork) = n.initial_parents
get_transition_graph(n::DynamicNetwork) = n.transition_parents

function assign_times(parents::Dict{Variable, Vector{Node}},
                      timed_varplaces::Dict{Node,Time})::Dict{Node,Time}
    # TODO be smarter. For now this is exhaustive at least, but exponentially slower than needed.
    # Unclear what the exact assumptions on model constraints should be
    # Easy improvement - select immediately identifiable vars with preference (i.e. FixedModel)

    untimed_vars = [var for var in parents.keys() if ~(var in keys(timed_varplaces))]
    if length(untimed_vars)==0
        return timed_varplaces
    end

    for untimed_var in untimed_vars
        new_time = assign_var(parents, timed_varplaces, untimed_var)
        if new_time
            sub_timed_varplaces = copy(timed_varplaces)
            sub_timed_varplaces[untimed_var] = new_time
            recursed = assign_times(parents, sub_timed_varplaces)
        else # Can't assign
            return false
        end

        if recursed
            return recursed
        end # Else try assigning another variable
    end

    return false
end


text:
using InterfaceTesting
using Base.Test
using Base.Test: Fail

type NoThrowTestSet <: Base.Test.AbstractTestSet
    results::Vector
    NoThrowTestSet(desc) = new([])
end
Base.Test.record(ts::NoThrowTestSet, t::Base.Test.Result) = (push!(ts.results, t); t)
Base.Test.finish(ts::NoThrowTestSet) = ts.results
num_fails(results) = sum(res->isa(res, Fail), results)

function includes_failing_string(results, substring)
    for res in results
        isa(res, Fail) || continue
        if contains(string(res), substring)
            return true
        end
    end
    false
end

immutable BrokenJuliaZeroFourIterator end
Base.start(::BrokenJuliaZeroFourIterator) = 1
Base.next(::BrokenJuliaZeroFourIterator, state) = state+1, state+1
Base.done(::BrokenJuliaZeroFourIterator, state) = state>5

immutable BrokenHasLengthIterator end
Base.start(::BrokenHasLengthIterator) = 1
Base.next(::BrokenHasLengthIterator, state) = state+1, state+1
Base.done(::BrokenHasLengthIterator, state) = state>5
Base.iteratorsize(::Type{BrokenHasLengthIterator}) = Base.HasLength #Return the Type not the instance
Base.length(::BrokenHasLengthIterator) = 5

immutable BrokenIteratorsizeIterator end
Base.start(::BrokenIteratorsizeIterator) = 1
Base.next(::BrokenIteratorsizeIterator, state) = state+1, state+1
Base.done(::BrokenIteratorsizeIterator, state) = false
Base.iteratorsize(::BrokenIteratorsizeIterator) = Base.IsInfinite()
#Set iteratorsize on the instance not the type

immutable BrokenLengthIterator end
Base.start(::BrokenLengthIterator) = 1
Base.next(::BrokenLengthIterator, state) = state+1, state+1
Base.done(::BrokenLengthIterator, state) = state>5
Base.iteratorsize(::Type{BrokenLengthIterator}) = Base.HasLength()
#Don't define the length

immutable UnknownLengthIterator end
Base.start(::UnknownLengthIterator) = 1
Base.done(::UnknownLengthIterator, state) = state>100
Base.iteratorsize(::Type{UnknownLengthIterator}) = Base.SizeUnknown()
function Base.next(::UnknownLengthIterator, state)
    ret = rand(1:10)
    ret, state+ret
end

@testset "partial_method_exists" begin
    @test partial_method_exists(length, (Vector,))
    @test !partial_method_exists(length, (Function,))

    @test partial_method_exists(rpad, (String,Int))
    @test partial_method_exists(rpad, (String,Int))
    @test !partial_method_exists(rpad, (Vector{Int},))
    @test partial_method_exists(rpad, (:, Int))
end

@testset "Iterators" begin
    test_iterator_interface(Vector)
    test_iterator_interface(Vector{Int})
    test_iterator_interface(typeof(zip([1,2,3],[2,23])))
    test_iterator_interface(typeof(zip([1,2,3],[2,23],[3,2,3])))
    test_iterator_interface(Matrix)
    test_iterator_interface(Int)

    test_iterator_interface(UnknownLengthIterator)

    BrokenIteratorsizeIterator_res = test_iterator_interface(BrokenIteratorsizeIterator; testset_type=NoThrowTestSet)
    @test 2==num_fails(BrokenIteratorsizeIterator_res)
    @test includes_failing_string(BrokenIteratorsizeIterator_res, "which(Base.iteratorsize")
    @test includes_failing_string(BrokenIteratorsizeIterator_res, "length")

    BrokenJuliaZeroFourIterator_res=test_iterator_interface(BrokenJuliaZeroFourIterator; testset_type=NoThrowTestSet)
    @test 1==num_fails(BrokenJuliaZeroFourIterator_res)
    @test includes_failing_string(BrokenJuliaZeroFourIterator_res, "length")

    BrokenLengthIterator_res=test_iterator_interface(BrokenLengthIterator; testset_type=NoThrowTestSet)
    @test 1==num_fails(BrokenLengthIterator_res)
    @test includes_failing_string(BrokenLengthIterator_res, "length")
end

@testset "Indexing" begin
    test_index_interface(Matrix)
    test_index_interface(SparseVector)
    test_index_interface(UnitRange{Int64}) #Logically this should *not* pass, since not setable, but that method as been coded to manually throw an error

    Dict_res = test_index_interface(Dict,  testset_type=NoThrowTestSet) #Not a proper indexable, no endof
    @test 0<num_fails(Dict_res)
end

@testset "Abstract Array" begin
    test_abstractarray_interface(Matrix)
    test_abstractarray_interface(SparseVector)
    test_abstractarray_interface(UnitRange{Int64}) #Logically this should *not* pass, since not setable, but that method as been coded to manually throw an error

    Dict_res = test_abstractarray_interface(Dict,  testset_type=NoThrowTestSet) #Not a proper indexable, no endof
    @test 0<num_fails(Dict_res)
end


text:
module Barcode

using DataFrames
using Images

include("code128_charset.jl")
include("encoding.jl")
include("pattern.jl")
include("utils.jl")

end


text:
using HTTP

base_url(year::Integer, day::Integer) = "https://adventofcode.com/$year/day/$day"

"""
    save_input(year, day, [session])

Retrieves input for a given day and year and saves it in e.g. "day2.in" for day 2.
Returns the filename of the stored input.
If no session is provided, the environment variable `AOC_SESSION` is used.
If the file already exists, no request will be made.
"""
function save_input(year, day, session::AbstractString)
    outfile = string("day", day, ".in")
    isfile(outfile) && return outfile

    res = HTTP.get(base_url(year, day) * "/input", cookies = Dict("session" => session))
    res.status == 200 ?  write(outfile, res.body) : error(res.body)
    outfile
end

save_input(year, day) = save_input(year, day, ENV["AOC_SESSION"])