text:
include("drawBlocks.jl")

function main()

    env = initCairo("out.svg",400.0,400.0)

    setLineWidthCairo(env ,1.0)

    #test a text box
    drawTextBox(env, "Text that needs \n to be split up to really \n this will a             lso just serve as a test",
                10.0,10.0)

    # draw diamond
    drawTextDiamond(env,"Text that needs to be split up",
                    100.0,200.0)

    destroyCairo(env)
end

main()

max_stars_repo_name: aklapatch/CodeVisualized

max_issues_repo_name: aklapatch/CodeVisualized

max_forks_repo_name: aklapatch/CodeVisualized


text:
module constants

export ME, RE, G, MM, RM, MCM, DISTANCE_TO_MOON, MOON_PERIOD, MOON_INITIAL_ANGLE, ORIGIN
export TOTAL_DURATION, MARKER_TIME, TOLERANCE, INITIAL_POSITION, INITIAL_VELOCITY

const ME = 5.97e24 # mass of earth in kg
const RE = 6.378e6 # radius of earth in m (at equator)
const G = 6.67e-11 # gravitational constant in m3 / kg s2
const MM = 7.35e22 # mass of moon in kg
const RM = 1.74e6 # radius of moon in m
const MCM = 5000. # mass of command module in kg
const DISTANCE_TO_MOON = 400.5e6 # m (actually, not at all a constant)
const MOON_PERIOD = 27.3 * 24.0 * 3600. # s
const MOON_INITIAL_ANGLE = pi / 180. * -61. # radians
const ORIGIN = [0., 0.] # Vector that represents the origin

const TOTAL_DURATION = 12. * 24. * 3600.# s
const MARKER_TIME = 0.5 * 3600. # (used in error correction) s
const TOLERANCE = 100000. # (used in error correction) m

const INITIAL_POSITION = [-6.701e6, 0.] # Initial vector for the spacecraft in m
const INITIAL_VELOCITY = [0., -10.818e3] # Initial velocity for spacecraft in m/s

end

ext: jl

max_stars_repo_name: rodrigolece/grupo-fisica-compu

max_issues_repo_name: rodrigolece/grupo-fisica-compu

max_forks_repo_name: rodrigolece/grupo-fisica-compu


text:
## path to TRPL folder --> set it before running the code!
path_to_trpl = "D:\\Projects\\TRPL" # some path to the TRPL folder (just an example here)

## activate package environment
cd(@__DIR__)
using Pkg
Pkg.activate("Project.toml")

## loading packages
using OrdinaryDiffEq, LinearAlgebra, StaticArrays, Random, NLopt, DelimitedFiles,
	  BlackBoxOptim, Printf, PoissonRandom, DataFrames, Plots, CSV

## backend quick plots
pyplot()
default(size = (800, 600))

## save output? (true or false)
out = false

## noise? (true or false)
noise = true

## experimental parameters (example for synthetic data set 1)
background = [1, 2, 12, 26, 18]
ND         = [1, 1,  2,  3,  4]
C_scaling  = 1e6

## parameters (example for synthetic data set 1)
# 		        τ_r0,  τ_nr0,   τ_d,  τ_0,   τ_1,  τ_l,     α,     β,    r,      C
parameters = [ 175.0,  350.0,  60.0, 20.0, 800.0, 250.0, 70.0, 110.0, 20.0, 1000.0]
eta0s      = [4.0e-5, 0.0004, 0.004, 0.04, 0.4]

## time (example for synthetic data set 1)
tspan = [0.0, 1750.0]
dt    = 0.25

## number of data set (example for synthetic data set 1)
idx_syn = 1

## rng seeds (to get the same noise when recalculating a data set: rng_seed = abs(rand(Int)))
rng_seeds = [4057480306083410532, 7737535138175179370, 3066106408424540935,
			 3829281771415197884, 1025260804200922037]

## useful variables
num_pars = 10
num_traj = length(eta0s)
dims =  num_pars + num_traj

## output name
name_syn_data = "syn_data_$(idx_syn)"

## definition rate equations
function rate_equations(du, u, p, t)
	du[1] = -((1 - u[2]) / p[3]) * u[1] - ((1 - u[3]) / p[4]) * u[1] + (1 / p[5]) * u[3] * 1 / p[7] -
			((1 + u[1]) / p[1]) * u[1] - ((1 + u[1]) / (p[2] * (1 + p[9] * u[1]))) * u[1]
	du[2] = ((1 - u[2]) / (p[3])) * p[8] * u[1]
	du[3] = ((1 - u[3]) / p[4]) * p[7] * u[1] - u[3] / p[5] - u[3] / p[6]
end

## ODEProblem definition with temporary u0 and parameters
u0 = [0.1, 0, 0]
t = [tspan[1]:dt:tspan[2];]
p_true = vcat(parameters,eta0s)
prob = ODEProblem(rate_equations, u0, tspan, p_true)

## function to generate the model for plotting
function generateipl(par)
    function prob_func(prob, i, repeat)
        remake(prob, u0 = [par[num_pars+i], 0, 0])
    end
    ensemble_simu = EnsembleProblem(prob, prob_func = prob_func)
    sol = Array(solve(ensemble_simu, Rosenbrock23(), p = par[1:num_pars-1], EnsembleThreads(),
		     		  trajectories = num_traj, saveat = t, save_idxs = [1],
		     		  reltol = 1e-6, abstol = 1e-8))[1,:,:]
    ipl = zeros(size(sol))
    for i = 1:num_traj
    	ipl[:,i] = ((C_scaling * par[num_pars] / (10 ^ ND[i])) .* sol[:,i] .* (sol[:,i] .+ 1)) .+ background[i]
    end
    return ipl
end

## get synthetic data with poisson noise if wanted
IPL = generateipl(p_true)
if noise
	IPL = pois_rand.(Random.seed!(rng_seeds[idx_syn]),IPL)
end

## get non-zero indicies
idx_IPL = zeros(Bool,size(IPL))
for i = 1:num_traj
	idx_IPL[:,i] = IPL[:,i] .> 0
end

## setup loss function to calculate loss goal
function loss(par)
	l = 0.0
	function prob_func(prob, i, repeat)
		remake(prob, u0 = [par[num_pars+i], 0, 0])
	end
	ensemble_simu = EnsembleProblem(prob, prob_func = prob_func)
	sol = Array(solve(ensemble_simu, Rosenbrock23(), p = par[1:num_pars-1], EnsembleThreads(),
					  trajectories = num_traj, saveat = t,save_idxs = [1],
					  reltol=1e-6,abstol=1e-8))[1,:,:]
	for i = 1:num_traj
		model = ((C_scaling * par[num_pars] / (10 ^ ND[i])) .* sol[:,i] .* (sol[:,i] .+ 1)) .+ background[i]
		idx_used = idx_IPL[:,i] .& (model .> 0)
		l += sum((model[idx_used]) .- (IPL[idx_used,i] .* log.(model[idx_used])))
	end
	return l
end

## calculate and print loss goal
l = loss(p_true)
println("loss goal = $l")

## plotting font size
fs = 16

## plot data
ylim = (10^0,10^(ceil(log10(maximum(IPL)))))
plt = plot(palette=palette(:seaborn_colorblind, num_traj),xlims=(t[1],t[end]),yscale=:log10,
		   ylims=ylim,title="Synthetic data set $idx_syn",titlefontsize=fs,grid=true,
		   xguidefontsize=fs,yguidefontsize=fs,dpi=300,framestyle=:box,xtickfontsize=fs,
		   ytickfontsize=fs,xlabe="t/ns",ylabel="Photon Counts",legendfontsize=fs-4,
		   leg=:topright,margin=10Plots.mm)
for i = num_traj:-1:1
	plot!(t, IPL[:,i],linealpha=0.75,linewidth=1.5,label=@sprintf("η₀ = %.2E", p_true[num_pars+i]))
end
display(plt)

## save output
if out
	# save input parameters
	io = open(joinpath(path_to_trpl,"Data\\Syn\\Params\\$name_syn_data.txt"), "w")
		println(io, "parameters")
		println(io, "τ_r0, τ_nr0, τ_2t, τ_1t,  τ_1e, τ_1d, α, β, r, C, η0s")
		println(io, p_true)
		println(io, "C scaling")
		println(io, C_scaling)
		println(io, "background")
		println(io, background)
		println(io, "ND")
		println(io, ND)
		println(io, "time span")
		println(io, tspan)
		println(io, "dt")
		println(io, dt)
		println(io, "loss goal")
		println(io, l)
		if noise
			println(io, "rng seed")
			println(io, rng_seed)
		end
	close(io)

	# save data for fit
	output = Array{Any}(undef,(size(IPL,1)+3,num_traj+1))
	output[1,2] = C_scaling
	output[2,2:end] = ND
	output[3,2:end] = background
	output[4:end,1] = t
	output[4:end,2:end] = IPL
	output[1:3,1] .= 0
	output[1,3:end] .= 0
	df = DataFrame(output,:auto)
	CSV.write(joinpath(path_to_trpl,"Data\\Syn\\CSV\\$name_syn_data.csv"), df)
end

max_stars_repo_name: MaxGrossmann/TRPL

max_issues_repo_name: MaxGrossmann/TRPL

max_forks_repo_name: MaxGrossmann/TRPL


text:
""" Tonti Diagrams
This file includes a framework for constructing and evaluating Tonti diagrams.
Tonti diagrams are stored as ACSets, and have an imperative interface for
describing physical variables and the relationships between them. This tooling
also lets a Tonti diagram be converted to a vectorfield, allowing for
simulation of physical systems through any traditional vectorfield solver.
"""

module Tonti
using Catlab.Present
using Catlab.Theories
using Catlab.CategoricalAlgebra
using Catlab.CategoricalAlgebra.FinSets
using Catlab.Programs
using Catlab.WiringDiagrams
using CombinatorialSpaces
import CombinatorialSpaces.DualSimplicialSets: EmbeddedDeltaDualComplex2D, ⋆

export TontiDiagram, AbstractTontiDiagram, addTransform!, addSpace!,
       addTime!, vectorfield, disj_union, gen_form, addBC!

# Helper Function for constructing dual complexes
function EmbeddedDeltaDualComplex2D(d::EmbeddedDeltaSet2D{O, D}) where {O, D}
  sd = EmbeddedDeltaDualComplex2D{O, eltype(D), D}(d)
  subdivide_duals!(sd, Barycenter())
  sd
end;

""" Construct a 0-form based on a scalar function

This operator accepts a scalar function and evaulates it at each point on the
simplex, returning a 0-form.
"""
function gen_form(s::EmbeddedDeltaSet2D, f::Function)
  map(f, point(s))
end

# This defines the structure of the Tonti Diagram ACSet
@present TheoryTontiDiagram(FreeSchema) begin
  Func::Data
  Label::Data

  Variable::Ob
  Corner::Ob
  Transform::Ob
  InParam::Ob
  OutParam::Ob

  tgt::Hom(InParam, Transform)
  src::Hom(OutParam, Transform)
  in_var::Hom(InParam, Variable)
  out_var::Hom(OutParam, Variable)

  func::Attr(Transform, Func)
  t_type::Attr(Transform, Label)

  corner::Hom(Variable, Corner)

  c_label::Attr(Corner, Label)
  v_label::Attr(Variable, Label)
end

# These functions are necessary for defining a TontiDiagram data structure from
# the ACSet
const AbstractTontiDiagram = AbstractACSetType(TheoryTontiDiagram)
const TontiDiagram = ACSetType(TheoryTontiDiagram,index=[:src, :tgt, :corner,
                                                         :in_var, :out_var],
                                                 unique_index=[:v_label,
                                                               :c_label])

# Define an interface for an OpenTontiDiagram which allows Tonti diagrams to
# be composed over their corners
const OpenTontiDiagramOb, OpenTontiDiagram = OpenACSetTypes(TontiDiagram,
                                                            :Corner)
Open(td::AbstractTontiDiagram) = begin
  OpenTontiDiagram{Function, Symbol}(td,
                                     FinFunction(collect(1:nparts(td, :Corner)),
                                                 nparts(td, :Corner)))
end

