text:
mutable struct NetworkRecorder
	simulationtime::Array{Int64,1}
	liveedges::Array{Int64,1}
	livesources::Array{Int64,1}
	weightedlivesources::Array{Float64,1}
	livesinks::Array{Int64,1}
	totalpower::Array{Float64,1}
	flowsequence::Array{Array{Float64,1},1}
end

function newrecorder(numedges::Int64, ns::Int64, P::Float64, nd::Int64, F::Float64, flowvec::Vector{Float64})
	r = NetworkRecorder(Int64[],Int64[],Int64[],Float64[],Int64[],Float64[],Float64[])
	push!(r.simulationtime, 0)
	push!(r.liveedges, numedges)
	push!(r.livesources, ns)
	push!(r.weightedlivesources, P)
	push!(r.livesinks, nd)
	push!(r.totalpower, F)
	push!(r.flowsequence, flowvec)
	return r
end

mutable struct NetworkRecorder2
	iterations::Array{Int64,1}
	liveedges::Array{Int64,1}
	livesources::Array{Int64,1}
	livesinks::Array{Int64,1}
	totalpower::Array{Float64,1}
end

function cascaderesults(info::NetworkRecorder)
	timetaken = info.simulationtime[end]
	fracedges = Float64(info.liveedges[end])/Float64(info.liveedges[1])
	fracsources = Float64(info.livesources[end])/Float64(info.livesources[1])
	fracweightedsources = info.weightedlivesources[end]/info.weightedlivesources[1]
	fracsinks = Float64(info.livesinks[end])/Float64(info.livesinks[1])
	fracpower = info.totalpower[end]/info.totalpower[1]
	return timetaken, fracedges, fracsources, fracweightedsources, fracsinks, fracpower, info.flowsequence
end

mutable struct CascadeNetwork
	network::Graph
	E::Matrix{Float64}
	sources::Vector{Int64}
	sinks::Vector{Int64}
	initialflow::Vector{Float64}
end

mutable struct CascadeNetworkSwing
	network::Graph
	E::Matrix{Float64}
	sources::Vector{Int64}
	sinks::Vector{Int64}
	initialflow::Vector{Float64}
	initialnodestate::Vector{Float64}
end

mutable struct CascadeNetworkRand
	network::Graph
	E::Matrix{Float64}
	sources::Vector{Int64}
	sinks::Vector{Int64}
	initialflow::Vector{Float64}
	alphanoise::Vector{Float64}
end

function triggerhist(F::Vector{Float64})
	m = length(F)
	println(sum(abs.(F)))
	probs = Array{Float64, 1}(undef, m)
	for i=1:m
		probs[i] = abs(F[i])/sum(abs.(F))
	end
	return probs
end

function triggerindex(F::Vector{Float64})
	m = length(F)
	while true
		z = rand(1:m)
		if rand() < abs(F[z])/sum(abs.(F))
			return z
		end
	end
end

function cascadebisectionrandomtrigger(supensemblesize::Int64, subensemblesize::Int64, tol::Float64, n::Int64, ns::Int64, nd::Int64, q::Float64, k::Int64, Ptot::Float64)

	#Make an ensemble of networks here, with node placements
	#and with initial flow pattern
	ensemble = Array{CascadeNetwork, 1}(undef, supensemblesize*subensemblesize)
	fmaxes = Array{Float64, 1}(undef, supensemblesize*subensemblesize)
	degmeans = Array{Float64, 1}(undef, supensemblesize*subensemblesize)
	triggervec = Array{Int64, 1}(undef, supensemblesize*subensemblesize)
	for i=1:supensemblesize

		#Construct each network
		net, E = smallworldnetwork(n, k, q)
		# net, E = scalefreenetwork(n, 2*n, 4)
		sources, sinks = sourcesinklocs(ns, nd, n)

		#Make the initial flow
		m = numedges(net)
		X = zeros(Float64, m) .+ 1.0
		P = sourcesinkvector(sources, sinks, n, Ptot)
		L = weightedlaplacianmatrix(E, X)
		theta = nodevoltages(L, P)
		F = edgecurrents(theta, E, X)
		flowmaxtemp = maximum(abs.(F))
		F = F./flowmaxtemp

		for j=1:subensemblesize
			ind = (i-1)*subensemblesize + j
			ensemblemember = CascadeNetwork(net, E, sources, sinks, F)
			ensemble[ind] = ensemblemember
			triggervec[ind] = triggerindex(F)
			degs = getdegrees(E)
			degmeans[ind] = mean(degs)
			fmaxes[ind] = flowmaxtemp
		end
	end

	alpha = 0.01
	stepsize = 0.3
	beneath0 = true
	beneath1 = true
	iterations = 1
	ensemblesize = supensemblesize*subensemblesize
	while abs(stepsize) > tol
		println(alpha)

		#For every member in the ensemble
		runningsurvivors = Array{Float64, 1}(undef, ensemblesize)
		for i=1:ensemblesize
			# println(i)

			#Unpack the network
			net = ensemble[i].network
			E = ensemble[i].E
			m = numedges(net)
			X = zeros(Float64, m) .+ 1.0
			sources = ensemble[i].sources
			sinks = ensemble[i].sinks
			F = ensemble[i].initialflow
			nodeset = collect(1:n)
			edgeset = collect(1:m)
			P = sourcesinkvector(sources, sinks, n, Ptot)

			#Reactivate all edges
			for j=1:length(net.edgelist)
				net.edgelist[j].active = true
			end

			#Knock out trigger
			d = triggervec[i]
			edgeset1 = filter(x->x≠d, edgeset)
			E1 = E[1:end, 1:end .≠ d]
			net.edgelist[d].active = false

			#Check to see if its broken
			Adj = adjacencymatrix(E1)
			components, table = connectedcomponents(Adj)

			#Do cascade
			flowvec = zeros(Float64, length(F)) .+ F
			flowvec[d] = NaN
			info = newrecorder(m, ns, Ptot, nd, sum(abs.(P))/2.0, flowvec)
			simtime = 0
			for j=1:components
				chunk = table[j]
				nodesubset = zeros(Int64, length(chunk))
				for l=1:length(chunk)
					nodesubset[l] = nodeset[chunk[l]]
				end
				info = fracture!(net, nodesubset, P, X, F, info, simtime, alpha, fmaxes[i], false)
			end

			#Record fraction of edges remaining
			T, fracedges, fracsources, fracweightedsources, fracsinks, fracpower, flowseq = cascaderesults(info)
			runningsurvivors[i] = fracedges

		end

		#If fraction has crossed the target then reverse and half the step size
		global avremaining = mean(runningsurvivors)
		if avremaining > 0.8
			beneath1 = false
		else
			beneath1 = true
		end
		if beneath1 ≠ beneath0
			stepsize = -stepsize/2.0
		end
		beneath0 = beneath1
		alpha += stepsize
		iterations += 1
	end

	return alpha, mean(fmaxes), iterations, degmeans

end

function sequencedstartbisection(tol::Float64, n::Int64, ns::Int64, nd::Int64, q::Float64, k::Int64)

	#Construct a network
	net, E = smallworldnetwork(n, k, q)
	# net, E = scalefreenetwork(n, 2*n, 3)
	sources, sinks = sourcesinklocs(ns, nd, n)
	# println("done")

	#Make the initial flow
	m = numedges(net)
	# println(m)
	X = zeros(Float64, m) .+ 1.0
	P = sourcesinkvector(sources, sinks, n, 1.0)
	L = weightedlaplacianmatrix(E, X)
	theta = nodevoltages(L, P)
	F = edgecurrents(theta, E, X)
	fmax = maximum(abs.(F))
	F = F./fmax
	nodeset = collect(1:n)
	edgeset = collect(1:m)

	alphacs = zeros(Float64, m)
	relativeload = zeros(Float64, m)

	#Iterate through triggering edges
	for z=1:m

		#Start bisection
		alpha = 0.01
		stepsize = 0.3
		beneath0 = true
		beneath1 = true
		iterations = 1

		relativeload[z] = abs(F[z])/mean(abs.(F))

		while abs(stepsize) > tol

			#Reactivate all edges
			for i=1:m
				net.edgelist[i].active = true
			end

			#Knock out trigger edge
			edgeset1 = filter(x->x≠z, edgeset)
			E1 = E[1:end, 1:end .≠ z]
			net.edgelist[z].active = false

			#Check to see if its broken
			Adj = adjacencymatrix(E1)
			components, table = connectedcomponents(Adj)

			#Do cascade
			flowvec = zeros(Float64, length(F)) .+ F
			flowvec[z] = NaN
			info = newrecorder(m, ns, 1.0, nd, sum(abs.(P))/2.0, flowvec)
			simtime = 0
			for j=1:components
				chunk = table[j]
				nodesubset = zeros(Int64, length(chunk))
				for l=1:length(chunk)
					nodesubset[l] = nodeset[chunk[l]]
				end
				info = fracture!(net, nodesubset, P, X, F, info, simtime, alpha, fmax, false)
			end

			#Record fraction of edges remaining
			T, fracedges, fracsources, fracweightedsources, fracsinks, fracpower, flowseq = cascaderesults(info)

			#If fraction has crossed the target then reverse and half the step size
			if fracedges > 0.5
				beneath1 = false
			else
				beneath1 = true
			end
			if beneath1 ≠ beneath0
				stepsize = -stepsize/2.0
			end
			beneath0 = beneath1
			alpha += stepsize
			iterations += 1

		end
		alphacs[z] = alpha
	end
	return relativeload, alphacs, fmax