""" Constructor for Tonti Diagram

This constructor generates an empty Tonti diagram of a given dimension. The
dimensionality determines which spatial elements will be included in the
data structure.
"""
function TontiDiagram(dimension = 3)
  td = TontiDiagram{Function, Symbol}()
  corners = ['P', 'L', 'S', 'V']

  for d in 1:(dimension+1)
    d2 = dimension - (d-1) + 1
    primal = add_parts!(td, :Corner, 2, c_label=[Symbol("I$(corners[d])"),
                                                  Symbol("T$(corners[d])")])
    dual = add_parts!(td, :Corner, 2, c_label=[Symbol("I$(corners[d2])2"),
                                                Symbol("T$(corners[d2])2")])
  end
  td
end

""" Constructor for Tonti Diagram with variables

This constructor accepts variables which are each associated with a corner.
"""
function TontiDiagram(dimension::Int, variables::Array{Pair{Symbol, Symbol}})
  td = TontiDiagram(dimension)
  add_parts!(td, :Variable, length(variables), v_label=first.(variables),
             corner=[first(incident(td, s, :c_label)) for s in last.(variables)])
  td
end

# Helper functions for working on Tonti diagrams
dimension(td::AbstractTontiDiagram) = nparts(td, :Corner) ÷ 4 - 1
var_corner(td::AbstractTontiDiagram, label::Symbol) = td[incident(td, label, :v_label), :corner]
get_var_ind(td::AbstractTontiDiagram, label::Symbol) = first(incident(td, label, :v_label))

""" Generates the appropriate hodge star operator between given corners

This operator will produce a scalar identity if the operator is between corners
where a hodge star is no defined.

TODO: This should throw an error instead of failing silently for invalid
inputs. Currently fails silently because the hodge star is added between all
transforms (this needs to be fixed in vectorfield before here)
"""
function gen_star(dom::Int64, codom::Int64, td::AbstractTontiDiagram, dual::EmbeddedDeltaDualComplex2D)
  corner2ind = Dict('P'=>0, 'L'=>1, 'S'=>2, 'V'=>3)
  # Check for which ⋆ operator to apply
  if "$(td[dom, :c_label])"[end] == '2'
    if "$(td[codom, :c_label])"[end] == '2'
      return 1
      #error("Transformation $(td[dom, :c_label]) -> $(td[codom, :c_label]) is not between primal and dual complexes")
    end
    inv(⋆(corner2ind["$(td[codom, :c_label])"[2]], dual))
  else
    if "$(td[codom, :c_label])"[end] != '2'
      return 1
      #error("Transformation $(td[dom, :c_label]) -> $(td[codom, :c_label]) is not between primal and dual complexes")
    end
    ⋆(corner2ind["$(td[dom, :c_label])"[2]], dual)
  end
end

""" Selects appropriate hodge star operator between given corners from a given
pre-computed array

TODO: This should also throw an error, but currently fails silently. See other
gen_star for details.
"""
function gen_star(dom::Int64, codom::Int64, td::AbstractTontiDiagram, star_arr)
  corner2ind = Dict('P'=>0, 'L'=>1, 'S'=>2, 'V'=>3)
  # Check for which ⋆ operator to apply
  if "$(td[dom, :c_label])"[end] == '2'
    if "$(td[codom, :c_label])"[end] == '2'
     return 1
    end
    inv(star_arr[corner2ind["$(td[codom, :c_label])"[2]]+1])
  else
    if "$(td[codom, :c_label])"[end] != '2'
      return 1
    end
    star_arr[corner2ind["$(td[dom, :c_label])"[2]]+1]
  end
end

""" Adds a transformation to the Tonti diagram between given variables

TODO: Remove this function, since hodge-star is not defined element-wise and
      cannot be added at this step.
"""
function addTransform!(td::AbstractTontiDiagram, complex::EmbeddedDeltaSet2D,
                       dom::Array{Symbol}, func::Function,
                       codom::Array{Symbol})
  # Check rules for transforms
  dual = EmbeddedDeltaDualComplex2D(complex)
  dom_sim = var_corner(td, first(dom))
  codom_sim = var_corner(td, first(codom))
  all(var_corner(td, d) == dom_sim for d in dom) || error("Domain is not consistently on the corner $dom_sim")
  all(var_corner(td, c) == codom_sim for c in codom) || error("Codomain is not consistently on the corner $codom_sim")

  star_op = gen_star(dom_sim, codom_sim, td, dual)
  new_func(x...) = begin
    func(x...)
  end
  # Add transformation
  tran = add_part!(td, :Transform, func=new_func, t_type=:constitutive)
  add_parts!(td, :InParam, length(dom), tgt=tran, in_var=[get_var_ind(td, v) for v in dom])
  add_parts!(td, :OutParam, length(codom), src=tran, out_var=[get_var_ind(td, v) for v in codom])
end

""" Adds a transformation to the Tonti diagram between given variables

This expects a Julia function which will take as input an array of values the
length of argument `dom` and will return an array the length of argument
`codom`. The order of variables passed to and recieved from `func` is
determined by the order in `dom` and `codom`.
"""
function addTransform!(td::AbstractTontiDiagram, dom::Array{Symbol}, func::Function, codom::Array{Symbol})
  # Check rules for transforms
  dom_sim = var_corner(td, first(dom))
  codom_sim = var_corner(td, first(codom))
  all(var_corner(td, d) == dom_sim for d in dom) || error("Domain is not consistently on the corner $dom_sim")
  all(var_corner(td, c) == codom_sim for c in codom) || error("Codomain is not consistently on the corner $codom_sim")

  # Add transformation
  tran = add_part!(td, :Transform, func=func, t_type=:constitutive)
  add_parts!(td, :InParam, length(dom), tgt=tran, in_var=[get_var_ind(td, v) for v in dom])
  add_parts!(td, :OutParam, length(codom), src=tran, out_var=[get_var_ind(td, v) for v in codom])
end

""" Adds a topological transformation to the Tonti diagram between given variables

The primary purpose of this function is to ensure that the variables are
appropriate for a topological transformation to exist between them.
"""
function addTopoTransform!(td::AbstractTontiDiagram, dom::Symbol, func::Function, codom::Symbol)
  v_dom = incident(td, dom, [:corner, :c_label])
  v_codom = incident(td, codom, [:corner, :c_label])

  (length(v_dom) <= 1 && length(v_codom) <= 1) ||
      error("Adding topological transformations with >1 variable per corner is not yet supported")

  if length(v_dom) == 0 || length(v_codom) == 0
    return false
  end

  tran = add_part!(td, :Transform, func=func, t_type=:topological)
  add_part!(td, :InParam, tgt=tran, in_var=v_dom[1])
  add_part!(td, :OutParam, src=tran, out_var=v_codom[1])
  return true
end

""" Adds a temporal transformation to the Tonti diagram between given variables

The primary purpose of this function is to ensure that the variables are
appropriate for a temporal transformation to exist between them.
"""
function addTempTransform!(td::AbstractTontiDiagram, dom::Symbol, func::Function, codom::Symbol)
  v_dom = incident(td, dom, [:corner, :c_label])
  v_codom = incident(td, codom, [:corner, :c_label])

  (length(v_dom) <= 1 && length(v_codom) <= 1) ||
      error("Adding topological transformations with >1 variable per corner is not yet supported")

  if length(v_dom) == 0 || length(v_codom) == 0
    return false
  end

  tran = add_part!(td, :Transform, func=func, t_type=:temporal)
  add_part!(td, :InParam, tgt=tran, in_var=v_dom[1])
  add_part!(td, :OutParam, src=tran, out_var=v_codom[1])
  return true
end

""" Adds a boundary condition to a Tonti Diagram

TODO: Extend this to beyond a mask over certain values. Currently this allows
for constant boundaries by applying a mask to the vectorfield result, but this
does not generalize to all mathematical boundary conditions.
"""
function addBC!(td::AbstractTontiDiagram, dom::Symbol, mask::Function)
  v_dom = incident(td, dom, [:v_label])

  tran = add_part!(td, :Transform, func=mask, t_type=:bc)
  add_part!(td, :InParam, tgt=tran, in_var=v_dom[1])
  add_part!(td, :OutParam, src=tran, out_var=v_dom[1])
end

""" Adds topological transformations to a 1D Tonti Diagram

Note that since the complex is not embedded, we are unable to use a true
dual derivative.
"""
function addSpace!(td::AbstractTontiDiagram, complex::AbstractDeltaSet1D)
  bound_1_0   = boundary(1,complex)
  cobound_0_1 = d(0,complex)
  addTopoTransform!(td, :IP, x->(cobound_0_1*x), :IL)
  addTopoTransform!(td, :TP, x->(cobound_0_1*x), :TL)

  # TODO: Add Hodge * operator instead of just swapping bound/cobound
  addTopoTransform!(td, :IP2, x->(bound_1_0*x), :IL2)
  addTopoTransform!(td, :TP2, x->(bound_1_0*x), :TL2)

  td
end

""" Adds topological transformations to a 2D Tonti Diagram
(from an embedded complex)
"""
function addSpace!(td::AbstractTontiDiagram, complex::EmbeddedDeltaSet2D)
  dual = EmbeddedDeltaDualComplex2D(complex)
  d_0_1 = d(0,complex)
  d_1_2 = d(1,complex)

  dd_0_1 = dual_derivative(0, dual)
  dd_1_2 = dual_derivative(1, dual)

  addTopoTransform!(td, :IP, x->(d_0_1*x), :IL)
  addTopoTransform!(td, :TP, x->(d_0_1*x), :TL)
  addTopoTransform!(td, :IL, x->(d_1_2*x), :IS)
  addTopoTransform!(td, :TL, x->(d_1_2*x), :TS)

  # TODO: Add Hodge * operator instead of just swapping bound/cobound
  addTopoTransform!(td, :IP2, x->(dd_0_1*x), :IL2)
  addTopoTransform!(td, :TP2, x->(dd_0_1*x), :TL2)
  addTopoTransform!(td, :IL2, x->(dd_1_2*x), :IS2)
  addTopoTransform!(td, :TL2, x->(dd_1_2*x), :TS2)
  td
end

""" Adds topological transformations to a 2D Tonti Diagram

Note that since the complex is not embedded, we are unable to use a true
dual derivative.

TODO: These should be able to be removed, but considerations must be made
for what it means for a particle system to exist on an embedded complex.
"""
function addSpace!(td::AbstractTontiDiagram, complex::AbstractDeltaSet2D)
  d_0_1 = d(0,complex)
  d_1_2 = d(1,complex)

  dd_0_1 = boundary(2, complex)
  dd_1_2 = boundary(1, complex)

  addTopoTransform!(td, :IP, x->(d_0_1*x), :IL)
  addTopoTransform!(td, :TP, x->(d_0_1*x), :TL)
  addTopoTransform!(td, :IL, x->(d_1_2*x), :IS)
  addTopoTransform!(td, :TL, x->(d_1_2*x), :TS)

  # TODO: Add Hodge * operator instead of just swapping bound/cobound
  addTopoTransform!(td, :IP2, x->(dd_0_1*x), :IL2)
  addTopoTransform!(td, :TP2, x->(dd_0_1*x), :TL2)
  addTopoTransform!(td, :IL2, x->(dd_1_2*x), :IS2)
  addTopoTransform!(td, :TL2, x->(dd_1_2*x), :TS2)

  td
end

""" Adds temporal transformations to a Tonti Diagram

Currently all temporal transformations simply apply Euler's method. When
attached to DifferentialEquations, more complex timestepping can be
accomplished.

TODO: Develop more time-stepping techniques for the Tonti diagram structure
"""
function addTime!(td::AbstractTontiDiagram; dt=1)
  addTempTransform!(td, :TP, x->dt*x, :IP)
  addTempTransform!(td, :TL, x->dt*x, :IL)
  addTempTransform!(td, :TP2, x->dt*x, :IP2)
  addTempTransform!(td, :TL2, x->dt*x, :IL2)

  if dimension(td) > 1
    addTempTransform!(td, :TS, x->dt*x, :IS)
    addTempTransform!(td, :TS2, x->dt*x, :IS2)
  end

  if dimension(td) > 2
    addTempTransform!(td, :TV, x->dt*x, :IV)
    addTempTransform!(td, :TV2, x->dt*x, :IV2)
  end
end

""" Generates a vectorfield for simulation of a system given a Tonti diagram

This function performs a topological sort to determine dependencies between
variables (forgetting the temporal transformations). It then determines the
necessary set of state variables to define the value of each variable for each
time step.

TODO: Convert dependencies to Graph in Catlab and take advantage of the
topological sort tooling from there. Will significantly reduce the size of this
function.
"""
function vectorfield(td::AbstractTontiDiagram, complex::Union{AbstractDeltaSet1D,
                                                              AbstractDeltaSet2D})
  # Initialize hodge-stars if the tonti diagram is embedded
  star_arr = Array{Union{Number,AbstractArray},1}(undef, 3)
  if(has_part(complex, :Point))
    dual = EmbeddedDeltaDualComplex2D(complex)
    star_arr .= [⋆(i, dual) for i in 0:2]
  else
    star_arr .= ones(Int64,3)
  end

  # Define order of evaluation for transformations
  var_deps = [Set(filter(t -> !(td[td[t,:src], :t_type] in [:temporal, :bc]), incident(td, i, :out_var)))
              for i in 1:nparts(td, :Variable)] # deps per variable

  tran_deps = [Set(td[incident(td, t, :src), :in_var])
               for t in 1:nparts(td, :Transform)] # deps per transforms

  var_to_bc = Dict(map(filter(t->(td[t, :t_type] == :bc), collect(1:nparts(td, :Transform)))) do t
                                td[incident(td, t,:src)[1], :out_var] => t
                              end
                             ) # Map from variables to their boundary conditions (if any)
  evaluated = [td[t, :t_type] == :bc for t in 1:nparts(td, :Transform)]
  to_evaluate = Array{Int64, 1}()
  to_eval_vars = Array{Int64, 1}()
  order = Array{Int64, 1}()
  time_vars = td[incident(td, :temporal, [:src, :t_type]), :out_var]

  append!(to_eval_vars, time_vars)
  for v in to_eval_vars
    if v in keys(var_to_bc)
      push!(order, var_to_bc[v])
    end
  end

  # Add any transformations which now have all variables satisfied
  for v in to_eval_vars
    for t in 1:length(tran_deps)
      delete!(tran_deps[t], v)
      if isempty(tran_deps[t]) && !evaluated[t]
        push!(to_evaluate, t)
        evaluated[t] = true
      end
    end
  end
  to_eval_vars = empty(to_eval_vars)

  # Evaluate transforms which have no more unsatisfied variable dependencies
  while !isempty(to_evaluate)
    cur_tran = pop!(to_evaluate)
    push!(order, cur_tran)

    outputs = incident(td, cur_tran, :src)
    out_vars = td[outputs, :out_var]

    # Update variable dependencies based on this transformation's eval
    for out in 1:length(outputs)
      delete!(var_deps[out_vars[out]], outputs[out])
      if isempty(var_deps[out_vars[out]])
        push!(to_eval_vars, out_vars[out])
        if out_vars[out] in keys(var_to_bc)
          push!(order, var_to_bc[out_vars[out]])
        end
      end
    end

    # Add transforms which have all variables satisfied
    for v in to_eval_vars
      for t in 1:length(tran_deps)
        delete!(tran_deps[t], v)
        if isempty(tran_deps[t]) && !evaluated[t]
          push!(to_evaluate, t)
          evaluated[t] = true
        end
      end
    end
    to_eval_vars = empty(to_eval_vars)
  end

  # Initialize memory
  data, v2ind = initData(td, complex)

  # Create maps from tonti diagram variables to the initialized memory
  timevar_to_ind = Dict{Int, Tuple{Int,Int}}()
  for v in 1:length(time_vars)
    cur_corner = td[time_vars[v], :corner]
    data[cur_corner] = vcat(data[cur_corner], zeros(1,size(data[cur_corner])[2]))
    timevar_to_ind[time_vars[v]] = (cur_corner,size(data[cur_corner])[1])
  end

  transforms = Array{Pair{Array{Tuple{Int, Int},1},Array{Tuple{Int, Int},1}},1}()
  masks = Dict{Pair{Int, Int}, Array{Number}}()
  star_op = Array{Union{Number,AbstractArray},1}()

  for t in 1:nparts(td, :Transform)
    in_vars = [(x ∈ time_vars) ? timevar_to_ind[x] : v2ind[x]
               for x in td[incident(td, t, :tgt), :in_var]]
    out_vars = v2ind[td[incident(td, t, :src), :out_var]]
    push!(star_op, gen_star(first(first(in_vars)), first(first(out_vars)), td, star_arr))
    push!(transforms, in_vars => out_vars)
  end

  # Vectorfield function definition
  function system(du, u, t, p)
    # Reset data to 0
    for i in 1:length(data)
      data[i] .= 0
    end

    # Update current state from u
    cur_ind = 1
    for v in time_vars
      c_ind = first(timevar_to_ind[v])
      v_ind = last(timevar_to_ind[v])
      v_len = size(data[c_ind])[2]
      data[c_ind][v_ind,:] .= view(u, cur_ind:(cur_ind+v_len-1))
      cur_ind += v_len
    end

    for t in order
      cur_t = transforms[t]
      dom_c = first(first(first(cur_t)))
      codom_c = first(first(last(cur_t)))
      cur_f = subpart(td, t, :func)
      if subpart(td, t, :t_type) == :topological
        # Topological transforms are a special case (and always 1:1)
        # Will likely want to deal with these uniquely in future implementations
        data[codom_c][last.(cur_t[2]), :] .+= cur_f(data[dom_c][last.(cur_t[1]),:]')'
      elseif subpart(td, t, :t_type) == :bc
        data[codom_c][last.(cur_t[2]), :] .*= cur_f()'
      else
        data[codom_c][last.(cur_t[2]), :] .+= mapslices(subpart(td, t, :func), view(data[dom_c], last.(cur_t[1]), :), dims=1) * star_op[t]
      end
    end
    du .= vcat([data[v2ind[v][1]][v2ind[v][2], :][:] for v in time_vars]...)
  end
  system, [td[v, :v_label] => size(data[v2ind[v][1]])[2] for v in time_vars]
end

# Helper function which intializes the data storage for vectorfield for a 1D
# system
function initData(td::AbstractTontiDiagram, complex::AbstractDeltaSet1D)
  data = Array{Array{Float64, 2}, 1}()
  corner_to_len = Dict(:IP => nv(complex), :TP => nv(complex), :IL2 => nv(complex), :TL2 => nv(complex),
                       :IL => ne(complex), :TL => ne(complex), :IP2 => ne(complex), :TP2 => ne(complex))

  v2ind = fill((0,0), nparts(td, :Variable))
  for c in 1:nparts(td, :Corner)
    cur_vars = incident(td, c, :corner)
    push!(data, zeros(length(cur_vars),corner_to_len[td[c, :c_label]]))
    v2ind[cur_vars] = [(c, v) for v in 1:length(cur_vars)]
  end
  data, v2ind
end

# Helper function which intializes the data storage for vectorfield for a 2D
# system
function initData(td::AbstractTontiDiagram, complex::AbstractDeltaSet2D)
  data = Array{Array{Float64, 2}, 1}()
  c_nv = nv(complex)
  c_ne = ne(complex)
  c_nt = ntriangles(complex)
  corner_to_len = Dict(:IP => c_nv, :TP => c_nv, :IS2 => c_nv, :TS2 => c_nv,
                       :IL => c_ne, :TL => c_ne, :IL2 => c_ne, :TL2 => c_ne,
                       :IS => c_nt, :TS => c_nt, :IP2 => c_nt, :TP2 => c_nt)

  v2ind = fill((0,0), nparts(td, :Variable))
  for c in 1:nparts(td, :Corner)
    cur_vars = incident(td, c, :corner)
    push!(data, zeros(length(cur_vars),corner_to_len[td[c, :c_label]]))
    v2ind[cur_vars] = [(c, v) for v in 1:length(cur_vars)]
  end
  data, v2ind
end

""" Constructs a Tonti diagram from two source Tonti diagrams, combining over
the corners.

The result Tonti diagram contains all information from the independent Tonti
diagrams and with variables on the same corners as before.

TODO: Make this accept a UWD so that tonti diagrams can be joined over shared
variables as well (also just provides a cleaner joining interface)
"""
function disj_union(td1::AbstractTontiDiagram, td2::AbstractTontiDiagram)
  o_td1 = Open(td1)
  o_td2 = Open(td2)

  td_merge = @relation (x,) begin
     td1(x)
     td2(x)
  end

  apex(oapply(td_merge, Dict(:td1=>o_td1, :td2=>o_td2)))
end

end


text:
# Exports
export Orbit
export interpolate
export Spacecraft
export Location
export Request
export load_requests
export set_collect_duration
export GroundStation
export load_stations
export Opportunity
export Done
export Noop
export Sunpoint
export Collect
export Contact
export SatPlanningProblem
export add_locations
export clear_locations
export clear_opportunities
export clear_all

#########
# Orbit #
#########

"""
Orbit provides a data structure which wraps SatelliteDynamics propagators

Arguments:
- `epc0::Union{Epoch, Real}` Initial Epoch of propagation
- `epcf::Union{Epoch, Real}` Final Epoch of propagation
- `eci0::Array{<:Real, 1}` Initial inertial state at initial Epoch [m; m/s]
- `id::Integer` Identification number for this orbit
- `timestep::Real` Timestep used for propgation output
- `dtmax::Real` Maximum allowable integrator step
- `mass::Real` Satellite mass [kg]
- `area_drag::Real` Wetted wind-facing area [m^2]
- `coef_drag::Real` Coefficient of drag
- `area_srp::Real` Area illuminated by solar radiation pressure [m^2]
- `coef_srp::Real` Coefficient of reflectivity
- `n_grav::Integer` Gravity model degree
- `m_grav::Integer` Gravity model order
- `drag::Bool` Enable drag perturbation in propgation
- `srp::Bool` Enable solar radiation pressure perturbation in propgation
- `moon::Bool` Enable thirdbody moon perturbation in propgation
- `sun::Bool` Enable thirdbody sun perturbation in propgation
- `relativity::Bool` Enable relativity perturbation in propgation

Attributes:
- `id::Int32` Identification number for the specific satellite
- `t::Array{Float64, 1}` Elapsed time of propagator output
- `epc::Array{Epoch, 1}` Epoch of propagator output
- `eci::Array{Float64, 2}` Earth Centered Inertial state of satellite output. State entries alightex with column output.
- `ecef::Array{Float64, 2}` Earth Centered Earth Fixed state of statellite output. State entries alightex with column output.
- `oe::Array{Float64, 2}` Osculating orbital element state of statellite output. State entries alightex with column output.
"""
mutable struct Orbit
    id::Int32
    t::Array{Float64, 1}
    epc::Array{Epoch, 1}
    eci::Array{Float64, 2}
    ecef::Array{Float64, 2}
    oe::Array{Float64, 2}

    function Orbit(epc0::Union{Epoch, Real},  epcf::Union{Epoch, Real}, eci0::Array{<:Real, 1};
                id=0::Integer, timestep=1.0::Real, dtmax=30.0::Real,
                mass=100.0::Real, area_drag=1.0::Real, coef_drag=2.3::Real,
                area_srp=1.0::Real, coef_srp=1.8::Real,
                n_grav=20::Integer, m_grav=20::Integer,
                drag=true::Bool, srp=true::Bool, moon=true::Bool,
                sun=true::Bool, relativity=true::Bool)

        # Initialize State Vector
        orb  = EarthInertialState(epc0, eci0, dt=timestep,
            mass=mass,
            area_drag=area_drag, coef_drag=coef_drag,
            area_srp=area_srp, coef_srp=coef_srp,
            n_grav=n_grav, m_grav=m_grav,
            drag=drag, srp=srp,
            moon=moon, sun=sun,
            relativity=relativity
        )

        # Propagate the orbit
        t, epc, eci = sim!(orb, epcf)

        ecef = zeros(Float64, 6, length(epc))
        oe   = zeros(Float64, 6, length(epc))

        for i in 1:length(t)
            ecef[:, i] = sECItoECEF(epc[i], eci[:, i])
            oe[:, i]   = sCARTtoOSC(eci[:, i], use_degrees=true)
        end

        new(id, t, epc, eci, ecef, oe)
    end
end

function Base.show(io::IO, orb::Orbit)

    s = @sprintf "Orbit(%d)" orb.id

    print(io, s)
end

# Helpfer function
function searchsortednearest(a,x)
    idx = searchsortedfirst(a,x)
    if (idx==1); return idx; end
    if (idx>length(a)); return length(a); end
    if (a[idx]==x); return idx; end
    if (abs(a[idx]-x) < abs(a[idx-1]-x))
       return idx
    else
       return idx-1
    end
 end

"""
Interpolate orbit information to a given `Epoch`.

This interpolation is not particularly accurate, care should be used to determine
that the use-case can handle the interprolation errors introduced.

Arguments:
- `orbit::Orbit` Orbit data to interpolate
- `epc::Union{Epoch, Real}` Epoch to interpolate state output to.

Returns:
- `eci::Arrray{Float64, 1}` Earth intertial state information interpolated to time of Epoch
"""
function interpolate(orbit::Orbit, epc::Union{Epoch, Real})
    # Check validity of input
    if epc < orbit.epc[1]
        throw(ArgumentError("Invaid interpolation time. Cannot interpolate before start of Orbit."))
    end

    if epc > orbit.epc[end]
        throw(ArgumentError("Invaid interpolation time. Cannot interpolate after end of Orbit."))
    end

    # Exit early if epoch is end of interpolation (corner case)
    if epc == orbit.epc[end]
        return orbit.eci[:, end]
    end

    # Get Uppler and lower time bounds
    idx_l = searchsortedfirst(orbit.epc, epc)
    idx_u = idx_l + 1

    # Epoch
    epc_l = orbit.epc[idx_l]
    epc_u = orbit.epc[idx_u]

    # Earth Intertial State
    eci_l = orbit.eci[:, idx_l]
    eci_u = orbit.eci[:, idx_u]

    # Interpolate
    eci = zeros(typeof(eci_l[1]), 6)
    for i in 1:length(eci)
        eci[i] = (eci_u[i] - eci_l[i])/(epc_u - epc_l)*(epc - epc_l) + eci_l[i]
    end

    return eci
end

##############
# Spacecraft #
##############

"""
Object to store spacecraft information

Model Parameters:
- `powergen_sunpoint::Float64` Power generated when sunpointed
- `powergen_image::Float64` Power generated when imaging (should be negative)
- `powergen_downlink::Float64` Power generated when downlinking (should be negative)
- `datagen_backorbit::Float64` Data generated in backorbit
- `datagen_image::Float64` Data generated when imaging
- `datagen_downlink::Float64` Data generated when downlinking (should be negative)
- `slew_rate::Float64` Spacecraft slew rate (degrees per second)
"""
@with_kw mutable struct Spacecraft
    # Core Parameters
    id::Integer
    tle::TLE

    # Resource Model
    powergen_sunpoint::Float64 = 0.0
    powergen_image::Float64 = 0.0
    powergen_downlink::Float64 = 0.0
    datagen_backorbit::Float64 = 0.0
    datagen_image::Float64 = 0.0
    datagen_downlink::Float64 = 0.0
    slew_rate::Float64 = 1.0
end

############
# Location #
############

abstract type Location end

function Base.isequal(ll::Location, rl::Location)
    return (
        (ll.id == rl.id) &&
        (ll.lat == rl.lat) &&
        (ll.lon == rl.lon) &&
        (ll.ecef == rl.ecef) && true
    )
end
Base.:(==)(ll::Location, rl::Location) = Base.isequal(ll, rl)

function Base.show(io::IO, loc::Location)

    s = @sprintf "%s(ID: %s, Lon: %.3f, Lat: %0.3f)" typeof(loc) string(loc.id) loc.lon loc.lat

    print(io, s)
end

###########
# Request #
###########

"""
Request represents a single area/point on Earth to be observed. It represents a location
for which it is possible to collect the entire area in a singel collect

Attributes:
- `lon::Float64` Longitude of request center [deg]
- `lat::Float64` Lattitude of iamge center [deg]
- `ecef::Array{Float64, 1}` Earth-Centered-Earth-Fixed Coordinates of request center [m]
- `look_angle_min::Float64` Minimum look angle of valid collect
- `look_angle_max::Float64` Maximum look angle of valid collect
- `collect_duration` Required time to well of request for a feasible collection
- `require_zero_doppler::Bool` Require that the collect occurs centered at zero Doppler offset
- `reward::Float64` Reward for for request collection
- `id::Union{UUID, Integer}` Unique request identifier
"""
mutable struct Request <: Location
    lon::Float64
    lat::Float64
    ecef::Array{Float64, 1}
    look_angle_min::Float64
    look_angle_max::Float64
    require_zero_doppler::Bool
    collect_duration::Float64
    reward::Float64
    id::Union{UUID, Integer}

    function Request(lon::Real, lat::Real ; id::Union{UUID, Integer}=uuid4(),
                    look_angle_min::Real=5.0, look_angle_max::Real=55.0,
                    require_zero_doppler::Bool=true, collect_duration::Real=1.0,
                    reward::Real=1.0, use_degrees::Bool=true)

        # Compute ECEF coordinates of scene center
        ecef = sGEODtoECEF([lon, lat], use_degrees=use_degrees::Real)

        # Generate new request object
        new(lon, lat, ecef, look_angle_min, look_angle_max,
            require_zero_doppler, collect_duration, reward, id)
    end
end

Base.hash(r::Request, h::UInt) = hash(r.id, hash(r.lon, hash(r.lat, hash(r.look_angle_min, hash(r.look_angle_max, hash(:Request, h))))))

function Base.show(io::IO, req::Request)

    s = @sprintf "Request(ID: %s, Lon: %.3f, Lat: %0.3f, Reward: %.2f, Duration: %.2f)" string(req.id) req.lon req.lat req.reward req.collect_duration

    print(io, s)
end

"""
Loads request data from a JSON file.

THe JSON file is expected to contain an array of JSON objects where each object
has the following format:

```json
{
    "lon": 22,
    "lat": 12.23,
    "reward": 42.0,
    "look_angle_min": 5.0,
    "look_angle_max": 55.0
}
```

Arguments:
- `file::String` Filepath to JSON file encoding request
"""
function load_requests(file::String; id_offset::Integer=0, require_zero_doppler::Bool=true)
    # Parse data file
    data = JSON.parsefile(file)

    requests = Array{Request, 1}(undef, length(data))

    for (i, req) in enumerate(data)
        lon = req["lon"]
        lat = req["lat"]
        look_angle_min = req["look_angle_min"]
        look_angle_max = req["look_angle_max"]
        id  = i + id_offset

        requests[i] = Request(lon, lat, id=id,
            look_angle_min=look_angle_min,
            look_angle_max=look_angle_max
        )
    end

    println("Loaded $(length(requests)) requests.")

    return requests
end

"""
Set the request collection duration for satellites in the constellation.
"""
function set_collect_duration(requests::Array{<:Request, 1}, collect_duration::Real)
    if collect_duration <= 0.0
        throw(ArgumentError("Invalid collection duration."))
    end

    for request in requests
        request.collect_duration = collect_duration
    end
end

#################
# GroundStation #
#################

"""
GroundStation represents a location which the satellite can use to communicate
with the ground.

Attributes:
- `lon::Float64` Longitude of request center [deg]
- `lat::Float64` Lattitude of iamge center [deg]
- `ecef::Array{Float64, 1}` Earth-Centered-Earth-Fixed Coordinates of request center [m]
- `elevation_min::Float64` Minimum look angle of valid collect
- `id::Union{UUID, Integer}` Unique station identifier
"""
mutable struct GroundStation <: Location
    lon::Float64
    lat::Float64
    ecef::Array{Float64, 1}
    elevation_min::Float64
    duration_min::Float64
    id::Union{UUID, Integer}

    function GroundStation(lon::Real, lat::Real ; id::Union{UUID, Integer}=uuid4(),
                    elevation_min::Real=5.0, duration_min::Real=0.0,
                    use_degrees::Bool=true)

        # Compute ECEF coordinates of scene center
        ecef = sGEODtoECEF([lon, lat], use_degrees=use_degrees)

        # Generate new request object
        new(lon, lat, ecef, elevation_min, duration_min, id)
    end
end

Base.hash(s::GroundStation, h::UInt) = hash(s.id, hash(s.lon, hash(s.lat, hash(s.elevation_min, hash(:GroundStation, h)))))

# function Base.isequal(ll::GroundStation, rl::GroundStation)
#     return (
#         (ll.id == rl.id) &&
#         (ll.lat == rl.lat) &&
#         (ll.lon == rl.lon) &&
#         (ll.ecef == rl.ecef) && true
#     )
# end
# Base.:(==)(ll::GroundStation, rl::GroundStation) = Base.isequal(ll, rl)

function Base.show(io::IO, sta::GroundStation)

    s = @sprintf "GroundStation(ID: %s, Lon: %.3f, Lat: %0.3f)" string(sta.id) sta.lon sta.lat

    print(io, s)
end

"""
Loads request data from a JSON file.

THe JSON file is expected to contain an array of JSON objects where each object
has the following format:

```json
{
    "lon": 22,
    "lat": 12.23,
    "elevation_min": 5.0
}
```

Arguments:
- `file::String` Filepath to JSON file encoding request
"""
function load_stations(file::String; id_offset::Integer=0)
    data = JSON.parsefile(file)

    # Initialize array of Requests
    stations = Array{GroundStation, 1}(undef, length(data))

    for (i, sta) in enumerate(data)
        lon = sta["lon"]
        lat = sta["lat"]
        elevation_min = sta["elevation_min"]
        id  = i + id_offset

        stations[i] = GroundStation(lon, lat, id=id, elevation_min=elevation_min)
    end

    println("Loaded $(length(stations)) stations.")

    return stations
end

"""
Set minimum access duration for all stations
"""
function set_station_duration_min(stations::Array{<:GroundStation, 1}, duration_min::Real)
    if duration_min <= 0.0
        throw(ArgumentError("Invalid minimum contact duration duration."))
    end

    for station in stations
        station.duration_min = duration_min
    end
end

###############
# Opportunity #
###############

"""
Opportunity represents an "opportunity" for action for a spacecraft

Common Attributes:
- `id::Integer` Unique Identifier for action
- `orbit::TLE` Orbit object associated with this collect
- `location::Location` Location object on the ground for interaction
- `t_start::Union{Epoch, Real}` Start of acquisition window
- `t_mid::Union{Epoch, Real}` Mid-time of acquisition window
- `t_end::Union{Epoch, Real}` End of possible acquisition window
- `duration::Real` Length of action
"""
abstract type Opportunity end

"""
"""
@with_kw mutable struct Done <: Opportunity
    t_start::Union{Epoch, Real} = Inf
end

"""
No-operation action. Advances state without changing other values
"""
@with_kw mutable struct Noop <: Opportunity
    id::Integer = 0
    t_start::Union{Epoch, Real} = 0.0
end

function Base.show(io::IO, noop::Noop)

    s = @sprintf "Noop(time: %s)" noop.t_start

    print(io, s)
end

@with_kw mutable struct Sunpoint <: Opportunity
    id::Integer = 0
    t_start::Union{Epoch, Real} = 0.0
end

function Base.show(io::IO, sp::Sunpoint)

    s = @sprintf "Sunpoint(time: %s)" sp.t_start

    print(io, s)
end

"""
Opportunity type for a Request collection.
"""
@with_kw mutable struct Collect <: Opportunity
    id::Integer
    spacecraft::Spacecraft
    location::Union{Request, Nothing}
    t_start::Union{Epoch, Real}
    t_mid::Union{Epoch, Real}
    t_end::Union{Epoch, Real}
    duration::Float64
    reward::Float64
    nr::Int

    function Collect(t_start::Union{Epoch, Real}, t_end::Union{Epoch, Real};
        id::Integer=0, spacecraft=nothing::Union{Spacecraft, Nothing},
        location=nothing::Union{Location, Nothing},
        nr::Int=0)

        duration = t_end - t_start
        t_mid    = t_start + duration/2.0
        new(id, spacecraft, location, t_start, t_mid, t_end, duration, location.reward, nr)
    end
end

Base.hash(c::Collect, h::UInt) = hash(c.id, hash(c.spacecraft.id, hash(c.duration, hash(:Collect, h))))

function JSON.lower(col::Collect)
    return Dict(
        "id" => col.id,
        "location" => col.location.id,
        "lon" => col.location.lon,
        "lat" => col.location.lat,
        "t_start" => string(col.t_start),
        "t_end"  => string(col.t_end),
        "duration" => col.duration
    )
end

function Base.isequal(cl::Collect, cr::Collect)
    return (
        (cl.id == cr.id) &&
        ((cl.spacecraft == nothing && cr.spacecraft == nothing) ||
         ((cl.spacecraft != nothing && cr.spacecraft != nothing &&
          (cl.spacecraft.id == cr.spacecraft.id)))) &&
        ((cl.location == nothing && cr.location == nothing) ||
         ((cl.location != nothing && cr.location != nothing &&
          (cl.location.id == cr.location.id)))) &&
        (cl.t_start == cr.t_start) &&
        (cl.t_end == cr.t_end) && true
    )
end
Base.:(==)(cl::Collect, cr::Collect) = Base.isequal(cl, cr)

function Base.show(io::IO, col::Collect)

    if typeof(col.t_start) == Epoch
        s = @sprintf "Collect(ID: %s, Spacecraft %s, Request: %s, t_start: %s, t_end: %s, Reward: %.2f, Remaining: %d)" string(col.id) string(col.spacecraft.id) string(col.location.id) col.t_start col.t_end col.reward col.nr
    else
        s = @sprintf "Collect(ID: %s, Spacecraft %s, Request: %s, t_start: %.3f, t_end: %.3f, Reward: %.2f, Remaining: %d)" string(col.id) string(col.spacecraft.id) string(col.location.id) col.t_start col.t_end col.reward col.nr
    end

    print(io, s)
end

"""
Opportunity type for a ground contact.
"""
mutable struct Contact <: Opportunity
    id::Integer
    spacecraft::Spacecraft
    location::Union{GroundStation, Nothing}
    t_start::Union{Epoch, Real}
    t_mid::Union{Epoch, Real}
    t_end::Union{Epoch, Real}
    duration::Float64
    nr::Int

    function Contact(t_start::Union{Epoch, Real}, t_end::Union{Epoch, Real};
        id::Integer=0, spacecraft=nothing::Union{Spacecraft, Nothing},
        location=nothing::Union{Location, Nothing}, nr::Int=0)

        duration = t_end - t_start
        t_mid    = t_start + duration/2.0
        new(id, spacecraft, location, t_start, t_mid, t_end, duration, nr)
    end
end

Base.hash(c::Contact, h::UInt) = hash(c.id, hash(c.spacecraft.id, hash(c.duration, hash(:Contact, h))))

function JSON.lower(con::Contact)
    return Dict(
        "id" => con.id,
        "location" => con.location.id,
        "lon" => con.location.lon,
        "lat" => con.location.lat,
        "t_start" => string(con.t_start),
        "t_end"  => string(con.t_end),
        "duration" => con.duration
    )
end

function Base.isequal(cl::Contact, cr::Contact)
    return (
        (cl.id == cr.id) &&
        ((cl.spacecraft == nothing && cr.spacecraft == nothing) ||
         ((cl.spacecraft != nothing && cr.spacecraft != nothing &&
          (cl.spacecraft.id == cr.spacecraft.id)))) &&
        ((cl.location == nothing && cr.location == nothing) ||
         ((cl.location != nothing && cr.location != nothing &&
          (cl.location.id == cr.location.id)))) &&
        (cl.t_start == cr.t_start) &&
        (cl.t_end == cr.t_end) && true
    )
end
Base.:(==)(cl::Contact, cr::Contact) = Base.isequal(cl, cr)

function Base.show(io::IO, con::Contact)

    minutes = floor(con.duration/60)
    seconds = con.duration - 60*minutes

    if typeof(con.t_start) == Epoch
        s = @sprintf "Contact(ID: %s,  Spacecraft %s, Station: %s, t_start: %s, t_end: %s, Duration: %02dm %02ds, Remaining: %d)" string(con.id) string(con.spacecraft.id) string(con.location.id) con.t_start con.t_end minutes seconds con.nr
    else
        s = @sprintf "Contact(ID: %s,  Spacecraft %s, Station: %s, t_start: %.3f, t_end: %.3f, Duration: %02dm %02ds, Remaining: %d)" string(con.id) string(con.spacecraft.id) string(con.location.id) con.t_start con.t_end minutes seconds con.nr
    end

    print(io, s)
end

#############
# MDP State #
#############

@with_kw struct SatMDPState
    time::Union{Epoch, Float64}
    last_cdo_action::Opportunity
    last_action::Opportunity
    request_ids::Array{Int, 1} = Int[]
    power::Float64 = 1.0
    data::Float64 = 0.0
end
Base.hash(s::SatMDPState, h::UInt) = hash(s.time, hash(s.last_action, hash(s.power, hash(s.data, hash(s.request_ids, hash(:SatMDPState, h))))))

function Base.show(io::IO, state::SatMDPState)

    s = @sprintf "SatMDPState(Time: %s, Last CDO Action: %s, Action: %s, Requests: %d, Power: %.2f, Data: %.2f)" string(state.time) string(state.last_cdo_action) string(state.last_action) length(state.request_ids) state.power state.data

    print(io, s)
end

function POMDPs.isequal(sl::SatMDPState, sr::SatMDPState)
    return (
        (sl.time == sr.time) &&
        (sl.last_action == sr.last_action) &&
        (sl.last_cdo_action == sr.last_cdo_action) &&
        (sl.power == sr.power) &&
        (sl.data == sr.data) &&
        (sl.request_ids == sr.request_ids)
    )
end

####################
# Planning Problem #
####################

@with_kw mutable struct SatPlanningProblem <: POMDPs.MDP{SatMDPState, Opportunity}
    # General Planning Settings
    t_start::Epoch
    t_end::Epoch

    # Spacecraft Settings
    spacecraft::Array{Spacecraft, 1} = Spacecraft[]
    constraints::Array{Function,1} = Function[]

    # Planning Locations
    locations::Array{<:Location, 1} = Location[]
    stations::Array{GroundStation, 1} = GroundStation[]
    requests::Array{Request, 1} = Request[]

    # Planning Opportunities
    opportunities::Array{<:Opportunity, 1} = Opportunity[]
    contacts::Array{Contact, 1} = Contact[]
    collects::Array{Collect, 1} = Collect[]

    # Lookup Tables
    lt_locations::Dict{Integer, Location} = Dict{Integer, Location}()
    lt_opportunities::Dict{Integer, Opportunity} = Dict{Integer, Opportunity}()
    lt_contacts::Dict{Integer, Contact} = Dict{Integer, Contact}()
    lt_collects::Dict{Integer, Collect} = Dict{Integer, Collect}()
    lt_loc_opps::Dict{Integer, AbstractVector{Integer}} = Dict{Integer, AbstractVector{Integer}}()

    # Action Lookup Table
    actions::AbstractVector{Opportunity} = Opportunity[]
    lt_actions::Dict{Integer, Opportunity} = Dict{Integer, Opportunity}()
    lt_feasible_actions::Dict{Tuple{Integer, Integer}, AbstractVector{Opportunity}} = Dict{Tuple{Opportunity, Opportunity}, AbstractVector{Opportunity}}()

    # Solver Parameters - General
    solve_allow_repeats::Bool = false
    solve_discount::Real   = 0.99
    solve_depth::Int    = 3
    solve_breadth::Int  = 0
    solve_horizon::Real = 90.0

    # # Solver Parameters - MDP
    # mdp_reward_scarcity::Bool = false

    # Solver Parameters - MCTS
    mcts_n_iterations::Int = 50
    mcts_exploration_constant::Real = 1.0
end

"""
Clear problem locations.
"""
function clear_locations(problem::SatPlanningProblem)
    # Reset Locations
    problem.locations = Location[]
    problem.stations = GroundStation[]
    problem.requests = Request[]
end

"""
Clear problem opportunities.
"""
function clear_opportunities(problem::SatPlanningProblem)
    # Reset Opportunities
    problem.opportunities = Opportunity[]
    problem.contacts = Contact[]
    problem.collects = Collect[]

    problem.lt_opportunities = Dict{Integer, Opportunity}()
    problem.lt_contacts      = Dict{Integer, Contact}()
    problem.lt_collects      = Dict{Integer, Collect}()
    problem.lt_loc_opps      = Dict{Integer, Array{Integer, 1}}()

    # Action Lookup Table
    lt_actions = Opportunity[]
    lt_feasible_actions = Dict{Tuple{Opportunity, Opportunity}, AbstractVector{Opportunity}}()
end

"""
Clear all problem opportunities
"""
function clear_all(problem::SatPlanningProblem)
    clear_locations(problem)
    clear_opportunities(problem)
end

"""
Add locations to planning problem.
"""
function add_locations(problem::SatPlanningProblem, locations::Array{<:Location})
    if length(locations) == 0
        @warn "No locations in array. Nothing to add..."
        return
    end

    # TODO: Check for duplicate IDs in requests (assumes already unique)

    # Add location to problem
    push!(problem.locations, locations...)

    # Create separate arrays of contacts and requests
    push!(problem.stations, filter(x -> typeof(x) == GroundStation, locations)...)
    push!(problem.requests, filter(x -> typeof(x) == Request, locations)...)

    # Update location lookup table
    for loc in problem.locations
        problem.lt_locations[loc.id] = loc
    end

    return
end


text:
using Base.Test, MathProgComplex

## Application of the Moment–SOS Approach to Global Optimization of the OPF Problem
# C. Josz, J. Maeght, P. Panciatici, and J. Ch. Gilbert
# https://arxiv.org/pdf/1311.6370.pdf
# Data from table 3, p. 13.

@testset "WB5 real formulation - dense - order 2" begin
    sols = SortedDict(-30.80 => (2, 945.83, 945.83),
                      -20.51 => (2, 1146.48, 954.82),
                      -10.22 => (2, 1209.11, 963.83),
                       00.07 => (2, 1267.79, 972.85),
                       10.36 => (2, 1323.86, 981.89),
                       20.65 => (2, 1377.97, 990.95),
                       30.94 => (2, 1430.54, 1005.13),
                       41.23 => (2, 1481.81, 1033.07),
                       51.52 => (2, 1531.97, 1070.39),
                       61.81 => (1, 1114.90, 1114.90))

    testfolder = joinpath(Pkg.dir("MathProgComplex"), "Mosek_runs", "tests")
    @testset "q5min=$q5min, rmeqs=$rmeqs" for (q5min, (dmax, obj_fixedphase, obj_rankrel)) in sols, rmeqs in Set([true])

        # info("Working on WB5, no eqs, q5min=$q5min, d=$d")
        problem = buildPOP_WB5(q5min=q5min, rmeqs=rmeqs)

        logpath = joinpath(testfolder, "WB5_q5min_$(q5min)_ordercv")
        ispath(logpath) && rm(logpath, recursive=true); mkpath(logpath)
        println("Saving file at $logpath")

        relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                            d = dmax)

        primobj, dualobj = run_hierarchy(problem, relax_ctx, logpath, save_pbs=true);
        @show primobj, dualobj, obj_fixedphase
        @test primobj ≈ obj_fixedphase atol=1e-2
        @test dualobj ≈ primobj atol=mosek_optgap*min(abs(primobj), abs(dualobj))
    end