end

function cascadebisection(ensemblesize::Int64, tol::Float64, n::Int64, ns::Int64, nd::Int64, q::Float64, k::Int64, Ptot::Float64)

	#Make an ensemble of networks here, with node placements
	#and with initial flow pattern
	ensemble = Array{CascadeNetwork, 1}(undef, ensemblesize)
	fmaxes = Array{Float64, 1}(undef, ensemblesize)
	degmeans = Array{Float64, 1}(undef, ensemblesize)
	for i=1:ensemblesize

		#Construct each network
		net, E = smallworldnetwork(n, k, q)
		# net, E = scalefreenetwork(n, 2*n, 4)
		degs = getdegrees(E)
		degmeans[i] = mean(degs)
		sources, sinks = sourcesinklocs(ns, nd, n)

		#Make the initial flow
		m = numedges(net)
		X = zeros(Float64, m) .+ 1.0
		P = sourcesinkvector(sources, sinks, n, Ptot)
		L = weightedlaplacianmatrix(E, X)
		theta = nodevoltages(L, P)
		F = edgecurrents(theta, E, X)
		flowmaxtemp = maximum(abs.(F))
		fmaxes[i] = flowmaxtemp
		F = F./flowmaxtemp

		#Store the ensemble member
		ensemblemember = CascadeNetwork(net, E, sources, sinks, F)
		ensemble[i] = ensemblemember
		# println(i)
	end

	alpha = 0.01
	stepsize = 0.3
	beneath0 = true
	beneath1 = true
	iterations = 1
	while abs(stepsize) > tol
		# println(alpha)

		#For every member in the ensemble
		runningsurvivors = Array{Float64, 1}(undef, ensemblesize)
		for i=1:ensemblesize

			#Unpack the network
			net = ensemble[i].network
			E = ensemble[i].E
			m = numedges(net)
			X = zeros(Float64, m) .+ 1.0
			sources = ensemble[i].sources
			sinks = ensemble[i].sinks
			F = ensemble[i].initialflow
			nodeset = collect(1:n)
			edgeset = collect(1:m)
			P = sourcesinkvector(sources, sinks, n, Ptot)

			#Reactivate all edges
			for j=1:length(net.edgelist)
				net.edgelist[j].active = true
			end

			#Knock out largest
			d = argmax(abs.(F))
			edgeset1 = filter(x->x≠d, edgeset)
			E1 = E[1:end, 1:end .≠ d]
			net.edgelist[d].active = false

			#Check to see if its broken
			Adj = adjacencymatrix(E1)
			components, table = connectedcomponents(Adj)

			#Do cascade
			flowvec = zeros(Float64, length(F)) .+ F
			flowvec[d] = NaN
			info = newrecorder(m, ns, Ptot, nd, sum(abs.(P))/2.0, flowvec)
			simtime = 0
			for j=1:components
				chunk = table[j]
				nodesubset = zeros(Int64, length(chunk))
				for l=1:length(chunk)
					nodesubset[l] = nodeset[chunk[l]]
				end
				info = fracture!(net, nodesubset, P, X, F, info, simtime, alpha, fmaxes[i], false)
			end

			#Record fraction of edges remaining
			T, fracedges, fracsources, fracweightedsources, fracsinks, fracpower, flowseq = cascaderesults(info)
			runningsurvivors[i] = fracedges

		end

		#If fraction has crossed the target then reverse and half the step size
		global avremaining = mean(runningsurvivors)
		if avremaining > 0.5
			beneath1 = false
		else
			beneath1 = true
		end
		if beneath1 ≠ beneath0
			stepsize = -stepsize/2.0
		end
		beneath0 = beneath1
		alpha += stepsize
		iterations += 1
	end

	return alpha, mean(fmaxes), iterations, degmeans

end

function cascadebisection2(ensemblesize::Int64, tol::Float64, n::Int64, ns::Int64, nd::Int64, q::Float64, k::Int64, Ptot::Float64)

	#Make an ensemble of networks here, with node placements
	#and with initial flow pattern
	ensemble = Array{CascadeNetwork, 1}(undef, ensemblesize)
	fmaxes = Array{Float64, 1}(undef, ensemblesize)
	degmeans = Array{Float64, 1}(undef, ensemblesize)
	for i=1:ensemblesize

		#Construct each network
		net, E = smallworldnetwork(n, k, q)
		# net, E = scalefreenetwork(n, 2*n, 3)
		degs = getdegrees(E)
		degmeans[i] = mean(degs)
		sources, sinks = sourcesinklocs(ns, nd, n)

		#Make the initial flow
		m = numedges(net)
		X = zeros(Float64, m) .+ 1.0
		P = sourcesinkvector(sources, sinks, n, Ptot)
		L = weightedlaplacianmatrix(E, X)
		theta = nodevoltages(L, P)
		F = edgecurrents(theta, E, X)
		flowmaxtemp = maximum(abs.(F))
		fmaxes[i] = flowmaxtemp
		F = F./flowmaxtemp

		#Store the ensemble member
		ensemblemember = CascadeNetwork(net, E, sources, sinks, F)
		ensemble[i] = ensemblemember
		# println(i)
	end

	alpha = 0.01
	stepsize = 0.3
	beneath0 = true
	beneath1 = true
	iterations = 1
	while abs(stepsize) > tol
		# println(alpha)

		#For every member in the ensemble
		runningsurvivors = Array{Float64, 1}(undef, ensemblesize)
		for i=1:ensemblesize

			#Unpack the network
			net = ensemble[i].network
			E = ensemble[i].E
			m = numedges(net)
			X = zeros(Float64, m) .+ 1.0
			sources = ensemble[i].sources
			sinks = ensemble[i].sinks
			F = ensemble[i].initialflow
			nodeset = collect(1:n)
			edgeset = collect(1:m)
			P = sourcesinkvector(sources, sinks, n, Ptot)

			#Reactivate all edges
			for j=1:length(net.edgelist)
				net.edgelist[j].active = true
			end

			#Knock out largest
			d = argmax(abs.(F))
			edgeset1 = filter(x->x≠d, edgeset)
			E1 = E[1:end, 1:end .≠ d]
			net.edgelist[d].active = false

			#Check to see if its broken
			Adj = adjacencymatrix(E1)
			components, table = connectedcomponents(Adj)

			#Do cascade
			flowvec = zeros(Float64, length(F)) .+ F
			flowvec[d] = NaN
			info = newrecorder(m, ns, Ptot, nd, sum(abs.(P))/2.0, flowvec)
			simtime = 0
			remainingedges = 0
			for j=1:components
				chunk = table[j]
				nodesubset = zeros(Int64, length(chunk))
				for l=1:length(chunk)
					nodesubset[l] = nodeset[chunk[l]]
				end
				remainingedges += fracture2!(net, nodesubset, P, X, F, simtime, alpha, fmaxes[i])
			end

			#Record fraction of edges remaining
			runningsurvivors[i] = remainingedges/Float64(m)

		end

		#If fraction has crossed the target then reverse and half the step size
		global avremaining = mean(runningsurvivors)
		if avremaining > 0.5
			beneath1 = false
		else
			beneath1 = true
		end
		if beneath1 ≠ beneath0
			stepsize = -stepsize/2.0
		end
		beneath0 = beneath1
		alpha += stepsize
		iterations += 1
	end

	return alpha, mean(fmaxes), iterations, degmeans
end