end


text:
clearconsole()
# using Base # not required
using FileIO
using MatlabCompat
using MatlabCompat.Io
using MatlabCompat.ImageTools
import MatlabCompat.MathTools
import MatlabCompat.MathTools: max
using Test
using Images
using FixedPointNumbers
# using Tk

const TEST_FOLDER = dirname(@__FILE__)
const JANUS_PATH = joinpath(TEST_FOLDER, "janus.m");
const TESTMAT_PATH = joinpath(TEST_FOLDER, "test.mat");

# basic test
println("foo test")
x = function ()
    return 1 == 1
end
@test x()

################################
#     ImageTools Module Tests
################################

# make a test image
array = [0.1  0.001  0  0;
0.1  0.401  0  0;
0.1  0.01   0  0;
0.1  0.901  0  0.6]; # should be Float : Float32, Float64, etc

# Any array can be treated as an Image. However, in graphical environments, only arrays with Colorant element types are automatically displayed as images.
#  Coloran typest=(Gray, RGB, ARGB, etc.)

# Syntax 1 - Array{Gray}
img1 =Gray.(array);
#  # calculates a grayscale representation of a color image using the Rec 601 luma. (https://juliaimages.org/latest/function_reference/#Color-conversion-1)
# imgMet=ImageMeta(img, time="after noon");
println(typeof(img1)) # gives Array{Gray{Float64},2}

# Syntax 2 - ReinterpretArray
img2 = colorview(Gray, array);
# returns a view of the numeric array A, interpreting successive elements of A as if they were channels of Colorant C. (https://juliaimages.org/latest/function_reference/#ImageCore.colorview)
imgMet=ImageMeta(img2, time="after noon"); # to add ImageMeta
println(typeof(img2)) # gives Base.ReinterpretArray{Gray{Float64},2,Float64,Array{Float64,2}}

# Syntax 3 - ReinterpretArray + ImageMeta
imgMet = ImageMeta(colorview(Gray, array), time="after noon")
 # returns an ImageMeta type, which holds both data() with Colorant element type, and properties() (Dict type)
properties(imgMet)["location"] = "university"; # to access and extract properties
img3=data(imgMet); # to access and extract data
println(typeof(img3)) # gives Base.ReinterpretArray{Gray{Float64},2,Float64,Array{Float64,2}}

# Syntax 4 - Array + ImageMeta
# this does not work with this module
imgMet = ImageMeta(array, time="after noon")
 # returns an ImageMeta type, which holds both data() but array has Float type instead of Colorant, and properties() (Dict type)
properties(imgMet)["location"] = "university"; # to access and extract properties
img4=data(imgMet); # to access and extract data
println(typeof(img4)) # gives Array{Float64,2}
# this synatx does