function cascadebisection3rand(ensemblesize::Int64, tol::Float64, n::Int64, ns::Int64, nd::Int64, q::Float64, k::Int64, Ptot::Float64)

	#Make an ensemble of networks here, with node placements
	#and with initial flow pattern
	ensemble = Array{CascadeNetworkRand, 1}(undef, ensemblesize)
	fmaxes = Array{Float64, 1}(undef, ensemblesize)
	degmeans = Array{Float64, 1}(undef, ensemblesize)
	for i=1:ensemblesize

		#Construct each network
		net, E = smallworldnetwork(n, k, q)
		# net, E = scalefreenetwork(n, 2*n, 3)
		degs = getdegrees(E)
		degmeans[i] = mean(degs)
		sources, sinks = sourcesinklocs(ns, nd, n)

		#Make the initial flow
		m = numedges(net)
		X = zeros(Float64, m) .+ 1.0
		P = sourcesinkvector(sources, sinks, n, Ptot)
		L = weightedlaplacianmatrix(E, X)
		theta = nodevoltages(L, P)
		F = edgecurrents(theta, E, X)
		flowmaxtemp = maximum(abs.(F))
		fmaxes[i] = flowmaxtemp
		F = F./flowmaxtemp

		#Make noise for the capacities
		d = Normal(0.0, 0.05)
		alphanoise = rand(d, m)

		#Store the ensemble member
		ensemblemember = CascadeNetworkRand(net, E, sources, sinks, F, alphanoise)
		ensemble[i] = ensemblemember
		# println(i)
	end

	alpha = 0.01
	stepsize = 0.3
	beneath0 = true
	beneath1 = true
	iterations = 1
	while abs(stepsize) > tol
		# println(alpha)

		#For every member in the ensemble
		runningsurvivors = Array{Float64, 1}(undef, ensemblesize)
		for i=1:ensemblesize

			#Unpack the network
			net = ensemble[i].network
			E = ensemble[i].E
			capnoise = ensemble[i].alphanoise
			m = numedges(net)
			X = zeros(Float64, m) .+ 1.0
			sources = ensemble[i].sources
			sinks = ensemble[i].sinks
			F = ensemble[i].initialflow
			nodeset = collect(1:n)
			edgeset = collect(1:m)
			P = sourcesinkvector(sources, sinks, n, Ptot)

			#Reactivate all edges
			for j=1:length(net.edgelist)
				net.edgelist[j].active = true
			end

			#Knock out largest
			d = argmax(abs.(F))
			edgeset1 = filter(x->x≠d, edgeset)
			E1 = E[1:end, 1:end .≠ d]
			net.edgelist[d].active = false

			#Check to see if its broken
			Adj = adjacencymatrix(E1)
			components, table = connectedcomponents(Adj)

			#Do cascade
			flowvec = zeros(Float64, length(F)) .+ F
			flowvec[d] = NaN
			info = newrecorder(m, ns, Ptot, nd, sum(abs.(P))/2.0, flowvec)
			simtime = 0
			remainingedges = 0
			for j=1:components
				chunk = table[j]
				nodesubset = zeros(Int64, length(chunk))
				for l=1:length(chunk)
					nodesubset[l] = nodeset[chunk[l]]
				end
				remainingedges += fracturerand!(net, nodesubset, P, X, F, capnoise, simtime, alpha, fmaxes[i])
			end

			#Record fraction of edges remaining
			runningsurvivors[i] = remainingedges/Float64(m)

		end

		#If fraction has crossed the target then reverse and half the step size
		global avremaining = mean(runningsurvivors)
		if avremaining > 0.5
			beneath1 = false
		else
			beneath1 = true
		end
		if beneath1 ≠ beneath0
			stepsize = -stepsize/2.0
		end
		beneath0 = beneath1
		alpha += stepsize
		iterations += 1
	end

	return alpha, mean(fmaxes), iterations, degmeans
end

function smcascadefixedtop(ensemblesize::Int64, n::Int64, ns::Int64, nd::Int64, q::Float64, alphares::Int64, alphamin::Float64, alphamax::Float64, k::Int64, Ptot::Float64, Xval::Float64)
	#A single network with an ensemble of different source-sink locations

	#Create a network topology
	net, E0 =  austriannetwork() #smallworldnetwork(n, k, q)
	m = numedges(net)
	nodeset = collect(1:n)
	edgeset = collect(1:m)
	X = zeros(Float64, m) .+ Xval

	alphas = collect(range(alphamin, length=alphares, stop=alphamax))
	edgemat = zeros(Float64, ensemblesize, alphares)
	Tmat = zeros(Float64, ensemblesize, alphares)

	#Loop over ensemble
	flowmax = zeros(Float64, ensemblesize)
	for z=1:ensemblesize
		# println(z)

		#Create a P vector
		sources, sinks = sourcesinklocsmod(net, ns, nd, n)
		P = noisysourcesinkvector(sources, sinks, n, Ptot)

		#Solve the DC flow equations
		L = weightedlaplacianmatrix(E0, 1.0./X)
		theta = nodevoltages(L, P)
		F = edgecurrents(theta, E0, X)
		flowmaxtemp = maximum(abs.(F))
		F = F./flowmaxtemp
		flowmax[z] = flowmaxtemp

		#Empty vectors to record data
		rT = zeros(Float64, alphares)
		redges = zeros(Float64, alphares)
		rsources = zeros(Float64, alphares)
		rsinks = zeros(Float64, alphares)
		rweightedsources = zeros(Float64, alphares)
		rpower = zeros(Float64, alphares)

		#Loop over alpha
		# println(z)
		for (ind, alpha) ∈ enumerate(alphas)
			# println(ind)

			#Reactive all the edges
			for i=1:length(net.edgelist)
				net.edgelist[i].active = true
			end

			#Knock out most loaded edge
			d = argmax(abs.(F))
			edgeset1 = filter(x->x≠d, edgeset)
			E1 = E0[1:end, 1:end .≠ d]
			net.edgelist[d].active = false

			#Check to see if its broken
			Adj = adjacencymatrix(E1)
			components, table = connectedcomponents(Adj)

			flowvec = zeros(Float64, length(F)) .+ F
			flowvec[d] = NaN
			info = newrecorder(m, ns, Ptot, nd, sum(abs.(P))/2.0, flowvec)
			simtime = 0

			for i=1:components
				chunk = table[i]
				nodesubset = zeros(Int64, length(chunk))
				for j=1:length(chunk)
					nodesubset[j] = nodeset[chunk[j]]
				end
				info = fracture!(net, nodesubset, P, X, F, info, simtime, alpha, flowmaxtemp, false)
			end

			#Build up vectors of rT etc
			T, fracedges, fracsources, fracweightedsources, fracsinks, fracpower, flowseq = cascaderesults(info)
			rT[ind] = T
			redges[ind] = fracedges
			rsources[ind] = fracsources
			rsinks[ind] = fracsinks
			rweightedsources[ind] = fracweightedsources
			rpower[ind] = fracpower
		end

		#Stack up vecs of rT etc.
		edgemat[z,:] = redges
		Tmat[z,:] = rpower
	end

	# p0 = 10
	# p1 = 15
	# p2 = 75

	#return the mean of the rT etc data
	return alphas, mean(edgemat, dims=1)[1,:], std(edgemat, dims=1)[1,:], mean(Tmat, dims=1)[1,:], std(Tmat, dims=1)[1,:], edgemat, Tmat, mean(flowmax)
end

function simplecascade(n::Int64, ns::Int64, nd::Int64, q::Float64, alpha::Float64, k::Int64, Ptot::Float64, Xval::Float64)

	#Create a network topology
	net, E0 = smallworldnetwork(n, k, q)
	m = numedges(net)
	nodeset = collect(1:n)
	edgeset = collect(1:m)
	X = zeros(Float64, m) .+ Xval

	#Create a P vector
	sources, sinks = sourcesinklocs(ns, nd, n)
	P = sourcesinkvector(sources, sinks, n, Ptot*n)

	#Solve the DC flow equations
	L = weightedlaplacianmatrix(E0, 1.0./X)
	theta = nodevoltages(L, P)
	F = edgecurrents(theta, E0, X)
	flowmaxtemp = maximum(abs.(F))
	F = F./flowmaxtemp
	meanflow = mean(abs.(F))

	#Reactive all the edges
	for i=1:length(net.edgelist)
		net.edgelist[i].active = true
	end

	d = argmax(abs.(F))
	edgeset1 = filter(x->x≠d, edgeset)
	E1 = E0[1:end, 1:end .≠ d]
	net.edgelist[d].active = false

	#Check to see if its broken
	Adj = adjacencymatrix(E1)
	components, table = connectedcomponents(Adj)

	fragments = 0
	for i=1:components
		chunk = table[i]
		nodesubset = zeros(Int64, length(chunk))
		for j=1:length(chunk)
			nodesubset[j] = nodeset[chunk[j]]
		end
		fragments = fracture3!(net, nodesubset, P, X, F, 0, alpha, flowmaxtemp)
		totalfragments += fragments
	end

	return totalfragments
end