img=img1; # choose img for test

# test the fuctions
println("setting true tests results");
trueThreshold = 0.11100752480522075;
trueBWSum= 3;

println("testing graythresh()");
# @test graythresh(img) == trueThreshold;
graythresh(img) == trueThreshold;

println("testing im2bw()");
# @test sum( reinterpret( Float32, float32( im2bw(img, trueThreshold) ) )[:] ) == trueBWSum;
sum( reinterpret( Float32, float32( im2bw(img, trueThreshold) ) )[:] ) == trueBWSum;

println("testing imread()");
img5 = MatlabCompat.imread("http://matlabcompat.github.io/img/example.tif");
expected = nothing
expected = "Images.Image{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2,Array{ColorTypes.Gray{FixedPointNumbers.UFixed{UInt8,8}},2}}"

@test "$(typeof(img5))" == expected;

################################
#     Support Module Tests
################################

println("testing im2mat()");
mat = im2mat(img)
@test "$(typeof(mat))" == "Array{Float32,2}";

println("testing mat2im()");
expected = nothing
if (false)
    expected = "Image{Float64,2,Array{Float64,2}}"
else
    expected = "Images.Image{Float64,2,Array{Float64,2}}"
end
@test "$(typeof(mat2im(img)))" == expected;

println("testing rosetta()");
@test "$(typeof(rosetta(JANUS_PATH)))" == "Array{UTF8String,1}";

################################
#     Io Module Tests
################################

println("testing load()");
println("$(typeof(MatlabCompat.load(TESTMAT_PATH)))")
@test "$(typeof(MatlabCompat.load(TESTMAT_PATH)))" == "Dict{ASCIIString,Any}";
# "$(typeof(A))" == "Dict{Any,Any}"

################################
#     StringTools Module Tests
################################
println("testing disp()");
@test disp("test") == true;

println("testing num2str()");
@test num2str(3) == "3";

println("testing strcat()");
@test strcat("1", "2", "3") == "123";

################################
#     MathTools Module Tests
################################
matrix = [1 2 3];

println("testing numel()");
@test numel(matrix[:]) == 3;

println("testing max()");

@test max(matrix) == maximum(matrix);


text:
struct Rxn{T<:Real}
    id::String
    S::Vector{T}
    mets::Vector
    c::T
    lb::T
    ub::T

    function Rxn{T}(id::String;
            mets = [],
            S = T[],
            c = zero(T),
            lb = zero(T),
            ub = zero(T)
        ) where {T<:Real}

        length(S) != length(mets) && error("'S' and 'mets' must have the same length")
        new{T}(id, S, mets, c, lb, ub)
    end
end

Rxn(id::String; kwargs...) = Rxn{Float64}(id; kwargs...)


text:
@testset "ElementaryPath" begin
    @testset "Transforming a solution" begin
        g = path_digraph(3)
        add_edge!(g, 1, 3)
        rewards = Dict(Edge(1, 2) => 1.0, Edge(2, 3) => 1.0, Edge(1, 3) => 1.0)
        weights = Dict(Edge(1, 2) => 1, Edge(2, 3) => 1, Edge(1, 3) => 4)

        i = MinimumBudget(ElementaryPathInstance(g, rewards, 1, 3), weights, 4)
        path = [Edge(1, 3)]
        solutions = Dict{Int, Vector{Edge{Int}}}(
            0 => [Edge(1, 2), Edge(2, 3)],
            1 => [Edge(1, 2), Edge(2, 3)],
            2 => [Edge(1, 2), Edge(2, 3)],
            3 => [Edge(1, 3)],
            4 => [Edge(1, 3)],
        )
        d = BudgetedElementaryPathSolution(i, path, solutions)

        warn_msg = "The requested maximum budget 5 is higher than the instance's minimum budget 4. Therefore, some values have not been computed and are unavailable."
        @test_logs (:warn, warn_msg) Kombinator.paths_all_budgets_as_tuples(
            d,
            5,
        )
        sol = Dict(
            0 => [(1, 2), (2, 3)],
            1 => [(1, 2), (2, 3)],
            2 => [(1, 2), (2, 3)],
            3 => [(1, 3)],
            4 => [(1, 3)],
        )
        @test Kombinator.paths_all_budgets_as_tuples(d, 4) == sol

        @test_logs (:warn, warn_msg) Kombinator.paths_all_budgets(d, 5)
        sol = Dict(
            0 => [Edge(1, 2), Edge(2, 3)],
            1 => [Edge(1, 2), Edge(2, 3)],
            2 => [Edge(1, 2), Edge(2, 3)],
            3 => [Edge(1, 3)],
            4 => [Edge(1, 3)],
        )
        @test Kombinator.paths_all_budgets(d, 4) == sol
    end

    @testset "Elementary paths" begin
        @testset "Interface" begin
            g = path_digraph(3)
            rewards = Dict(Edge(1, 2) => 1.0, Edge(2, 3) => 1.0)

            @testset "Global interface" begin
                i = ElementaryPathInstance(g, rewards, 1, 2)
                @test objective(i) == Maximise()
            end

            @testset "Source not in graph" begin
                @test_throws ErrorException ElementaryPathInstance(g, rewards, 50, 2)
            end

            @testset "Destination not in graph" begin
                @test_throws ErrorException ElementaryPathInstance(g, rewards, 1, 50)
            end

            @testset "Source is destination" begin
                @test_throws ErrorException ElementaryPathInstance(g, rewards, 1, 1)
            end

            @testset "Source is destination" begin
                @test_throws ErrorException ElementaryPathInstance(g, rewards, 1, 1)
            end

            @testset "An edge has a reward but is not in the graph" begin
                rewards = Dict(Edge(1, 2) => 1.0, Edge(2, 3) => 1.0, Edge(1, 3) => 1.0)
                @test_throws ErrorException ElementaryPathInstance(g, rewards, 1, 2)
            end

            @testset "Positive-reward cycle" begin
                g = path_digraph(3)
                add_edge!(g, 2, 1)
                rewards = Dict(
                    Edge(1, 2) => 1.0,
                    Edge(2, 3) => 3.0,
                    Edge(2, 1) => 25.0,
                )

                i = ElementaryPathInstance(g, rewards, 1, 3)
                @test_logs (
                    :warn,
                    "The graph contains a positive-cost cycle around edge 2 -> 1.",
                ) solve(i, BellmanFordAlgorithm())
            end

            @testset "Copy" begin
                g = path_digraph(3)
                add_edge!(g, 2, 1)
                rewards = Dict(Edge(1, 2) => 1.0, Edge(2, 3) => 3.0, Edge(2, 1) => 25.0)

                i = ElementaryPathInstance(g, rewards, 1, 3)
                i2 = copy(i)
                @test i.graph == i2.graph
                @test i.rewards == i2.rewards
                @test i.src == i2.src
                @test i.dst == i2.dst
                @test i.objective == i2.objective
            end

            @testset "Approximation" begin
                i = ElementaryPathInstance(g, rewards, 1, 3)

                @test approximation_ratio(i, DynamicProgramming()) == 1.0
                @test approximation_term(i, DynamicProgramming()) == 0.0

                @test approximation_ratio(i, DefaultLinearFormulation(Cbc.Optimizer)) == 1.0
                @test approximation_term(i, DefaultLinearFormulation(Cbc.Optimizer)) == 0.0
            end

            @testset "Make solution" begin
                i = ElementaryPathInstance(g, rewards, 1, 3)
                d = Dict(Edge(1, 3) => 0.4, Edge(1, 2) => 0.8, Edge(2, 3) => 0.7)
                s = make_solution(i, d)

                @test s.instance == i
                @test length(s.variables) == 2
                @test Edge(1, 3) ∉ s.variables
                @test Edge(1, 2) ∈ s.variables
                @test Edge(2, 3) ∈ s.variables
            end

            @testset "Solution helpers" begin
                i = ElementaryPathInstance(g, rewards, 1, 3)

                # Infeasible solution.
                s = ElementaryPathSolution(i, Edge{Int}[])
                @test value(s) == -Inf

                # Valid solution.
                d = Dict(Edge(1, 3) => 0.4, Edge(1, 2) => 0.8, Edge(2, 3) => 0.7)
                s = make_solution(i, d)

                @test value(s) == 4
            end
        end

        @testset "Basic" begin
            @testset "Directed path graph" begin
                g = path_digraph(3)
                rewards = Dict(Edge(1, 2) => 1.0, Edge(2, 3) => 1.0)
                i = ElementaryPathInstance(g, rewards, 1, 3)

                @testset "Bellman-Ford" begin
                    s = solve(i, BellmanFordAlgorithm())
                    @test s.instance == i
                    @test s.variables == [Edge(1, 2), Edge(2, 3)]
                end

                if !is_travis
                    @testset "Linear programming" begin
                        s = solve(
                            i,
                            DefaultLinearFormulation(Gurobi.Optimizer), # Cbc unsupported.
                        )
                        @test s.instance == i
                        @test s.variables == [Edge(1, 2), Edge(2, 3)]
                    end
                end
            end

            @testset "High-weight edge" begin
                g = path_digraph(3)
                add_edge!(g, 1, 3)
                rewards = Dict(
                    Edge(1, 2) => 1.0,
                    Edge(2, 3) => 1.0,
                    Edge(1, 3) => 4.0,
                )
                weights =
                    Dict(Edge(1, 2) => 1, Edge(2, 3) => 1, Edge(1, 3) => 4)

                i = ElementaryPathInstance(g, rewards, 1, 3)

                @testset "Bellman-Ford" begin
                    s = solve(i, BellmanFordAlgorithm())
                    @test s.instance == i
                    @test s.variables == [Edge(1, 3)]
                end

                if !is_travis
                    @testset "Linear programming" begin
                        s = solve(
                            i,
                            DefaultLinearFormulation(Gurobi.Optimizer),
                        ) # Cbc unsupported.
                        @test s.instance == i
                        @test s.variables == [Edge(1, 3)]
                    end
                end
            end
        end
    end

    @testset "Budgeted elementary paths" begin
        @testset "Interface" begin
            graph = path_digraph(3)
            rewards = Dict(Edge(1, 2) => 1.0, Edge(2, 3) => 1.0)
            weights = Dict(Edge(1, 2) => 1, Edge(2, 3) => 1)

            @testset "Global interface" begin
                i = MinimumBudget(ElementaryPathInstance(graph, rewards, 1, 2), weights)
                @test objective(i) == Maximise()
            end

            @testset "More edges have rewards than weights" begin
                g = path_digraph(4)
                r = Dict(Edge(1, 2) => 1.0, Edge(2, 3) => 1.0, Edge(4, 1) => 1.0)
                w = Dict(Edge(1, 2) => 1, Edge(2, 3) => 1)
                @test_throws ErrorException MinimumBudget(ElementaryPathInstance(g, r, 1, 2), w)
            end

            @testset "Negative weight or budget" begin
                w = Dict(Edge(1, 2) => 1, Edge(2, 3) => 1)
                @test_throws ErrorException MinimumBudget(ElementaryPathInstance(graph, rewards, 1, 2), w, -1)

                w[Edge(4, 2)] = 1
                @test_throws ErrorException MinimumBudget(ElementaryPathInstance(graph, rewards, 1, 2), w)
                w[Edge(4, 2)] = -1
                @test_throws ErrorException MinimumBudget(ElementaryPathInstance(graph, rewards, 1, 2), w)
            end

            @testset "Copy" begin
                i = MinimumBudget(ElementaryPathInstance(graph, rewards, 1, 2), weights)
                i2 = copy(i)
                @test i.weights == i2.weights
                @test i.min_budget == i2.min_budget
                @test i.compute_all_values == i2.compute_all_values
                @test i.instance.graph == i2.instance.graph
                @test i.instance.rewards == i2.instance.rewards
                @test i.instance.src == i2.instance.src
                @test i.instance.dst == i2.instance.dst
                @test i.instance.objective == i2.instance.objective
            end

            @testset "Approximation" begin
                i = MinimumBudget(ElementaryPathInstance(graph, rewards, 1, 2), weights)

                @test approximation_ratio(i, DynamicProgramming()) == 1.0
                @test approximation_term(i, DynamicProgramming()) == 0.0

                @test approximation_ratio(i, DefaultLinearFormulation(Cbc.Optimizer)) == 1.0
                @test approximation_term(i, DefaultLinearFormulation(Cbc.Optimizer)) == 0.0
            end

            @testset "Make solution" begin
                i = MinimumBudget(ElementaryPathInstance(graph, rewards, 1, 2), weights)
                d = Dict(Edge(1, 3) => 0.4, Edge(1, 2) => 0.8, Edge(2, 3) => 0.7)
                s = make_solution(i, d)

                @test s.instance == i
                @test length(s.variables) == 2
                @test Edge(1, 3) ∉ s.variables
                @test Edge(1, 2) ∈ s.variables
                @test Edge(2, 3) ∈ s.variables
            end

            @testset "Solution helpers" begin
                i = MinimumBudget(ElementaryPathInstance(graph, rewards, 1, 2), weights)

                # Infeasible solution.
                s = BudgetedElementaryPathSolution(i, Edge{Int}[])
                @test value(s) == -Inf

                # Valid solution.
                d = Dict(Edge(1, 3) => 0.4, Edge(1, 2) => 0.8, Edge(2, 3) => 0.7)
                s = make_solution(i, d)

                @test value(s) == 2
            end
        end

        @testset "Basic" begin
            @testset "Directed path graph" begin
                g = path_digraph(3)
                rewards = Dict(Edge(1, 2) => 1.0, Edge(2, 3) => 1.0)
                weights = Dict(Edge(1, 2) => 1, Edge(2, 3) => 1)

                i = MinimumBudget(
                    ElementaryPathInstance(g, rewards, 1, 3),
                    weights,
                    4,
                )

                @testset "Bellman-Ford" begin
                    d = solve(i, BellmanFordAlgorithm())

                    @test d.instance == i
                    @test d.variables == [] # No path with a total weight of at least 4.

                    for β in [0, 1, 2]
                        @test d.solutions[β] == [Edge(1, 2), Edge(2, 3)]
                    end
                    for β in [3, 4]
                        @test d.solutions[β] == []
                    end
                end

                if !is_travis
                    @testset "Linear programming" begin
                        s = solve(
                            i,
                            DefaultLinearFormulation(Gurobi.Optimizer),
                        ) # Cbc unsupported.
                        @test s.instance == i
                        @test s.variables == []
                    end
                end
            end

            @testset "High-weight edge" begin
                g = path_digraph(3)
                add_edge!(g, 1, 3)
                rewards = Dict(
                    Edge(1, 2) => 1.0,
                    Edge(2, 3) => 1.0,
                    Edge(1, 3) => 1.0,
                )
                weights =
                    Dict(Edge(1, 2) => 1, Edge(2, 3) => 1, Edge(1, 3) => 4)

                i = MinimumBudget(
                    ElementaryPathInstance(g, rewards, 1, 3),
                    weights,
                    4,
                )

                @testset "Bellman-Ford" begin
                    s = solve(i, BellmanFordAlgorithm())
                    @test s.instance == i
                    @test s.variables == [Edge(1, 3)]
                end

                if !is_travis
                    @testset "Linear programming" begin
                        s = solve(
                            i,
                            DefaultLinearFormulation(Gurobi.Optimizer),
                        ) # Cbc unsupported.
                        @test s.instance == i
                        @test s.variables == [Edge(1, 3)]
                    end
                end
            end
        end

        @testset "Conformity" begin
            # More advanced tests to ensure the algorithm works as expected.
            g = complete_digraph(3)
            rewards = Dict(
                Edge(1, 2) => 1.0,
                Edge(3, 1) => -1.0,
                Edge(3, 2) => -1.0,
                Edge(2, 3) => 1.0,
                Edge(2, 1) => -1.0,
                Edge(1, 3) => 0.0,
            )
            weights = Dict(
                Edge(1, 2) => 0,
                Edge(3, 1) => 0,
                Edge(3, 2) => 0,
                Edge(2, 3) => 0,
                Edge(2, 1) => 0,
                Edge(1, 3) => 2,
            )
            i = MinimumBudget(
                ElementaryPathInstance(g, rewards, 1, 3),
                weights,
                2,
                compute_all_values=true,
            )

            @testset "Bellman-Ford" begin
                warn_msg = "The graph contains a positive-cost cycle around edge 3 -> 1."
                d = @test_logs (:warn, warn_msg) solve(
                    i,
                    BellmanFordAlgorithm(),
                )

                @test d.variables == [Edge(1, 3)]

                @test d.solutions[0] == [Edge(1, 2), Edge(2, 3)]
                @test d.solutions[1] == [Edge(1, 3)]
                @test d.solutions[2] == [Edge(1, 3)]
            end

            if !is_travis
                @testset "Linear programming" begin
                    l = solve(
                        i,
                        DefaultLinearFormulation(Gurobi.Optimizer),
                    ) # Cbc unsupported.
                    # All solutions have the same destination! Only paths from 1 to 3, unlike DP.
                    @test l.solutions[0] == [Edge(1, 2), Edge(2, 3)]
                    @test l.solutions[1] == [Edge(1, 3)]
                    @test l.solutions[2] == [Edge(1, 3)]
                end
            end
        end
    end
end


text:
"""
jupyter_main(vcf_filename,saving_options,variant_filters,sample_selection,plotting_options)

filters, plots visualization, and saves as figure.
utilizes all global variables set in first cell of jupyter notebook
"""
function jupyter_main(vcf_filename,saving_options,variant_filters,sample_selection,plotting_options)

    genomic_range=variant_filters[1]
    pass_filter=variant_filters[2]
    positions_list=variant_filters[3]
    group_samples=sample_selection[1]
    select_samples=sample_selection[2]
    heatmap_options=plotting_options[1]
    read_depth_scatter_plots=plotting_options[2]
    y_axis_label_option=plotting_options[3]
    x_axis_labels=plotting_options[4]
    heatmap_title=saving_options[1]
    save_format=saving_options[2]
    output_directory=saving_options[3]
    num_array=saving_options[4]
    remote_option=saving_options[5]

    #create vcf reader object
    println("Reading $vcf_filename")
    println()
    reader = VCF.Reader(open(vcf_filename, "r"))

    #make sample names object
    sample_names = get_sample_names(reader)

    #store save format and create output directory if it doesn't exist yet

    #save_ext = save_format
    VariantVisualization.checkfor_outputdirectory(output_directory)
    output_directory=output_directory

    save_ext = save_format
    #output_directory=output_directory

    #store plot label options
    if occursin("true",x_axis_labels)
        x_axis_label_option = true
    else
        x_axis_label_option = false
    end

    #check for filters and apply then show stats again

    #no filter
    if !occursin(".csv",positions_list) && !occursin("chr",genomic_range) && !occursin("true",pass_filter)
        println("no filters applied. Large vcf files will take a long time to process and heatmap visualizations will not be useful at this scale.")
        println()
        println("Loading VCF file into memory for visualization")

        sub = Array{Any}(undef,0)
        for record in reader
            push!(sub,record)
        end

        number_rows=size(sub,1)
    end

    #all filter combos

    #pass_filter and genomic_range and list
    if occursin("true",pass_filter) && occursin("chr",genomic_range) && occursin(".csv",positions_list)
        sig_list =  load_siglist(positions_list)
        sub = VariantVisualization.pass_genomic_range_siglist_filter(vcf_filename, sig_list, genomic_range)
        number_rows = size(sub,1)
        println("Selected $number_rows variants with Filter status: PASS, that match list of chromosome positions of interest, and are within chromosome range: $genomic_range")
    end

    #pass_filter and genomic_range
    if occursin("true",pass_filter) && occursin("chr",genomic_range) && !occursin(".csv",positions_list) && !occursin("true",pass_filter)
        sub = VariantVisualization.pass_genomic_range_filter(reader, genomic_range,vcf_filename)
        number_rows = size(sub,1)
        println("Selected $number_rows variants with Filter status: PASS and are within chromosome range: $genomic_range")
    end

    #pass_filter and list
    if occursin("true",pass_filter) && occursin(".csv",positions_list) && !occursin("chr",genomic_range)
        sig_list =  load_siglist(positions_list)
        sub = VariantVisualization.pass_siglist_filter(vcf_filename, sig_list)
        number_rows = size(sub,1)
        println("Selected $number_rows variants with Filter status: PASS and that match list of chromosome positions of interest")
    end

    #genomic_range and list
    if occursin("chr",genomic_range) && occursin(".csv",positions_list) && !occursin("true",pass_filter)
        sig_list =  load_siglist(positions_list)
        sub = VariantVisualization.genomic_range_siglist_filter(vcf_filename, sig_list, genomic_range)
        number_rows = size(sub,1)
        println("Selected $number_rows variants that are within chromosome range: $genomic_range and that match list of chromosome positions of interest")
    end

    #pass_filter
    if occursin("true",pass_filter) && !occursin("chr",genomic_range) && !occursin(".csv",positions_list)
        println("Only pass filter is applied. Large vcf files with many PASS variants will take a long time to process and heatmap visualizations will lose resolution at this scale unless viewed in interactive html for zooming.")
        println()
        sub = VariantVisualization.io_pass_filter(vcf_filename)
        number_rows = size(sub,1)
        println("Selected $number_rows variants with Filter status: PASS")
        heatmap_input = "pass_filtered"
    end

    #genomic_range
    if occursin("chr",genomic_range) && !occursin("true",pass_filter) && !occursin(".csv",positions_list)
        sub = VariantVisualization.io_genomic_range_vcf_filter(genomic_range,vcf_filename)
        number_rows = size(sub,1)
        println("Selected $number_rows variants within chromosome range: $genomic_range")
        heatmap_input = "range_filtered"
    end

    #list
    if occursin(".csv",positions_list) && !occursin("chr",genomic_range) && !occursin("true",pass_filter)
        sig_list =  load_siglist(positions_list)
        sub = VariantVisualization.io_sig_list_vcf_filter(sig_list,vcf_filename)
        number_rows = size(sub,1)
        println("Selected $number_rows variants that match list of chromosome positions of interest")
        heatmap_input = "positions_filtered"
    end

    println()
    println("Finished Filtering.")
    println("_______________________________________________")

    if occursin("genotype",heatmap_options) && !occursin("read_depth",heatmap_options)
        gt_num_array,gt_chromosome_labels = combined_all_genotype_array_functions(sub)

        if heatmap_title != ""
            title = "Genotype_$(heatmap_title)"
        else
            bn = Base.Filesystem.basename(vcf_filename)
            title = "Genotype_$bn"
        end

        chr_pos_tuple_list = generate_chromosome_positions_for_hover_labels(gt_chromosome_labels)

        chrom_label_info = VariantVisualization.chromosome_label_generator(gt_chromosome_labels[:,1])

        if length(group_samples) == 2

            if select_samples != ""

                println("Make sure that select_samples list and phenotype matrix contain the same sample names.")
                println("Selecting samples listed in $select_samples")

                gt_num_array,col_selectedcolumns = select_columns(select_samples,
                                              gt_num_array,
                                              sample_names)

                sample_names = col_selectedcolumns
            end

            group_trait_matrix_filename=group_samples[1]
            trait_to_group_by = group_samples[2]
            println()
            println("Grouping samples by $trait_to_group_by")
            println()

            ordered_num_array,group_label_pack,pheno,id_list,trait_labels = sortcols_by_phenotype_matrix(group_trait_matrix_filename, trait_to_group_by, gt_num_array, sample_names)

            if occursin("true",num_array)
            save_numerical_array(ordered_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
            end

            pheno_num_array,trait_label_array,chrom_label_info=add_pheno_matrix_to_gt_data_for_plotting(ordered_num_array,pheno,trait_labels,chrom_label_info,number_rows)

            graphic = VariantVisualization.genotype_heatmap_with_groups(pheno_num_array,title,chrom_label_info,group_label_pack,id_list,chr_pos_tuple_list,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)
            graphic

        else

            if select_samples != ""

                println("Selecting samples listed in $select_samples")

                gt_num_array,col_selectedcolumns = select_columns(select_samples,
                                              gt_num_array,
                                              sample_names)

                sample_names = col_selectedcolumns
            end

            if occursin("true",num_array)
            save_numerical_array(gt_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
            end

            #graphic = VariantVisualization.genotype_heatmap2(gt_num_array,title,chrom_label_info,sample_names,chr_pos_tuple_list,y_axis_label_option,x_axis_label_option)
            graphic = VariantVisualization.genotype_heatmap2_new_legend(gt_num_array,title,chrom_label_info,sample_names,chr_pos_tuple_list,y_axis_label_option,x_axis_label_option)
            graphic

        end

        println("Saving genotype heatmap")

        save_graphic(graphic,output_directory,save_ext,title,remote_option)

    end

    if occursin("read_depth",heatmap_options) && !occursin("genotype",heatmap_options)

        dp_num_array,dp_chromosome_labels = combined_all_read_depth_array_functions(sub)

        if heatmap_title != ""
            title = "Read_Depth_$heatmap_title"
        else
            bn = Base.Filesystem.basename(vcf_filename)
            title = "Read_Depth_$bn"
        end

        chr_pos_tuple_list = generate_chromosome_positions_for_hover_labels(dp_chromosome_labels)

        chrom_label_info = VariantVisualization.chromosome_label_generator(dp_chromosome_labels[:,1])

        if length(group_samples) == 2

            if select_samples != ""

                println("Make sure that select_samples list and phenotype matrix contain the same sample names.")
                println("Selecting samples listed in $select_samples")

                gt_num_array,col_selectedcolumns = select_columns(select_samples,
                                              dp_num_array,
                                              sample_names)

                sample_names = col_selectedcolumns
            end

            group_trait_matrix_filename=group_samples[1]
            trait_to_group_by = group_samples[2]
            println()
            println("Grouping samples by $trait_to_group_by")
            println()

            ordered_dp_num_array,group_label_pack,pheno,id_list,trait_labels = sortcols_by_phenotype_matrix(group_trait_matrix_filename, trait_to_group_by, dp_num_array, sample_names)
            dp_num_array_limited=read_depth_threshhold(ordered_dp_num_array)

             if occursin("true",num_array)
            save_numerical_array(ordered_dp_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
            end

            pheno_num_array,trait_label_array,chrom_label_info = add_pheno_matrix_to_dp_data_for_plotting(dp_num_array_limited,pheno,trait_labels,chrom_label_info,number_rows)

            graphic = VariantVisualization.dp_heatmap2_with_groups(pheno_num_array,title,chrom_label_info,group_label_pack,id_list,chr_pos_tuple_list,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)
            graphic

        else

                       if select_samples != ""

                           println("Selecting samples listed in $select_samples")

                           dp_num_array,col_selectedcolumns = select_columns(select_samples,
                                                         gt_num_array,
                                                         sample_names)

                           sample_names = col_selectedcolumns
                       end

            if occursin("true",num_array)
            save_numerical_array(dp_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
            end

            dp_num_array_limited=read_depth_threshhold(dp_num_array)

            graphic = VariantVisualization.dp_heatmap2(dp_num_array, title, chrom_label_info, sample_names,chr_pos_tuple_list,y_axis_label_option,x_axis_label_option)
           graphic

        end

        println("Saving read depth heatmap")

        save_graphic(graphic,output_directory,save_ext,title,remote_option)
    end

    if occursin("read_depth",heatmap_options) && occursin("genotype",heatmap_options)

        gt_num_array,gt_chromosome_labels = combined_all_genotype_array_functions(sub)

        if heatmap_title != ""
            title = "Genotype_$(heatmap_title)"
        else
            bn = Base.Filesystem.basename(vcf_filename)
            title = "Genotype_$bn"
        end

        chr_pos_tuple_list = generate_chromosome_positions_for_hover_labels(gt_chromosome_labels)

        chrom_label_info = VariantVisualization.chromosome_label_generator(gt_chromosome_labels[:,1])

        if length(group_samples) == 2

            if select_samples != ""

                println("Make sure that select_samples list and phenotype matrix contain the same sample names.")
                println("Selecting samples listed in $select_samples")

                gt_num_array,col_selectedcolumns = select_columns(select_samples,
                                              gt_num_array,
                                              sample_names)

                sample_names = col_selectedcolumns
            end

            group_trait_matrix_filename=group_samples[1]
            trait_to_group_by = group_samples[2]
            println()
            println("Grouping samples by $trait_to_group_by")
            println()

            ordered_num_array,group_label_pack,pheno,id_list,trait_labels = sortcols_by_phenotype_matrix(group_trait_matrix_filename, trait_to_group_by, gt_num_array, sample_names)

            if occursin("true",num_array)
            save_numerical_array(ordered_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
            end

            pheno_num_array,trait_label_array,chrom_label_info=add_pheno_matrix_to_gt_data_for_plotting(ordered_num_array,pheno,trait_labels,chrom_label_info,number_rows)

            #graphic = VariantVisualization.genotype_heatmap_with_groups(pheno_num_array,title,chrom_label_info,group_label_pack,id_list,chr_pos_tuple_list,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)
            graphic = VariantVisualization.genotype_heatmap_with_groups(pheno_num_array,title,chrom_label_info,group_label_pack,id_list,chr_pos_tuple_list,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)
            graphic

        else

            if select_samples != ""

                println("Selecting samples listed in $select_samples")

                gt_num_array,col_selectedcolumns = select_columns(select_samples,
                                              gt_num_array,
                                              sample_names)

                sample_names = col_selectedcolumns
            end

            if occursin("true",num_array)
            save_numerical_array(gt_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
            end

            #graphic = VariantVisualization.genotype_heatmap2(gt_num_array,title,chrom_label_info,sample_names,chr_pos_tuple_list,y_axis_label_option,x_axis_label_option)
            graphic = VariantVisualization.genotype_heatmap2_new_legend(gt_num_array,title,chrom_label_info,sample_names,chr_pos_tuple_list,y_axis_label_option,x_axis_label_option)

        end

        println("Saving genotype heatmap")

        save_graphic(graphic,output_directory,save_ext,title,remote_option)

        if heatmap_title != ""
            title = "Read_Depth_$heatmap_title"
        else
            bn = Base.Filesystem.basename(vcf_filename)
            title = "Read_Depth_$bn"
        end

        dp_num_array,dp_chromosome_labels = combined_all_read_depth_array_functions(sub)

        chr_pos_tuple_list = generate_chromosome_positions_for_hover_labels(dp_chromosome_labels)

        chrom_label_info = VariantVisualization.chromosome_label_generator(dp_chromosome_labels[:,1])

        if length(group_samples) == 2

            if select_samples != ""

                println("Make sure that select_samples list and phenotype matrix contain the same sample names.")
                println("Selecting samples listed in $select_samples")

                gt_num_array,col_selectedcolumns = select_columns(select_samples,
                                              dp_num_array,
                                              sample_names)

                sample_names = col_selectedcolumns
            end

            group_trait_matrix_filename=group_samples[1]
            trait_to_group_by = group_samples[2]
            println()
            println("Grouping samples by $trait_to_group_by")
            println()

            ordered_dp_num_array,group_label_pack,pheno,id_list,trait_labels = sortcols_by_phenotype_matrix(group_trait_matrix_filename, trait_to_group_by, dp_num_array, sample_names)
            dp_num_array_limited=read_depth_threshhold(ordered_dp_num_array)

            if occursin("true",num_array)
            save_numerical_array(ordered_dp_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
            end

            pheno_num_array,trait_label_array,chrom_label_info = add_pheno_matrix_to_dp_data_for_plotting(dp_num_array_limited,pheno,trait_labels,chrom_label_info,number_rows)

            graphic = VariantVisualization.dp_heatmap2_with_groups(pheno_num_array,title,chrom_label_info,group_label_pack,id_list,chr_pos_tuple_list,y_axis_label_option,trait_label_array,x_axis_label_option,number_rows)
            graphic

        else

           if select_samples != ""

               println("Selecting samples listed in $select_samples")

               dp_num_array,col_selectedcolumns = select_columns(select_samples,
                                             gt_num_array,
                                             sample_names)

               sample_names = col_selectedcolumns
           end

            if occursin("true",num_array)
            save_numerical_array(dp_num_array,sample_names,chr_pos_tuple_list,title,output_directory)
            end

            dp_num_array_limited=read_depth_threshhold(dp_num_array)

            graphic = VariantVisualization.dp_heatmap2(dp_num_array, title, chrom_label_info, sample_names,chr_pos_tuple_list,y_axis_label_option,x_axis_label_option)
            graphic

            end

        println("Saving read depth heatmap")

        save_graphic(graphic,output_directory,save_ext,title,remote_option)

    end

end


text:
module Cumulants
  using SymmetricTensors
  using Combinatorics
  using Distributions
  using Distributed
  import SymmetricTensors: pyramidindices, ind2range, sizetest, getblockunsafe
  import Distributions: moment
  if VERSION >= v"1.3"
   using CompilerSupportLibraries_jll
 end

  #calculates moments and cumulants using block structures (SymmetricTensors)
  include("cumulant.jl")

  #naive implementation
  include("naivecumulants.jl")

  export moment, cumulants, naivecumulant, naivemoment
end


text:
@testset "Objective" begin
    T = 3
    num_state = 2
    num_action = 1
    num_parameter = 0
    ot = (x, u, w) -> dot(x, x) + 0.1 * dot(u, u)
    oT = (x, u, w) -> 10.0 * dot(x, x)
    ct = Cost(ot, num_state, num_action, num_parameter=num_parameter)
    cT = Cost(oT, num_state, 0, num_parameter=num_parameter)
    objective = [[ct for t = 1:T-1]..., cT]

    J = [0.0]
    gradient = zeros((T - 1) * (num_state + num_action) + num_state)
    idx_xu = [collect((t - 1) * (num_state + num_action) .+ (1:(num_state + (t == T ? 0 : num_action)))) for t = 1:T]
    x1 = ones(num_state)
    u1 = ones(num_action)
    w1 = zeros(num_parameter)
    X = [x1 for t = 1:T]
    U = [t < T ? u1 : zeros(0) for t = 1:T]
    W = [w1 for t = 1:T]

    ct.evaluate(ct.evaluate_cache, x1, u1, w1)
    ct.gradient(ct.gradient_cache, x1, u1, w1)
    @test ct.evaluate_cache[1] ≈ ot(x1, u1, w1)
    @test norm(ct.gradient_cache - [2.0 * x1; 0.2 * u1]) < 1.0e-8

    cT.evaluate(cT.evaluate_cache, x1, u1, w1)
    cT.gradient(cT.gradient_cache, x1, u1, w1)
    @test cT.evaluate_cache[1] ≈ oT(x1, u1, w1)
    @test norm(cT.gradient_cache - 20.0 * x1) < 1.0e-8

    @test DTO.cost(objective, X, U, X) - sum([ot(X[t], U[t], W[t]) for t = 1:T-1]) - oT(X[T], U[T], W[T]) ≈ 0.0
    DTO.gradient!(gradient, idx_xu, objective, X, U, W)
    @test norm(gradient - vcat([[2.0 * x1; 0.2 * u1] for t = 1:T-1]..., 20.0 * x1)) < 1.0e-8

    # info = @benchmark DTO.cost($objective, $X, $U, $W)
    # info = @benchmark DTO.gradient!($gradient, $idx_xu, $objective, $X, $U, $W)
end


text:
function linreg_loss(w, b, x, y)
    yhat = sum(w .* x) + b
    loss = sum(abs2(yhat - y))
end

inputs = [:w => [1.0, 1.0], :b => 2.0, :x => [3.0, 3.0], :y => 4.0]
vals = ([v for (k,v) in inputs]...)
dlinreg_loss = xdiff(linreg_loss; inputs...)
dvals = dlinreg_loss(vals...)
@test dvals == (16.0, [96.0, 96.0], 32.0, [32.0, 32.0], -32.0)

# ReverseDiff doesn't support scalars yet
# See: https://github.com/JuliaDiff/ReverseDiff.jl/issues/40
# linreg_loss2(W, b, x, y) = sum(abs2.(y - (W * x .+ b)))
# test_compare(linreg_loss2; W=randn(1,11), b = rand(), x=randn(11,1000), y=randn(1,1000))