function smcascadedist(ensemblesize::Int64, n::Int64, ns::Int64, nd::Int64, q::Float64, alpha::Float64, k::Int64, Ptot::Float64, Xval::Float64)
	#A single network with an ensemble of different source-sink locations

	#Create a network topology
	net, E0 = smallworldnetwork(n, k, q)
	m = numedges(net)
	nodeset = collect(1:n)
	edgeset = collect(1:m)
	X = zeros(Float64, m) .+ Xval

	#Create a P vector
	sources, sinks = sourcesinklocs(ns, nd, n)
	P = sourcesinkvector(sources, sinks, n, Ptot*n)

	#Solve the DC flow equations
	L = weightedlaplacianmatrix(E0, 1.0./X)
	theta = nodevoltages(L, P)
	F = edgecurrents(theta, E0, X)
	flowmaxtemp = maximum(abs.(F))
	F = F./flowmaxtemp
	meanflow = mean(abs.(F))

	relativeload = zeros(Float64, ensemblesize)
	survivingedges = zeros(Float64, ensemblesize)
	survivingpower = zeros(Float64, ensemblesize)

	#Loop over ensemble
	for z=1:ensemblesize
		# println(z)

		#Reactive all the edges
		for i=1:length(net.edgelist)
			net.edgelist[i].active = true
		end

		#Knock out most loaded edge
		d = rand(1:m)
		# d = argmax(abs.(F))
		edgeset1 = filter(x->x≠d, edgeset)
		E1 = E0[1:end, 1:end .≠ d]
		net.edgelist[d].active = false

		#Compute and save relative load of knocked out edge
		relativeload[z] = abs(F[d])/meanflow

		#Check to see if its broken
		Adj = adjacencymatrix(E1)
		components, table = connectedcomponents(Adj)

		flowvec = zeros(Float64, length(F)) .+ F
		flowvec[d] = NaN
		info = newrecorder(m, ns, Ptot, nd, sum(abs.(P))/2.0, flowvec)
		simtime = 0

		for i=1:components
			chunk = table[i]
			nodesubset = zeros(Int64, length(chunk))
			for j=1:length(chunk)
				nodesubset[j] = nodeset[chunk[j]]
			end
			info = fracture!(net, nodesubset, P, X, F, info, simtime, alpha, flowmaxtemp, false)
		end

		#Build up vectors of rT etc
		T, fracedges, fracsources, fracweightedsources, fracsinks, fracpower, flowseq = cascaderesults(info)

		survivingedges[z] = fracedges
		survivingpower[z] = fracpower

	end

	#return the mean of the rT etc data
	return relativeload, survivingedges, survivingpower
end

function smcascadeclock(n::Int64, pos::Int64, q::Float64, alphares::Int64, alphamin::Float64, alphamax::Float64, k::Int64, Ptot::Float64, Xval::Float64)
	#A single network with with a fixed source-sink pair location
	#Create a network topology

	ensemblesize = 1
	net, E0 = smallworldnetwork(n, k, q)
	m = numedges(net)
	nodeset = collect(1:n)
	edgeset = collect(1:m)
	X = zeros(Float64, m) .+ Xval

	alphas = collect(range(alphamin, length=alphares, stop=alphamax))
	edgemat = zeros(Float64, ensemblesize, alphares)
	Tmat = zeros(Float64, ensemblesize, alphares)

	#Loop over ensemble
	for z=1:ensemblesize

		#Create a P vector
		sources, sinks = sourcessinkclockface(n, pos)
		P = sourcesinkvector(sources, sinks, n, Ptot)

		#Solve the DC flow equations
		L = weightedlaplacianmatrix(E0, 1.0./X)
		theta = nodevoltages(L, P)
		F = edgecurrents(theta, E0, X)
		flowmaxtemp = maximum(abs.(F))
		F = F./flowmaxtemp
		meanflow = mean(abs.(F))

		#Empty vectors to record data
		rT = zeros(Float64, alphares)
		redges = zeros(Float64, alphares)
		rsources = zeros(Float64, alphares)
		rsinks = zeros(Float64, alphares)
		rweightedsources = zeros(Float64, alphares)
		rpower = zeros(Float64, alphares)

		#Loop over alpha
		for (ind, alpha) ∈ enumerate(alphas)

			#Reactive all the edges
			for i=1:length(net.edgelist)
				net.edgelist[i].active = true
			end

			#Knock out most loaded edge
			d = argmax(abs.(F))
			edgeset1 = filter(x->x≠d, edgeset)
			E1 = E0[1:end, 1:end .≠ d]
			net.edgelist[d].active = false

			#Check to see if its broken
			Adj = adjacencymatrix(E1)
			components, table = connectedcomponents(Adj)

			flowvec = zeros(Float64, length(F)) .+ abs.(F)
			flowvec[d] = NaN
			info = newrecorder(m, 2, Ptot, 62, sum(abs.(P))/2.0, flowvec)
			simtime = 0

			for i=1:components
				chunk = table[i]
				nodesubset = zeros(Int64, length(chunk))
				for j=1:length(chunk)
					nodesubset[j] = nodeset[chunk[j]]
				end
				info = fracture!(net, nodesubset, P, X, F, info, simtime, alpha, flowmaxtemp, false)
			end

			#Build up vectors of rT etc
			T, fracedges, fracsources, fracweightedsources, fracsinks, fracpower, flowseq = cascaderesults(info)
			rT[ind] = T
			redges[ind] = fracedges
			rsources[ind] = fracsources
			rsinks[ind] = fracsinks
			rweightedsources[ind] = fracweightedsources
			rpower[ind] = fracpower
		end

		#Stack up vecs of rT etc.
		edgemat[z,:] = redges
		Tmat[z,:] = rT
	end

	#return the mean of the rT etc data
	return alphas, mean(edgemat, dims=1)[1,:], std(edgemat, dims=1)[1,:], mean(Tmat, dims=1)[1,:], std(Tmat, dims=1)[1,:], edgemat
end

function smcascadeclock2(n::Int64, pos::Int64, q::Float64, alphares::Int64, alphamin::Float64, alphamax::Float64, k::Int64, Ptot::Float64, Xval::Float64)
	#A single network with with a fixed source-sink pair location
	#Create a network topology

	ensemblesize = 1
	net, E0 = smallworldnetwork(n, k, q)
	m = numedges(net)
	nodeset = collect(1:n)
	edgeset = collect(1:m)
	X = zeros(Float64, m) .+ Xval

	alphas = collect(range(alphamin, length=alphares, stop=alphamax))
	edgemat = zeros(Float64, alphares)

	#Create a P vector
	sources, sinks = sourcessinkclockface(n, pos)
	P = sourcesinkvector(sources, sinks, n, Ptot)

	#Solve the DC flow equations
	L = weightedlaplacianmatrix(E0, 1.0./X)
	theta = nodevoltages(L, P)
	F = edgecurrents(theta, E0, X)
	flowmaxtemp = maximum(abs.(F))
	F = F./flowmaxtemp
	meanflow = mean(abs.(F))

	#Loop over alpha
	for (ind, alpha) ∈ enumerate(alphas)

		#Reactive all the edges
		for i=1:length(net.edgelist)
			net.edgelist[i].active = true
		end

		#Knock out most loaded edge
		d = argmax(abs.(F))
		edgeset1 = filter(x->x≠d, edgeset)
		E1 = E0[1:end, 1:end .≠ d]
		net.edgelist[d].active = false

		#Check to see if its broken
		Adj = adjacencymatrix(E1)
		components, table = connectedcomponents(Adj)

		flowvec = zeros(Float64, length(F)) .+ abs.(F)
		simtime = 0

		for i=1:components
			chunk = table[i]
			nodesubset = zeros(Int64, length(chunk))
			for j=1:length(chunk)
				nodesubset[j] = nodeset[chunk[j]]
			end
			edgemat[ind] = fracture2!(net, nodesubset, P, X, F, simtime, alpha, flowmaxtemp)
		end
	end

	#return the mean of the rT etc data
	return edgemat
end

function smcascadeclockrand(n::Int64, pos::Int64, q::Float64, alphares::Int64, alphamin::Float64, alphamax::Float64, k::Int64, Ptot::Float64, Xval::Float64)
	#A single network with with a fixed source-sink pair location
	#Create a network topology

	ensemblesize = 1
	net, E0 = smallworldnetwork(n, k, q)
	m = numedges(net)
	nodeset = collect(1:n)
	edgeset = collect(1:m)
	X = zeros(Float64, m) .+ Xval

	alphas = collect(range(alphamin, length=alphares, stop=alphamax))
	edgemat = zeros(Float64, alphares)

	#Make edge noise
	d = Normal(0.0, 0.05)
	alphanoise = rand(d, m)

	#Create a P vector
	sources, sinks = sourcessinkclockface(n, pos)
	P = sourcesinkvector(sources, sinks, n, Ptot)

	#Solve the DC flow equations
	L = weightedlaplacianmatrix(E0, 1.0./X)
	theta = nodevoltages(L, P)
	F = edgecurrents(theta, E0, X)
	flowmaxtemp = maximum(abs.(F))
	F = F./flowmaxtemp
	meanflow = mean(abs.(F))

	#Loop over alpha
	for (ind, alpha) ∈ enumerate(alphas)

		#Reactive all the edges
		for i=1:length(net.edgelist)
			net.edgelist[i].active = true
		end

		#Knock out most loaded edge
		d = argmax(abs.(F))
		edgeset1 = filter(x->x≠d, edgeset)
		E1 = E0[1:end, 1:end .≠ d]
		net.edgelist[d].active = false

		#Check to see if its broken
		Adj = adjacencymatrix(E1)
		components, table = connectedcomponents(Adj)

		flowvec = zeros(Float64, length(F)) .+ abs.(F)
		simtime = 0

		for i=1:components
			chunk = table[i]
			nodesubset = zeros(Int64, length(chunk))
			for j=1:length(chunk)
				nodesubset[j] = nodeset[chunk[j]]
			end
			edgemat[ind] = fracturerand!(net, nodesubset, P, X, F, alphanoise, simtime, alpha, flowmaxtemp)
		end
	end

	#return the mean of the rT etc data
	return edgemat
end

# function cascadesequence(net::Graph, E::Matrix{Float64}, sources::Vector{Int64}, sinks::Vector{Int64}, alpha::Float64, Ptot::Float64, Xval::Float64)

# 	n = size(E,1)
# 	m = size(E,2)
# 	nodeset = collect(1:n)
# 	edgeset = collect(1:m)
# 	X = zeros(Float64, m) .+ Xval
# 	ns = length(sources)
# 	nd = length(sinks)
# 	P = sourcesinkvector(sources, sinks, n, Ptot*n)

# 	#Solve the DC flow equations
# 	L = weightedlaplacianmatrix(E, 1.0./X)
# 	theta = nodevoltages(L, P)
# 	F = edgecurrents(theta, E, X)

# 	#Active all the edges
# 	for i=1:m
# 		net.edgelist[i].active = true
# 	end

# 	#Knock out most loaded edge
# 	d = argmax(abs.(F))
# 	edgeset1 = filter(x->x≠d, edgeset)
# 	E1 = E[1:end, 1:end .≠ d]
# 	net.edgelist[d].active = false

# 	#Check to see if its broken
# 	Adj = adjacencymatrix(E1)
# 	components, table = connectedcomponents(Adj)

# 	flowvec = zeros(Float64, length(F)) .+ abs.(F)
# 	println("Initial flow:")
# 	println(flowvec)
# 	flowvec[d] = NaN
# 	info = newrecorder(m, ns, Ptot, nd, sum(abs.(P))/2.0, flowvec)
# 	simtime = 0

# 	for i=1:components
# 		chunk = table[i]
# 		nodesubset = zeros(Int64, length(chunk))
# 		for j=1:length(chunk)
# 			nodesubset[j] = nodeset[chunk[j]]
# 		end
# 		info = fracture!(net, nodesubset, P, X, F, info, simtime, alpha, true)
# 	end

# 	T, fracedges, fracsources, fracweightedsources, fracsinks, fracpower, flowseq = cascaderesults(info)
# 	return T, fracedges, flowseq
# end

# function singlecascade(net::Graph, E0::Matrix{Float64}, ensemblesize::Int64, n::Int64, ns::Int64, nd::Int64, q::Float64, alphares::Int64, alphamin::Float64, alphamax::Float64, k::Int64, Ptot::Float64, Xval::Float64)
# 	#A single network with an ensemble of different source-sink locations

# 	#Create a network topology
# 	m = numedges(net)
# 	nodeset = collect(1:n)
# 	edgeset = collect(1:m)
# 	X = zeros(Float64, m) .+ Xval

# 	alphas = collect(range(alphamin, length=alphares, stop=alphamax))
# 	edgemat = zeros(Float64, ensemblesize, alphares)
# 	Tmat = zeros(Float64, ensemblesize, alphares)

# 	#Loop over ensemble
# 	flowmax = zeros(Float64, ensemblesize)
# 	for z=1:ensemblesize
# 		# println(z)

# 		#Create a P vector
# 		sources, sinks = sourcesinklocs(ns, nd, n)
# 		P = sourcesinkvector(sources, sinks, n, Ptot*n)

# 		#Solve the DC flow equations
# 		L = weightedlaplacianmatrix(E0, 1.0./X)
# 		theta = nodevoltages(L, P)
# 		F = edgecurrents(theta, E0, X)
# 		flowmax[z] = maximum(abs.(F))

# 		#Empty vectors to record data
# 		rT = zeros(Float64, alphares)
# 		redges = zeros(Float64, alphares)
# 		rsources = zeros(Float64, alphares)
# 		rsinks = zeros(Float64, alphares)
# 		rweightedsources = zeros(Float64, alphares)
# 		rpower = zeros(Float64, alphares)

# 		#Loop over alpha
# 		# println(z)
# 		for (ind, alpha) ∈ enumerate(alphas)
# 			# println(ind)

# 			#Reactive all the edges
# 			for i=1:length(net.edgelist)
# 				net.edgelist[i].active = true
# 			end

# 			#Knock out most loaded edge
# 			d = argmax(abs.(F))
# 			edgeset1 = filter(x->x≠d, edgeset)
# 			E1 = E0[1:end, 1:end .≠ d]
# 			net.edgelist[d].active = false

# 			#Check to see if its broken
# 			Adj = adjacencymatrix(E1)
# 			components, table = connectedcomponents(Adj)

# 			flowvec = zeros(Float64, length(F)) .+ F
# 			flowvec[d] = NaN
# 			info = newrecorder(m, ns, Ptot, nd, sum(abs.(P))/2.0, flowvec)
# 			simtime = 0

# 			for i=1:components
# 				chunk = table[i]
# 				nodesubset = zeros(Int64, length(chunk))
# 				for j=1:length(chunk)
# 					nodesubset[j] = nodeset[chunk[j]]
# 				end
# 				info = fracture!(net, nodesubset, P, X, F, info, simtime, alpha, false)
# 			end

# 			#Build up vectors of rT etc
# 			T, fracedges, fracsources, fracweightedsources, fracsinks, fracpower, flowseq = cascaderesults(info)
# 			rT[ind] = T
# 			redges[ind] = fracedges
# 			rsources[ind] = fracsources
# 			rsinks[ind] = fracsinks
# 			rweightedsources[ind] = fracweightedsources
# 			rpower[ind] = fracpower
# 		end

# 		#Stack up vecs of rT etc.
# 		edgemat[z,:] = redges
# 		Tmat[z,:] = rT
# 	end

# 	# p0 = 10
# 	# p1 = 15
# 	# p2 = 75

# 	#return the mean of the rT etc data
# 	return alphas, mean(edgemat, dims=1)[1,:], std(edgemat, dims=1)[1,:], mean(Tmat, dims=1)[1,:], std(Tmat, dims=1)[1,:], edgemat, Tmat, mean(flowmax)
# end

max_stars_repo_name: o-smith/GridResilience

max_issues_repo_name: o-smith/GridResilience

max_forks_repo_name: o-smith/GridResilience


text:
# test of code to remove redundante cycles after extracting quartet
# Claudia May 2015

include("/Users/Clauberry/Documents/phylo/software/CFimplementation/julia/git_laptop/CFnetworks/types.jl")
include("/Users/Clauberry/Documents/phylo/software/CFimplementation/julia/git_laptop/CFnetworks/functions.jl")

tree = "(((((((1,2),3),4),5),(6,7)),(8,9)),10);"

seed = 2738
currT0 = readTopologyUpdate(tree);
srand(seed)
besttree = deepcopy(currT0);
success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(besttree);
success
success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(besttree);
success
printEdges(besttree)
writeTopologyLevel1(besttree,true)
net = deepcopy(besttree);

q1 = Quartet(1,["1","2","3","4"],[0.5,0.4,0.1]);
extractQuartet!(net,q1);
qnet=deepcopy(q1.qnet);
printEdges(qnet)
printNodes(qnet)

## redundantCycle!(qnet,qnet.node[14])
## redundantCycle!(qnet,qnet.node[8])
## printEdges(qnet)
## printNodes(qnet)
identifyQuartet!(qnet)
eliminateHybridization!(qnet)
updateSplit!(qnet)
updateFormula!(qnet)
calculateExpCF!(qnet)

ext: jl


text:
abstract type AbstractScenarioProblems{S <: AbstractScenario} end

abstract type AbstractBlockStructure{N} <: AbstractStochasticStructure{N} end

# MOI #
# ========================== #
function MOI.is_valid(structure::AbstractBlockStructure, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer)
    stage == 1 && error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return MOI.is_valid(scenarioproblems(structure, stage), index, scenario_index)
end

function MOI.delete(structure::AbstractBlockStructure{N}, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    MOI.delete(scenarioproblems(structure, stage), index, scenario_index)
    return nothing
end
function MOI.delete(structure::AbstractBlockStructure{N}, indices::Vector{MOI.VariableIndex}, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    MOI.delete(scenarioproblems(structure, stage), indices, scenario_index)
    return nothing
end

# JuMP #
# ========================== #
function JuMP.objective_function_type(structure::AbstractBlockStructure{N}, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return objective_function_type(scenarioproblems(structure, stage), scenario_index)
end

function JuMP.objective_function(structure::AbstractBlockStructure{N},
                                 proxy::JuMP.Model,
                                 stage::Integer,
                                 scenario_index::Integer,
                                 FunType::Type{<:AbstractJuMPScalar}) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    objective::FunType = objective_function(scenarioproblems(structure, stage), proxy, scenario_index, FunType)
    return objective
end

function JuMP._moi_optimizer_index(structure::AbstractBlockStructure, index::VI, scenario_index::Integer)
    return JuMP._moi_optimizer_index(scenarioproblems(structure), index, scenario_index)
end
function JuMP._moi_optimizer_index(structure::AbstractBlockStructure, ci::CI, scenario_index::Integer)
    return JuMP._moi_optimizer_index(scenarioproblems(structure), ci, scenario_index)
end

function DecisionRef(proxy::JuMP.Model, structure::AbstractBlockStructure, index::VI, scenario_index::Integer)
    return DecisionRef(proxy, index)
end

# Getters #
# ========================== #
function scenarioproblems(structure::AbstractBlockStructure{N}, stage::Integer) where N
    1 < stage <= N || error("Stage $stage not in range 2 to $N.")
    stage == 1 && error("Stage 1 does not have scenario problems.")
    N == 2 && (stage == 2 || error("Stage $stage not available in two-stage model."))
    return structure.scenarioproblems[stage-1]
end
function scenarioproblems(structure::AbstractBlockStructure{2})
    return scenarioproblems(structure, 2)
end
function scenario_types(structure::AbstractBlockStructure{N}) where N
    return ntuple(Val{N-1}()) do s
        scenario_type(scenarioproblems(structure), s + 1)
    end
end
function decision(structure::AbstractBlockStructure{N}, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer) where N
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return decision(scenarioproblems(structure, stage), index, scenario_index)
end
function scenario(structure::AbstractBlockStructure, stage::Integer, scenario_index::Integer)
    scenario(scenarioproblems(structure, stage), scenario_index)
end
function scenario(structure::AbstractBlockStructure{2}, scenario_index::Integer)
    scenario(scenarioproblems(structure, 2), scenario_index)
end
function scenarios(structure::AbstractBlockStructure, stage::Integer)
    scenarios(scenarioproblems(structure, stage))
end
function scenarios(structure::AbstractBlockStructure{2})
    scenarios(scenarioproblems(structure, 2))
end
function expected(structure::AbstractBlockStructure, stage::Integer)
    return expected(scenarioproblems(structure, stage))
end
function expected(structure::AbstractBlockStructure{2})
    return expected(scenarioproblems(structure, 2))
end
function scenario_type(structure::AbstractBlockStructure, stage::Integer)
    return scenario_type(scenarioproblems(structure, stage))
end
function scenario_type(structure::AbstractBlockStructure{2})
    return scenario_type(scenarioproblems(structure, 2))
end
function probability(structure::AbstractBlockStructure, stage::Integer, scenario_index::Integer)
    return probability(scenarioproblems(structure, stage), scenario_index)
end
function probability(structure::AbstractBlockStructure{2}, scenario_index::Integer)
    return probability(scenarioproblems(structure, 2), scenario_index)
end
function stage_probability(structure::AbstractBlockStructure, stage::Integer)
    return probability(scenarioproblems(structure, stage))
end
function stage_probability(structure::AbstractBlockStructure{2})
    return probability(scenarioproblems(structure, 2))
end
function subproblem(structure::AbstractBlockStructure, stage::Integer, scenario_index::Integer)
    return subproblem(scenarioproblems(structure, stage), scenario_index)
end
function subproblem(structure::AbstractBlockStructure{2}, scenario_index::Integer)
    return subproblem(scenarioproblems(structure, 2), scenario_index)
end
function subproblems(structure::AbstractBlockStructure, stage::Integer)
    return subproblems(scenarioproblems(structure, stage))
end
function subproblems(structure::AbstractBlockStructure{2})
    return subproblems(scenarioproblems(structure, 2))
end
function num_subproblems(structure::AbstractBlockStructure, stage::Integer)
    return num_subproblems(scenarioproblems(structure, stage))
end
function num_subproblems(structure::AbstractBlockStructure{2})
    return num_subproblems(scenarioproblems(structure, 2))
end
function num_scenarios(structure::AbstractBlockStructure, stage::Integer)
    return num_scenarios(scenarioproblems(structure, stage))
end
function num_scenarios(structure::AbstractBlockStructure{2})
    return num_scenarios(scenarioproblems(structure, 2))
end
deferred(structure::AbstractBlockStructure{N}) where N = deferred(structure, Val(N))
deferred(structure::AbstractBlockStructure, ::Val{1}) = deferred_first_stage(structure)
function deferred(structure::AbstractBlockStructure, ::Val{N}) where N
    return deferred_stage(structure, N) || deferred(structure, Val(N-1))
end
deferred_first_stage(structure::AbstractBlockStructure) = false
function deferred_stage(structure::AbstractBlockStructure{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 && return deferred_first_stage(structure)
    num_subproblems(structure, stage) < num_scenarios(structure, stage)
end
function distributed(structure::AbstractBlockStructure{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage == 1 && return false
    return distributed(scenarioproblems(structure, stage))
end

# ========================== #

# Setters
# ========================== #
function update_known_decisions!(structure::AbstractBlockStructure, change::DecisionModification, stage::Integer, scenario_index::Integer)
    update_known_decisions!(scenarioproblems(structure, stage), change, scenario_index)
    return nothing
end
function add_scenario!(structure::AbstractBlockStructure, stage::Integer, scenario::AbstractScenario)
    add_scenario!(scenarioproblems(structure, stage), scenario)
    return nothing
end
function add_worker_scenario!(structure::AbstractBlockStructure, stage::Integer, scenario::AbstractScenario, w::Integer)
    add_scenario!(scenario(structure, stage), scenario, w)
    return nothing
end
function add_scenario!(scenariogenerator::Function, structure::AbstractBlockStructure, stage::Integer)
    add_scenario!(scenariogenerator, scenarioproblems(structure, stage))
    return nothing
end
function add_worker_scenario!(scenariogenerator::Function, structure::AbstractBlockStructure, stage::Integer, w::Integer)
    add_scenario!(scenariogenerator, scenarioproblems(structure, stage), w)
    return nothing
end
function add_scenarios!(structure::AbstractBlockStructure, stage::Integer, scenarios::Vector{<:AbstractScenario})
    add_scenarios!(scenarioproblems(structure, stage), scenarios)
    return nothing
end
function add_worker_scenarios!(structure::AbstractBlockStructure, stage::Integer, scenarios::Vector{<:AbstractScenario}, w::Integer)
    add_scenarios!(scenarioproblems(structure, stage), scenarios, w)
    return nothing
end
function add_scenarios!(scenariogenerator::Function, structure::AbstractBlockStructure, stage::Integer, n::Integer)
    add_scenarios!(scenariogenerator, scenarioproblems(structure, stage), n)
    return nothing
end
function add_worker_scenarios!(scenariogenerator::Function, structure::AbstractBlockStructure, stage::Integer, n::Integer, w::Integer)
    add_scenarios!(scenariogenerator, scenarioproblems(structure, stage), n, w)
    return nothing
end
function sample!(structure::AbstractBlockStructure, stage::Integer, sampler::AbstractSampler, n::Integer)
    sample!(scenarioproblems(structure, stage), sampler, n)
    return nothing
end
# ========================== #

# Includes
# ========================== #
include("scenarioproblems.jl")
include("vertical.jl")
include("horizontal.jl")


text:
#################################################################################################
# demo 3: Thresholding with OpenCV.jl
#
# http://docs.opencv.org/trunk/modules/imgproc/doc/miscellaneous_transformations.html#threshold
# http://docs.opencv.org/trunk/doc/tutorials/imgproc/threshold/threshold.html#basic-threshold
#################################################################################################

# Julia filename with full path
inputfile = joinpath(Pkg.dir("OpenCV"), "./test/images/lena.png")
outfile = joinpath(Pkg.dir("OpenCV"), "./test/images/lena.jpeg")

#Main code
cxx"""
// global variables
int threshold_value = 0;
int const max_value = 255;

void imThreshold(const char *inputfile, const char *outfile) {

  const std::string ifname = inputfile;
  const std::string ofname = outfile;

  cv::Mat src = cv::imread(ifname);
  cv::Mat img_gray, img_thresh;

  if (src.empty()) {
    std::cout << "Can not open the file!" << std::endl;
    exit(0);
  }

  // Create windows to display images
  cv::namedWindow("Original", cv::WINDOW_AUTOSIZE);
  cv::namedWindow("Thresholded", cv::WINDOW_AUTOSIZE);
  cv::namedWindow("Blur", cv::WINDOW_AUTOSIZE);

  // Convert the image to Gray
  cv::cvtColor(src, img_gray, cv::COLOR_RGB2GRAY);

  // Gaussian filtering
  cv::GaussianBlur(img_gray, img_gray, cv::Size(5,5), 0);

  // Display original and blurred images
  cv::imshow("Original", src);
  cv::imshow("Blur", img_gray);

       // cv::THRESH_BINARY
       // cv::THRESH_OTSU
       // cv::THRESH_BINARY_INV
       // cv::THRESH_TRUNC
       // cv::THRESH_TOZERO
       // cv::THRESH_TOZERO_INV

  cv::threshold(img_gray, img_thresh, threshold_value, max_value, cv::THRESH_BINARY + cv::THRESH_OTSU);

  cv::imshow("Thresholded", img_thresh);

  cv::waitKey(0);
  cv::imwrite(ofname, img_thresh);
  cv::destroyAllWindows();
}
"""

@cxx imThreshold(pointer(inputfile), pointer(outfile))


text:
# Adapted from https://github.com/tpapp/TransformVariables.jl/blob/master/src/special_arrays.jl

using LinearAlgebra
const CorrCholeskyUpper = CorrCholeskyFactor

export CorrCholeskyLower

"""
    CorrCholeskyLower(n)
Lower Cholesky factor of a correlation matrix of size `n`.
Transforms ``n(n-1)/2`` real numbers to an ``n×n`` lower-triangular matrix `L`, such that
`L*L'` is a correlation matrix (positive definite, with unit diagonal).

# Notes
If
- `z` is a vector of `n` IID standard normal variates,
- `σ` is an `n`-element vector of standard deviations,
- `L` is obtained from `CorrCholeskyLower(n)`,
then `Diagonal(σ) * L * z` is a zero-centered multivariate normal variate with the standard deviations `σ` and
correlation matrix `L * L'`.
"""
struct CorrCholeskyLower <: TV.VectorTransform
    n::Int
end

TV.dimension(t::CorrCholeskyLower) = TV.dimension(CorrCholeskyUpper(t.n))

# TODO: For now we just transpose the CorrCholeskyUpper result. We should
# consider whether it can help performance to implement this directly for the
# lower triangular case
function TV.transform_with(flag::TV.LogJacFlag, t::CorrCholeskyLower, x::AbstractVector, index)
    U, ℓ, index = TV.transform_with(flag, CorrCholeskyUpper(t.n), x, index)
    return U', ℓ, index
end

TV.inverse_eltype(t::CorrCholeskyLower, L::LowerTriangular) = TV.extended_eltype(L)

function TV.inverse_at!(x::AbstractVector, index, t::CorrCholeskyLower, L::LowerTriangular)
    return TV.inverse_at!(x, index, CorrCholeskyUpper(t.n), L')
end


text:
printstyled(stdout,"  extended methods\n", color=:light_green)

printstyled(stdout,"    getindex\n", color=:light_green)
(S,T) = mktestseis()
# @test findid(T, S) == [0, 0, 4, 5]

printstyled(stdout,"    getindex + Int on SeisData ==> SeisChannel\n", color=:light_green)
C = S[1]
@test S[firstindex(S)] == C
i_targ = 3
C = S[i_targ]
test_fields_preserved(C,S,i_targ)

printstyled(stdout,"    getindex + Range on SeisData ==> SeisData\n", color=:light_green)
D = S[i_targ:i_targ+1]
test_fields_preserved(D,S,2,i_targ+1)

printstyled(stdout,"    setindex!\n", color=:light_green)
A = SeisData(3)
setindex!(A, C, 3)
test_fields_preserved(C, A, 3)
A[1:2]=D
test_fields_preserved(A, S, 2, i_targ+1)
test_fields_preserved(C, S, 3)

printstyled(stdout,"    in\n", color=:light_green)
@test ("XX.TMP01.00.BHZ" in S.id)

printstyled(stdout,"      findid\n", color=:light_green)
@test ≈(findid("CC.LON..BHZ",S), findid(S,"CC.LON..BHZ"))
@test ≈(findid(S,"CC.LON..BHZ"), 4)
@test ≈(findid(S,C), 3)

printstyled(stdout,"    isempty\n", color=:light_green)
D = SeisData()
@test isempty(D)

printstyled(stdout,"    append!\n", color=:light_green)
(S,T) = mktestseis()
append!(S, T)
sizetest(S, 9)

C = deepcopy(S[4])
deleteat!(S, 4)
sizetest(S, 8)
@test ≈(length(findall(S.name.==C.name)),1)

C = deepcopy(S[3:4])
deleteat!(S,3:4)
nt = 6
@test ≈(S.n, nt)
@test ≈(maximum([length(getfield(S,i)) for i in datafields]), nt)
@test ≈(minimum([length(getfield(C,i)) for i in datafields]), 2)
@test length(findall(S.name.==C.name[1])).*length(findall(S.id.==C.id[1])) == 0
@test length(findall(S.name.==C.name[2])).*length(findall(S.id.==C.id[2])) == 1

s = "CC.LON..BHZ"
delete!(S, s)
sizetest(S, 5)

s = r"EH"
# @test_throws BoundsError S - s
delete!(S, s, exact=false)
sizetest(S, 2)

# untested methods in SeisData
for i = 1:5
  S = randSeisData()
  @test sizeof(S) > 0

  r = S.id[1]
  U = S - r
  @test sizeof(U) < sizeof(S)
  @test S.n == U.n + 1

  T = pull(S,r)
  @test isa(T, SeisChannel)
  @test(U == S)

  S = randSeisData(5)
  C = S[1]
  @test firstindex(S) == 1
  @test S[firstindex(S)] == S[1] == C
  @test length(S) == S.n == lastindex(S)

  push!(S, SeisChannel())
  T = prune(S)
  prune!(S)
  @test S == T
  @test S.n == T.n == 5
end

H = randSeisHdr()
@test sizeof(H) > 0
clear_notes!(H)
@test length(H.notes) == 1

(S,T) = mktestseis()
U = S-T
sizetest(S,5)
sizetest(T,4)
sizetest(U,3)

(S,T) = mktestseis()
U = deepcopy(S)
deleteat!(U, 1:3)
@test (S - [1,2,3]) == U
sizetest(S,5)

(S,T) = mktestseis()
@test in("UW.SEP..EHZ",S)
U = S[3:S.n]
V = deepcopy(S)
deleteat!(S, 1:2)
deleteat!(V, [1,2])
@test S == U == V
delete!(S, "CC.", exact=false)
delete!(U,V)

(S,T) = mktestseis()
X = deepcopy(T)
U = pull(S, 5)
@test U.id == "UW.SEP..EHZ"
sizetest(S, 4)
@test findid("UW.SEP..EHZ", S.id) == 0

U = pull(S, 3:4)
sizetest(S, 2)
@test findid(U.id[1], S.id) == 0
@test findid(U.id[2], S.id) == 0

V = pull(T, [2,3,4])
sizetest(T, 1)
@test findid(V.id[1], T.id) == 0
@test findid(V.id[2], T.id) == 0
@test findid(V.id[3], T.id) == 0
V = sort(V)

deleteat!(X, 1)
# @test findid(V,X) == [2,3,1]
Y = sort(X)
@test V == Y

# added 2019-02-23
S = SeisData(randSeisData(5), SeisChannel(), SeisChannel(),
    SeisChannel(id="UW.SEP..EHZ", name="Darth Exploded",
    loc=GeoLoc(lat=46.1967, lon=-122.1875, el=1440.0), t=[1 0; 1024 0], x=rand(1024)))
prune!(S)
@test (S.n == 6)
J = findchan("EHZ",S)
@test (6 in J)

printstyled("  SeisChannel methods\n", color=:light_green)

id = "UW.SEP..EHZ"
name = "Darth Exploded"

Ch = randSeisChannel()
Ch.id = id
Ch.name = name
S = SeisData(Ch)

@test in(id, Ch) == true
@test isempty(Ch) == false
@test convert(SeisData, Ch) == SeisData(Ch)
@test findid(Ch, S) == 1
@test sizeof(Ch) > 0
@test lastindex(S) == 1

Ch.gain = 1.0
Ch.fs = 0.0
@test isempty(Ch) == false


text:
using Cambrian
import YAML

cfg = get_config("test.yaml")

@testset "Random fitness" begin
    e = Evolution{FloatIndividual}(cfg)
    @show e
    evaluate!(e::Evolution) = random_evaluate!(e)
    populate!(e::Evolution) = oneplus_populate!(e)
    @test length(e.population) == prod(cfg.population.size)
    for i in e.population
        @test all(i.fitness .== -Inf)
    end

    evaluate!(e)
    fits = [i.fitness[1] for i in e.population]
    @show fits
    evaluate!(e)
    # random fitness, all values should change
    for i in eachindex(e.population)
        @show i
        @test fits[i] != e.population[i].fitness[1]
    end

    fits = copy([i.fitness[:] for i in e.population])
    step!(e)
    # random fitness, all values should change
    for i in eachindex(e.population)
        @test fits[i] != e.population[i].fitness[1]
    end
end


text:
"""
Provides calculations for orbit maneuvers.
"""
module Maneuvers

using Reexport
@reexport using ..AstrodynamicsCore

include("../Misc/DocStringExtensions.jl")
include("../Misc/UnitfulAliases.jl")

using ..NBody
using ..TwoBody
using LinearAlgebra: norm, cross, ×, dot, ⋅

export AbstractManeuver, TwoBodyManeuver, ConstantManeuver
export escape_radius, escape_velocity, escape_time, escape_path_length

include("maneuver_types.jl")
include("twobody_maneuver_calculations.jl")

end


text:
using Distributed
addprocs(16)
@everywhere include("get_srb.jl")
get_dist()

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
using SparseArrays

# iterate through units in "reading order":
#   if enemy adjacent, attack, else move
# move:
#   find (reachable) squares in range of enemy
#   find closest by path length (break ties by reading order)
#   take first step
#
# _COULD_ optimize this a bit: do like a BFS/Dijkstra's from the start,
# terminate when you find an enemy...but then how to handle the actual
# pathfinding part?  I guess if you build the minimum distance _tree_ at the
# same time...

# keep track of
# * map (empty squares)
# * list of units
# * array of units (inverse index location -> unit)

mutable struct Board
    openmap::BitArray{2}
    units::SparseMatrixCSC{Int,Int}
    turns::Int
    attacks::Tuple{Int,Int}
end

function Base.show(io::IO, board::Board)
    xs, ys = axes(board.openmap)
    for y in ys
        hps = []
        for x in xs
            unit = board.units[x,y]
            if unit == 0
                print(io, board.openmap[x,y] ? '.' : '#')
            else
                ch = unit > 0 ? 'E' : 'G'
                print(io, ch)
                push!(hps, "$ch ($(abs(unit)))")
            end
        end
        print(io, "   ")
        println(io, join(hps, " "))
    end
end

test_input = """
#######
#.G...#
#...EG#
#.#.#G#
#..G#E#
#.....#
#######
"""

parseinput(input) = reduce(hcat, collect.(split(strip(input, '\n'))))

Board(input::String, args...) = Board(parseinput(input), args...)
function Board(input::Matrix{Char}, attacks=(3,-3))
    openmap = input .!= '#'

    elves = sparse(200 .* (input .== 'E'))
    goblins = sparse(200 .* (input .== 'G'))
    units = elves .- goblins
    Board(openmap, units, 0, attacks)
end

neighbors(x) = (x + d for d in CartesianIndex.(((0,-1), (-1,0), (1,0), (0,1))))

# keep track of a frontier (nodes that are waiting to be visited).  in a queue
# but ... don't need to sort because the distances are the same... so you work
# through all the distance = 1 first, and each of those pushes a d=2, then you
# hit those and they each push a d=3.
#
# but...what do we need to keep track of in order to choose the MOVE??  keep
# track of the ancestor for each visited location

function next_step(start, board::Board)

    openmap = board.openmap
    units = board.units

    goodness = sign(units[start])

    queue = CartesianIndex{2}[start]
    ancestors = zeros(CartesianIndex{2}, size(openmap))

    while !isempty(queue)
        current = popfirst!(queue)
        # println("current: $(current.I)")
        # current tile is not done.  check each neighbor for valid
        for neighbor in neighbors(current)
            if sign(units[neighbor]) * goodness == -1
                # we're on a neighboring tile to opposite badness, so done.
                # println("enemy found at $(neighbor.I) (from $(current.I))")
                # walk back to position where current is ancestor:
                # print("  walking back: $(current.I)")
                while ancestors[current] != start && !iszero(ancestors[current])
                    current = ancestors[current]
                    # print(" -> $(current.I)")
                end
                # println()
                return current
            elseif openmap[neighbor] && units[neighbor] == 0 && iszero(ancestors[neighbor])
                # println("  neighbor $(neighbor.I): added to queue")
                # neighbor is valid and unvisited. record ancestor and add to queue
                ancestors[neighbor] = current
                push!(queue, neighbor)
            else
                # println("invalid!")
            end
        end
    end
    # no reachable enemy found so don't move, return start location as "next step"
    return start
end

function turn!(board, unit)
    units = board.units
    next = next_step(unit, board)
    if iszero(board.units[next])
        # println("moving: $(unit.I) -> $(next.I)")
        # move
        units[unit], units[next], unit = units[next], units[unit], next
    end

    enemies = [(abs(units[n]), n)
               for n
               in neighbors(unit)
               if sign(units[n]*units[unit])==-1]
    if !isempty(enemies)
        hp, location = minimum(enemies)
        attack = units[unit] > 0 ? board.attacks[1] : board.attacks[2]
        if hp <= abs(attack)
            units[location] = 0 # dead
        else
            units[location] += attack
        end
    end
end

done(board::Board) = all(u>=0 for u in board.units) || all(u<=0 for u in board.units)

function turn!(board::Board)
    for unit in findall(!iszero, board.units)
        done(board) && return
        # println("unit at $(unit.I)")
        turn!(board, unit)
        # println(board)
    end
    board.turns += 1
    return board
end

board = Board(test_input)

for _ in 1:23
    turn!(board)
end

for _ in 24:47
    turn!(board)
end

test_input2 = """
#######
#G..#E#
#E#E.E#
#G.##.#
#...#E#
#...E.#
#######
"""

score(board::Board) = board.turns * abs(sum(board.units))

function star1(input)
    board = Board(input)
    while turn!(board) !== nothing
    end
    return score(board), board
end

star1(test_input)
star1(test_input2)

"""
#########
#G......#
#.E.#...#
#..##..G#
#...##..#
#...#...#
#.G...G.#
#.....G.#
#########
""" |> star1

"""
#######
#.E...#
#.#..G#
#.###.#
#E#G#G#
#...#G#
#######
""" |> star1

read("15.input", String) |> star1

################################################################################
# star 2: boost attack power of elves until no elves die
n_elves(board::Board) = sum(u > 0 for u in board.units)

function star2(input)
    for attack = 4:100
        board = Board(input, (attack, -3))
        n0 = n_elves(board)
        while n_elves(board)==n0 && turn!(board) !== nothing end
        println("attack $attack, score $score, board\n$board)")
        if n_elves(board) == n0
            return board.turns * abs(sum(board.units)), board
        end
    end
end

"""
#######
#.G...#
#...EG#
#.#.#G#
#..G#E#
#.....#
#######
""" |> star2

"""
#######
#E..EG#
#.#G.E#
#E.##E#
#G..#.#
#..E#.#
#######
""" |> star2

"""
#######
#E.G#.#
#.#G..#
#G.#.G#
#G..#.#
#...E.#
#######
""" |> star2

"""
#########
#G......#
#.E.#...#
#..##..G#
#...##..#
#...#...#
#.G...G.#
#.....G.#
#########
""" |> star2

read("15.input", String) |> star2


text:
module ComplementaryDistributions

using Distributions
using Random
using SpecialFunctions
using StatsFuns: twoπ

export GeneralizedInverseGaussian, PolyaGamma, LaplaceTransformDistribution
include("generalizedinversegaussian.jl")
include("polyagamma.jl")
include("lap_transf_dist.jl")

end


text:
function find_datasets(root::String, query::String)
    dataset_folder = root
    for entry in readdir(dataset_folder)
        path = joinpath(dataset_folder, entry)
        isdir(path) || continue
        if entry == query
            return path
        end
    end
    return nothing
end

"""
Get the appropriate dataset from anywhere we can find.
Available options: ModelNet10, ModelNet40
"""
function dataset(which::String, root::String)
    if which === "ModelNet10PCloud"
        path = find_datasets(root, "modelnet40_normal_resampled")
        if path == nothing
            download("ModelNet", root)
            path = find_datasets(root, "modelnet40_normal_resampled")
        end
        return path
    elseif which == "ModelNet40PCloud"
        path = find_datasets(root, "modelnet40_normal_resampled")
        if path == nothing
            download("ModelNet", root)
            path = find_datasets(root, "modelnet40_normal_resampled")
        end
        return path
    else
        error("Autodetection not supported for $(which)")
    end
end

function download(which::String, root::String)
    if which == "ModelNet"
        local_path = joinpath(root, "modelnet40_normal_resampled.zip")
        dir_path = joinpath(root) #TODO check if dir_path exists and if not then create one.
        if(!isdir(joinpath(dir_path, "modelnet40_normal_resampled")))
            if(!isfile(local_path))
                # dataset prepared by authors of pointnet2
                Base.download("https://shapenet.cs.stanford.edu/media/modelnet40_normal_resampled.zip", local_path)
            end
            run(`unzip $local_path -d $dir_path`)
        end
    else
        error("Download not supported for $(which)")
    end
end