text:
isempty(line::AbstractString) = !ismatch(r"[^\s\t\n\r]+", line)

iscomment(line::AbstractString) = ismatch(r"^#", line)

istitle(line::AbstractString) = ismatch(r"^\*\*\*(.*)MODEL (STATES|PARAMETERS|ODES)", line)

is_ode_name(name::AbstractString) = ismatch(r"^d/dt\((.*)\)$", name)

function get_title_type(line::AbstractString)
  if ismatch(r"^\*\*\*(.*)MODEL STATES(.*)", line)
    return "state"
  elseif ismatch(r"^\*\*\*(.*)MODEL PARAMETERS(.*)", line)
    return "parameter"
  elseif ismatch(r"^\*\*\*(.*)MODEL ODES(.*)", line)
    return "ode"
  else
    throw(ParseError("Unknown title \""*line*"\" type in model specification"))
  end
end

function parse_line(line::AbstractString, linetype::AbstractString, model::OdeModel)
  if linetype == "state"
    nTokens = parse_state_line(line, model.states)
  elseif linetype == "parameter"
    nTokens = parse_parameter_line(line, model.parameters)
  elseif linetype == "ode"
    nTokens = parse_ode_line(line, model.odes)
  else
    throw(ParseError("Unknown title \""*line*"\" type in model specification"))
  end

  return nTokens
end

function parse_state_line(line::AbstractString, dict::Dict{AbstractString, Float64})
  tokens = split(line, '=')
  nTokens = length(tokens)

  if nTokens != 2
    throw(ParseError("Line \""*line*"\" has wrong format"))
  end

  name, value = [strip(i) for i in tokens]

  haskey(dict, name) ? throw(KeyError("\""*name*"\" specified more than once")) : (dict[name] = Float64(parse(value)))

  return nTokens
end

function parse_parameter_line(line::AbstractString, dict::Dict{AbstractString, Float64})
  tokens = split(line, '=')
  nTokens = length(tokens)

  if nTokens == 1
    dict[strip(tokens)] = NaN
  elseif nTokens == 2
    name, value = [strip(i) for i in tokens]

    haskey(dict, name) ? throw(KeyError("\""*name*"\" specified more than once")) : (dict[name] = Float64(parse(value)))
  else
    throw(ParseError("Line \""*line*"\" has wrong format"))
  end

  return nTokens
end

function parse_ode_line(line::AbstractString, dict::Dict{AbstractString, NSE})
  tokens = split(line, '=')
  nTokens = length(tokens)

  if nTokens != 2
    throw(ParseError("Line \""*line*"\" has wrong format"))
  end

  name, value = [strip(i) for i in tokens]

  is_ode_name(name) ? (name = name[6:end-1]) : throw(KeyError("Time derivative \""*name*"\" has wrong format"))

  haskey(dict, name) ? throw(KeyError("\""*name*"\" specified more than once")) : (dict[name] = parse(value))

  return nTokens
end

function parse_model(file::AbstractString)
  odeModel = OdeModel(file, Dict{AbstractString, Float64}(), Dict{AbstractString, Float64}(), Dict{AbstractString, NSE}())

  s = open(file)

  while !eof(s)
    line = readline(s)
    if isempty(line) || iscomment(line)
      continue
    elseif istitle(line)
      title = get_title_type(line)
    else
      parse_line(line, title, odeModel)
    end
  end

  close(s)

  if sort(collect(keys(odeModel.states))) != sort(collect(keys(odeModel.odes)))
    throw(KeyError("States and their time derivatives in the ODE model are inconsistent"))
  end

  return odeModel
end


text:
module UnifiedPseudopotentialFormat

include("name.jl")
include("file.jl")
include("PSlibrary.jl")

end

ext: jl


text:
println("Enter the array size ")
n = parse(Int32,readline())
a=[]
println("Enter the array values ")
for i = 1:n
          temp =  parse(Int32,readline())
          append!(a, temp)
       end
println("Enter the element to search ")
s =  parse(Int32,readline())
flag = 0
for i = 1:n
          if s == a[i]
              println(i)
              global flag += 1
              break
            end
       end
if flag == 0
    println("Not found")
    end

max_stars_repo_name: PRIYANSHU-AMBITION/Linear-Search

max_issues_repo_name: PRIYANSHU-AMBITION/Linear-Search

max_forks_repo_name: PRIYANSHU-AMBITION/Linear-Search


text:
using Test

function get_module_2_excersise_1_author_name()
    return "PUT YOUR NAME HERE"
end

function calculate_sin_series_term(i, x)
    return 0 # implement your sin series term here
end

function calculate_series_sin(x)
    series_expansion_value = 0
    if x > 2*pi
        throw(error("x must be less than 2π"))
    end
    for i in 0:100
        # implement the loop here
        i+= 1
    end
    return series_expansion_value
end


text:
using ITensors, Test

@testset "Fermions" begin
  @testset "parity_sign function" begin

    # Full permutations
    p1 = [1, 2, 3]
    @test ITensors.parity_sign(p1) == +1
    p2 = [2, 1, 3]
    @test ITensors.parity_sign(p2) == -1
    p3 = [2, 3, 1]
    @test ITensors.parity_sign(p3) == +1
    p4 = [3, 2, 1]
    @test ITensors.parity_sign(p4) == -1

    ## Partial permutations
    p5 = [2, 7]
    @test ITensors.parity_sign(p5) == +1
    p6 = [5, 3]
    @test ITensors.parity_sign(p6) == -1
    p7 = [1, 9, 3, 10]
    @test ITensors.parity_sign(p7) == -1
    p8 = [1, 12, 9, 3, 11]
    @test ITensors.parity_sign(p8) == +1
  end
end


text:
package java.lang;

public class NullPointerException extends RuntimeException {
    public NullPointerException() { }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
"""
    test_chan_ext(file::String, fmt::String, id::String, fs::Float64, i::Int64, t0::Int64)

Test for Issue 34
"""

function test_chan_ext(file::String, fmt::String, id::String, fs::Float64, i::Int64, t0::Int64)

  # Comparison struct
  S = verified_read_data(fmt, file)
  j = findid(id, S)
  if j == 0
    println("Test failing; cannot find ID!")
    println(rpad("File: ", 12), file)
    println(rpad("Format: ", 12), fmt)
    println("Channel ID: ", id)
    println("Output IDs: ", S.id)
    throw("Test failed; ID not found!")
  end

  # Create SeisData struct with data present using correct headers
  ny = max(100, round(Int64, 10*fs))
  y = rand(Float32, ny)
  S1 = SeisData(
    SeisChannel(
      id = id,
      fs = fs,
      gain = S.gain[j],
      units = S.units[j],
      loc = deepcopy(S.loc[j]),
      resp = deepcopy(S.resp[j]),
      t = [1 rand(0:10000000); ny 0],
      x = copy(y)
      )
    )

  # Read into S1
  verified_read_data!(S1, fmt, file)
  t1 = t_expand(S1.t[i], S.fs[i])[ny+1]
  t2 = S.t[j][1,2]

  # Check that read preserves time, data
  @test S1.x[i][1:ny] == y
  @test S1.x[i][ny+1:end] == S.x[j]
  @test t0 == t1 == t2
  return nothing
end


text:
immutable SimulationPeriod <: Period
  value :: Float64
  function SimulationPeriod(value::Number=0)
    new(value)
  end
end

immutable SimulationInstant <: Dates.Instant
  periods :: SimulationPeriod
end

immutable SimulationTime <: TimeType
  instant :: SimulationInstant
end

function SimulationTime(value::Number=0)
  SimulationTime(SimulationInstant(SimulationPeriod(value)))
end

(==)(x::SimulationTime, y::SimulationTime) = x.instant.periods.value == y.instant.periods.value

function isless(t1::SimulationTime, t2::SimulationTime) :: Bool
  t1.instant.periods.value < t2.instant.periods.value
end

(+)(t::SimulationTime, p::SimulationPeriod)=SimulationTime(t.instant.periods.value + p.value)

function show(io::IO, t::SimulationTime)
  print(io, "$(t.instant.periods.value)")
end


text:
function sinematrix(s::Vector{Float64}, c::Vector{Float64})
	return s.*c' .- c.*s'
end

function cosmatrix(s::Vector{Float64}, c::Vector{Float64})
	return c.*c' + s.*s'
end

function cmatrix(κ::Float64, A::Matrix{Float64}, S::Matrix{Float64})
	return κ.*A.*S
end

function omegadot(ω::Vector{Float64}, θ::Vector{Float64}, A::Matrix{Float64}, P::Vector{Float64}, I::Float64, D::Float64, κ::Float64)
	s = sin.(θ)
	c = cos.(θ)
	smat = sinematrix(s, c)
	C = cmatrix(κ, A, smat)
	return (P .- D.*ω .- C*ones(length(ω)))./I
end

function fswing(ψ::Vector{Float64}, A::Matrix{Float64}, P::Vector{Float64}, I::Float64, D::Float64, κ::Float64)
	m = length(ψ)
	n = div(m,2)
	ω = ψ[1:n]
	θ = ψ[n+1:end]
	ω_dot = omegadot(ω, θ, A, P, I, D, κ)
	return vcat(ω_dot, ω)
end

function fsteadystate(θ::Vector{Float64}, A::Matrix{Float64}, P::Vector{Float64}, κ::Float64)
	n = length(θ)
	s = sin.(θ)
	c = cos.(θ)
	smat = sinematrix(s, c)
	C = cmatrix(κ, A, smat)
	return P .- C*ones(n)
end

function mainjac(θ::Vector{Float64}, E::Matrix{Float64}, I::Float64, κ::Float64)
	weights = κ*cos.(E'*θ)
	W = Diagonal(weights)
	return -E*W*E'
end

function dfswing(ψ::Vector{Float64}, A::Matrix{Float64}, I::Float64, D::Float64, κ::Float64)
	m = length(ψ)
	n = div(m,2)
	θ = ψ[n+1:end]
	DF = zeros(m,m)
	DF[1:n,1:n] = D.*Matrix{Float64}(I, n, n)./I
	DF[1:n,n+1:end] = mainjac(θ, A, I, κ)./I
	DF[n+1:end,1:n] = Matrix{Float64}(I, n, n)
	return DF
end

function edgepower(θ::Vector{Float64}, E::Matrix{Float64}, κ::Float64)
	Δθ = transpose(E)*θ
	return κ.*sin.(Δθ)
end

function centralityvstransient(n::Int64, ns::Int64, nd::Int64, k::Int64, q::Float64, κ::Float64, I::Float64, D::Float64)

	#Set up network
	net, E = smallworldnetwork(n, k, q)
	ψ = rand(2*n)
	A = adjacencymatrix(E)
	sources, sinks = sourcesinklocs(ns, nd, n)
	P = sourcesinkvector(sources, sinks, n, 1.0)
	obj1 = NetStore(A, P, κ)

	#Find initial steady state
	f(u, p, t) = fswing(u, obj1.adjmat, P, I, D, κ)
	tspan = (0.0, 500.0)
	prob = ODEProblem(f, ψ, tspan)
	sol = solve(prob,reltol=1e-8,abstol=1e-8,save_everystep=false,dense=false,save_end=true)
	ψ0 = sol[:,end]
	θ0 = sol[n+1:end,end]
	flow = edgepower(θ0, E, κ)
	meanflow = mean(abs.(flow))

	#Set up callback function
	function conditiontemp(u1, t1, integrator, ob::NetStore)
		n = div(length(u1),2)
		θ = u1[n+1:end,end]
		resid = fsteadystate(θ, ob.adjmat, ob.pvec, ob.kappa)
		if norm(resid,2) < 0.001
			return true
		else
			return false
		end
	end
	condition(u,t,integrator) = conditiontemp(u,t,integrator,obj1)
	affect!(integrator) = terminate!(integrator)

	m = numedges(net)
	relflow = zeros(m)
	transienttimes = zeros(m)
	distances = zeros(m)
	for d=1:m

		#Delete edge and find its load
		E1 = E[1:end, 1:end .≠ d]
		A1 = adjacencymatrix(E1)
		obj1.adjmat = A1
		relflow[d] = abs(flow[d])/meanflow
		flow_old = flow[1:end .≠ d]

		#Time step on reduced network
		u0 = ψ0
		tspan2 = (0.0, 5000.0)
		prob2 = ODEProblem(f, u0, tspan2)
		cb = DiscreteCallback(condition, affect!)
		sol2 = solve(prob2, reltol=1e-8,abstol=1e-8,save_everystep=false,dense=false,save_end=true,callback=cb)

		#Find new flow and measure its distance
		theta_new = sol2[n+1:end,end]
		flow_new = edgepower(theta_new, E1, κ)
		sep = abs.(flow_new .- flow_old)
		transienttimes[d] = sol2.t[end]
		distances[d] = norm(sep, 2)
	end
	return relflow, transienttimes, distances

end

function swingfracture!(net::Graph, ψ::Vector{Float64}, nodeset::Vector{Int64}, P::Vector{Float64}, synctol::Float64, alpha::Float64, kappa::Float64, maxflow::Float64)

	tol = 1e-5

	#Check the status of P, see if has any sources or sinks. If not, halt recursion
	nodecardinality = length(nodeset)
	P1 = zeros(Float64, nodecardinality)
	sourcecounter = 0
	sinkcounter = 0
	for i=1:nodecardinality
		vi = nodeset[i]
		P1[i] = P[vi]
		if P1[i] > tol
			sourcecounter += 1
		end
		if P1[i] < -tol
			sinkcounter += 1
		end
	end
	if (sourcecounter == 0) || (sinkcounter == 0)
		return 0
	end

	#Balance P1 homogeneously
	delta = sum(P1)/2.0
	individualdelta_sink = delta/Float64(sinkcounter)
	individualdelta_source = delta/Float64(sourcecounter)
	for i=1:nodecardinality
		if P1[i] < -tol
			P1[i] -= individualdelta_sink
		end
		if P1[i] > tol
			P1[i] -= individualdelta_source
		end
	end

	# #Balance P1 - if there is power excess, increase loads
	# delta = sum(P1)
	# individualdelta = delta/Float64(sinkcounter)
	# for i=1:nodecardinality
	# 	if P1[i] < -tol
	# 		P1[i] -= individualdelta
	# 	end
	# end

	#Find which edges are connected to the node-set
	edgeset = Array{Int64, 1}([])
	for vi ∈ nodeset
		v = net.nodelist[vi]
		for e ∈ v.outedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
		for e ∈ v.inedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
	end
	edgecardinality = length(edgeset)

	#Construct E1
	E1 = zeros(Float64, nodecardinality, edgecardinality)
	for i=1:nodecardinality
		vi = nodeset[i]
		v = net.nodelist[vi]
		for j=1:edgecardinality
			if edgeset[j] ∈ v.outedges
				E1[i,j] = 1.0
			end
			if edgeset[j] ∈ v.inedges
				E1[i,j] = -1.0
			end
		end
	end

	#Set up the storage object
	A1 = adjacencymatrix(E1)
	obj1 = NetStore2(A1, P1, E1, kappa, synctol, alpha, maxflow, true, false)

	#Set up callback function
	function conditiontemp(u1, t1, integrator, ob::NetStore2)
		n = div(length(u1),2)
		ω = u1[1:n]
		θ = u1[n+1:end]

		if (t1 > 0.0)
			#Has network synchronised?
			if norm(ω,2) > ob.synctol
				ob.sync = false
				return true
			end

			#Have any edges tripped the capacity
			flow = edgepower(θ, ob.incidence, ob.kappa)./ob.maxflow
			for i=1:length(flow)
				if abs(flow[i]) > ob.alpha
					ob.trip = true
					return true
				end
			end

			#Has the solution converged
			resid = fsteadystate(θ, ob.adjmat, ob.pvec, kappa)
			if norm(resid,2) < 1e-6
				return true
			end
		end

		#If the network is still syncing and no edge trips, continue
		return false
	end
	condition(u,t,integrator) = conditiontemp(u,t,integrator,obj1)
	affect!(integrator) = terminate!(integrator)

	#Time step
	f(u, p, t) = fswing(u, obj1.adjmat, P1, 1.0, 1.0, kappa)
	tspan = (0.0, 500.0)
	prob = ODEProblem(f, ψ, tspan)
	cb = DiscreteCallback(condition, affect!)
	sol = solve(prob,reltol=1e-8,abstol=1e-8,save_everystep=false,dense=false,save_end=true,callback=cb)

	#Get the state vector, split into theta and omega
	ψ_1 = sol[:,end]
	ω_1 = ψ_1[1:nodecardinality]
	θ_1 = ψ_1[nodecardinality+1:end]

	#If the solution has not converged, then return
	if obj1.sync == false
		return 0
	end

	#If the solution has converged with no further overloads, then return
	if (obj1.trip == false) && (obj1.sync == true)
		return edgecardinality
	end

	#If there are overloads, but omega still in acceptable range,
	#then compute number of new connected components and continue
	#recursion on each of them
	#First construct new edgeset and incidence matrix
	overloadcounter = 0
	survivors = Array{Int64, 1}([])
	flow = edgepower(θ_1, E1, kappa)./maxflow
	for i=1:edgecardinality
		ei = edgeset[i]
		if abs(flow[i]) > alpha #*abs(Finit[ei])
			net.edgelist[ei].active = false
			overloadcounter += 1
		else
			push!(survivors, i)
		end
	end
	newedgecard = length(survivors)
	E2 = zeros(Float64, nodecardinality, newedgecard)
	for j=1:newedgecard
		E2[:,j] = E1[:,survivors[j]]
	end

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

	#Clear any superfluous data
	E1 = nothing
	E2 = nothing
	Adj = nothing
	obj1 = nothing

	#Call fracture function on each of the components, passing down info
	descendent_edges = 0
	for i=1:components
		chunk = table[i]
		n_subset = length(chunk)
		nodesubset = zeros(Int64, length(chunk))
		ω_subset = zeros(Float64, length(chunk))
		θ_subset = zeros(Float64, length(chunk))
		ψ_subset = zeros(Float64, 2*n_subset)
		for j=1:length(chunk)
			nodesubset[j] = nodeset[chunk[j]]
			ω_subset[j] = ω_1[chunk[j]]
			θ_subset[j] = θ_1[chunk[j]]
		end
		ψ_subset[1:n_subset] = ω_subset
		ψ_subset[n_subset+1:end] = θ_subset
		descendent_edges += swingfracture!(net, ψ_subset, nodesubset, P, synctol, alpha, kappa, maxflow)
	end

	return descendent_edges
end

function swingfracturewithtime!(net::Graph, ψ::Vector{Float64}, nodeset::Vector{Int64}, P::Vector{Float64}, synctol::Float64, alpha::Float64, kappa::Float64, maxflow::Float64, D::Float64)

	tol = 1e-5

	#Check the status of P, see if has any sources or sinks. If not, halt recursion
	nodecardinality = length(nodeset)
	P1 = zeros(Float64, nodecardinality)
	sourcecounter = 0
	sinkcounter = 0
	for i=1:nodecardinality
		vi = nodeset[i]
		P1[i] = P[vi]
		if P1[i] > tol
			sourcecounter += 1
		end
		if P1[i] < -tol
			sinkcounter += 1
		end
	end
	if (sourcecounter == 0) || (sinkcounter == 0)
		return 0, 0.0
	end

	#Balance P1 homogeneously
	delta = sum(P1)/2.0
	individualdelta_sink = delta/Float64(sinkcounter)
	individualdelta_source = delta/Float64(sourcecounter)
	for i=1:nodecardinality
		if P1[i] < -tol
			P1[i] -= individualdelta_sink
		end
		if P1[i] > tol
			P1[i] -= individualdelta_source
		end
	end

	# #Balance P1 - if there is power excess, increase loads
	# delta = sum(P1)
	# individualdelta = delta/Float64(sinkcounter)
	# for i=1:nodecardinality
	# 	if P1[i] < -tol
	# 		P1[i] -= individualdelta
	# 	end
	# end

	#Find which edges are connected to the node-set
	edgeset = Array{Int64, 1}([])
	for vi ∈ nodeset
		v = net.nodelist[vi]
		for e ∈ v.outedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
		for e ∈ v.inedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
	end
	edgecardinality = length(edgeset)

	#Construct E1
	E1 = zeros(Float64, nodecardinality, edgecardinality)
	for i=1:nodecardinality
		vi = nodeset[i]
		v = net.nodelist[vi]
		for j=1:edgecardinality
			if edgeset[j] ∈ v.outedges
				E1[i,j] = 1.0
			end
			if edgeset[j] ∈ v.inedges
				E1[i,j] = -1.0
			end
		end
	end

	#Set up the storage object
	A1 = adjacencymatrix(E1)
	obj1 = NetStore2(A1, P1, E1, kappa, synctol, alpha, maxflow, true, false)

	#Set up callback function
	function conditiontemp(u1, t1, integrator, ob::NetStore2)
		n = div(length(u1),2)
		ω = u1[1:n]
		θ = u1[n+1:end]

		if (t1 > 0.0)
			#Has network synchronised?
			if norm(ω,2) > ob.synctol
				ob.sync = false
				return true
			end

			#Have any edges tripped the capacity
			flow = edgepower(θ, ob.incidence, ob.kappa)./ob.maxflow
			for i=1:length(flow)
				if abs(flow[i]) > ob.alpha
					ob.trip = true
					return true
				end
			end

			#Has the solution converged
			resid = fsteadystate(θ, ob.adjmat, ob.pvec, kappa)
			if norm(resid,2) < 1e-6
				return true
			end
		end

		#If the network is still syncing and no edge trips, continue
		return false
	end
	condition(u,t,integrator) = conditiontemp(u,t,integrator,obj1)
	affect!(integrator) = terminate!(integrator)

	#Time step
	f(u, p, t) = fswing(u, obj1.adjmat, P1, 1.0, D, kappa)
	tspan = (0.0, 500.0)
	prob = ODEProblem(f, ψ, tspan)
	cb = DiscreteCallback(condition, affect!)
	sol = solve(prob,reltol=1e-8,abstol=1e-8,save_everystep=false,dense=false,save_end=true,callback=cb)

	#Get the state vector, split into theta and omega
	ψ_1 = sol[:,end]
	ω_1 = ψ_1[1:nodecardinality]
	θ_1 = ψ_1[nodecardinality+1:end]

	#If the solution has not converged, then return
	if obj1.sync == false
		return 0, 0.0
	end

	#If the solution has converged with no further overloads, then return
	if (obj1.trip == false) && (obj1.sync == true)
		return edgecardinality, sol.t[end]
	end

	#If there are overloads, but omega still in acceptable range,
	#then compute number of new connected components and continue
	#recursion on each of them
	#First construct new edgeset and incidence matrix
	overloadcounter = 0
	survivors = Array{Int64, 1}([])
	flow = edgepower(θ_1, E1, kappa)./maxflow
	for i=1:edgecardinality
		ei = edgeset[i]
		if abs(flow[i]) > alpha #*abs(Finit[ei])
			net.edgelist[ei].active = false
			overloadcounter += 1
		else
			push!(survivors, i)
		end
	end
	newedgecard = length(survivors)
	E2 = zeros(Float64, nodecardinality, newedgecard)
	for j=1:newedgecard
		E2[:,j] = E1[:,survivors[j]]
	end

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

	#Clear any superfluous data
	E1 = nothing
	E2 = nothing
	Adj = nothing
	obj1 = nothing

	#Call fracture function on each of the components, passing down info
	descendent_edges = 0
	times = []
	for i=1:components
		chunk = table[i]
		n_subset = length(chunk)
		nodesubset = zeros(Int64, length(chunk))
		ω_subset = zeros(Float64, length(chunk))
		θ_subset = zeros(Float64, length(chunk))
		ψ_subset = zeros(Float64, 2*n_subset)
		for j=1:length(chunk)
			nodesubset[j] = nodeset[chunk[j]]
			ω_subset[j] = ω_1[chunk[j]]
			θ_subset[j] = θ_1[chunk[j]]
		end
		ψ_subset[1:n_subset] = ω_subset
		ψ_subset[n_subset+1:end] = θ_subset
		xz, tz = swingfracturewithtime!(net, ψ_subset, nodesubset, P, synctol, alpha, kappa, maxflow, D)
		descendent_edges += xz
		push!(times, tz)
	end
	tprime = sol.t[end] + maximum(times)
	return descendent_edges, tprime
end

function swingfracturewithbreakdown!(net::Graph, ψ::Vector{Float64}, nodeset::Vector{Int64}, P::Vector{Float64}, synctol::Float64, alpha::Float64, kappa::Float64, maxflow::Float64)

	tol = 1e-5

	#Check the status of P, see if has any sources or sinks. If not, halt recursion
	nodecardinality = length(nodeset)
	P1 = zeros(Float64, nodecardinality)
	sourcecounter = 0
	sinkcounter = 0
	for i=1:nodecardinality
		vi = nodeset[i]
		P1[i] = P[vi]
		if P1[i] > tol
			sourcecounter += 1
		end
		if P1[i] < -tol
			sinkcounter += 1
		end
	end

	#Find which edges are connected to the node-set
	edgeset = Array{Int64, 1}([])
	for vi ∈ nodeset
		v = net.nodelist[vi]
		for e ∈ v.outedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
		for e ∈ v.inedges
			if (e ∉ edgeset) && (net.edgelist[e].active == true)
				push!(edgeset, e)
			end
		end
	end
	edgecardinality = length(edgeset)
	if (sourcecounter == 0) || (sinkcounter == 0)
		return 0, 0, 0, edgecardinality
	end

	#Balance P1 homogeneously
	delta = sum(P1)/2.0
	individualdelta_sink = delta/Float64(sinkcounter)
	individualdelta_source = delta/Float64(sourcecounter)
	for i=1:nodecardinality
		if P1[i] < -tol
			P1[i] -= individualdelta_sink
		end
		if P1[i] > tol
			P1[i] -= individualdelta_source
		end
	end

	#Construct E1
	E1 = zeros(Float64, nodecardinality, edgecardinality)
	for i=1:nodecardinality
		vi = nodeset[i]
		v = net.nodelist[vi]
		for j=1:edgecardinality
			if edgeset[j] ∈ v.outedges
				E1[i,j] = 1.0
			end
			if edgeset[j] ∈ v.inedges
				E1[i,j] = -1.0
			end
		end
	end

	#Set up the storage object
	A1 = adjacencymatrix(E1)
	obj1 = NetStore2(A1, P1, E1, kappa, synctol, alpha, maxflow, true, false)

	#Set up callback function
	function conditiontemp(u1, t1, integrator, ob::NetStore2)
		n = div(length(u1),2)
		ω = u1[1:n]
		θ = u1[n+1:end]

		if (t1 > 0.0)
			#Has network synchronised?
			if norm(ω,2) > ob.synctol
				ob.sync = false
				return true
			end

			#Have any edges tripped the capacity
			flow = edgepower(θ, ob.incidence, ob.kappa)./ob.maxflow
			for i=1:length(flow)
				if abs(flow[i]) > ob.alpha
					ob.trip = true
					return true
				end
			end

			#Has the solution converged
			resid = fsteadystate(θ, ob.adjmat, ob.pvec, kappa)
			if norm(resid,2) < 1e-6
				return true
			end
		end

		#If the network is still syncing and no edge trips, continue
		return false
	end
	condition(u,t,integrator) = conditiontemp(u,t,integrator,obj1)
	affect!(integrator) = terminate!(integrator)

	#Time step
	f(u, p, t) = fswing(u, obj1.adjmat, P1, 1.0, 1.0, kappa)
	tspan = (0.0, 500.0)
	prob = ODEProblem(f, ψ, tspan)
	cb = DiscreteCallback(condition, affect!)
	sol = solve(prob,reltol=1e-8,abstol=1e-8,save_everystep=false,dense=false,save_end=true,callback=cb)

	#Get the state vector, split into theta and omega
	ψ_1 = sol[:,end]
	ω_1 = ψ_1[1:nodecardinality]
	θ_1 = ψ_1[nodecardinality+1:end]

	#If the solution has not converged, then return
	if obj1.sync == false
		return 0, 0, edgecardinality, 0
	end

	#If the solution has converged with no further overloads, then return
	if (obj1.trip == false) && (obj1.sync == true)
		return edgecardinality, 0, 0, 0
	end

	#If there are overloads, but omega still in acceptable range,
	#then compute number of new connected components and continue
	#recursion on each of them
	#First construct new edgeset and incidence matrix
	overloadcounter = 0
	survivors = Array{Int64, 1}([])
	flow = edgepower(θ_1, E1, kappa)./maxflow
	for i=1:edgecardinality
		ei = edgeset[i]
		if abs(flow[i]) > alpha
			net.edgelist[ei].active = false
			overloadcounter += 1
		else
			push!(survivors, i)
		end
	end
	newedgecard = length(survivors)
	E2 = zeros(Float64, nodecardinality, newedgecard)
	for j=1:newedgecard
		E2[:,j] = E1[:,survivors[j]]
	end

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

	#Clear any superfluous data
	E1 = nothing
	E2 = nothing
	Adj = nothing
	obj1 = nothing

	#Call fracture function on each of the components, passing down info
	descendent_edges = 0
	descendent_overloads = 0
	descendent_desyncs = 0
	descendent_unbalances = 0
	for i=1:components
		chunk = table[i]
		n_subset = length(chunk)
		nodesubset = zeros(Int64, length(chunk))
		ω_subset = zeros(Float64, length(chunk))
		θ_subset = zeros(Float64, length(chunk))
		ψ_subset = zeros(Float64, 2*n_subset)
		for j=1:length(chunk)
			nodesubset[j] = nodeset[chunk[j]]
			ω_subset[j] = ω_1[chunk[j]]
			θ_subset[j] = θ_1[chunk[j]]
		end
		ψ_subset[1:n_subset] = ω_subset
		ψ_subset[n_subset+1:end] = θ_subset
		xz, oz, dz, uz = swingfracturewithbreakdown!(net, ψ_subset, nodesubset, P, synctol, alpha, kappa, maxflow)
		descendent_edges += xz
		descendent_overloads += oz
		descendent_desyncs += dz
		descendent_unbalances += uz
	end
	return descendent_edges, descendent_overloads, descendent_desyncs, descendent_unbalances
end

function swingcascadeandtvalpha(ensemblesize::Int64, n::Int64, ns::Int64, nd::Int64, q::Float64,
								k::Int64, alphares::Int64, alphamin::Float64, alphamax::Float64,
								I::Float64, D::Float64, κ::Float64)

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

	println("Starting")
	for z=1:ensemblesize
		println("Ensemble = ", z)

		net, E = smallworldnetwork(n, k, q)
		# net, E = scalefreenetwork(n, 2*n, 3)
		sources, sinks = sourcesinklocs(ns, nd, n)
		P = sourcesinkvector(sources, sinks, n, 1.0)
		m = numedges(net)
		nodeset = collect(1:n)
		edgeset = collect(1:m)
		A = adjacencymatrix(E)
		ψ = rand(2*n)

		#Time step to find steady state
		f(u, p, t) = fswing(u, A, P, I, D, κ)
		tspan = (0.0, 100.0)
		prob = ODEProblem(f, ψ, tspan)
		sol = solve(prob,reltol=1e-8,abstol=1e-8,save_everystep=false,dense=false,save_end=true)

		#Check for convergence
		ψ = sol[:,end]
		ω = ψ[1:n]
		θ = ψ[n+1:end]
		x = fsteadystate(θ, A, P, κ)
		if norm(x,2) > 1e-3
			println("Warning: no steady state found.")
		end

		#Compute the flow and normalise
		flow_swing = edgepower(θ, E, κ)
		flowmax_swing = maximum(abs.(flow_swing))
		fmaxes[z] = flowmax_swing
		flow_swing = flow_swing./flowmax_swing

		redges = zeros(Float64, alphares)
		rtimes = zeros(Float64, alphares)

		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.(flow_swing))
			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)

			totedgessurviving = 0
			timetakenvec = []
			for i=1:components
				chunk = table[i]
				n_subset = length(chunk)
				nodesubset = zeros(Int64, length(chunk))
				ω_subset = zeros(Float64, length(chunk))
				θ_subset = zeros(Float64, length(chunk))
				ψ_subset = zeros(Float64, 2*n_subset)
				for j=1:length(chunk)
					nodesubset[j] = nodeset[chunk[j]]
					ω_subset[j] = ω[chunk[j]]
					θ_subset[j] = θ[chunk[j]]
				end
				ψ_subset[1:n_subset] = ω_subset
				ψ_subset[n_subset+1:end] = θ_subset
				edgessurviving, timetaken = swingfracturewithtime!(net, ψ_subset, nodesubset, P, 3.0, alpha, κ, flowmax_swing, D)
				totedgessurviving += edgessurviving
				push!(timetakenvec, timetaken)
			end
			redges[ind] = totedgessurviving/m
			rtimes[ind] = maximum(timetakenvec)
		end

		#Find the critical value
		aftertransition = findall(x->x>0.5, redges)
		if length(aftertransition) < 1
			critvalues[z] = NaN
		else
			ac = alphas[aftertransition[1]]
			critvalues[z] = ac
		end

		#Store S
		edgemat[z,:] = redges
		tmat[z,:] = rtimes
	end
	return alphas, mean(edgemat, dims=1)[1,:], mean(tmat, dims=1)[1,:] ,mean(fmaxes), critvalues
end

function svt(ensemblesize::Int64, n::Int64, ns::Int64, nd::Int64, q::Float64,
								k::Int64, alpha::Float64,
								I::Float64, D::Float64, κ::Float64)

	println("Starting")
	vedges = zeros(Float64, ensemblesize)
	vtimes = zeros(Float64, ensemblesize)
	for z=1:ensemblesize
		println("Ensemble = ", z)

		net, E = smallworldnetwork(n, k, q)
		sources, sinks = sourcesinklocs(ns, nd, n)
		P = sourcesinkvector(sources, sinks, n, 1.0)
		m = numedges(net)
		nodeset = collect(1:n)
		edgeset = collect(1:m)
		A = adjacencymatrix(E)
		ψ = rand(2*n)

		#Time step to find steady state
		f(u, p, t) = fswing(u, A, P, I, D, κ)
		tspan = (0.0, 100.0)
		prob = ODEProblem(f, ψ, tspan)
		sol = solve(prob,reltol=1e-8,abstol=1e-8,save_everystep=false,dense=false,save_end=true)

		#Check for convergence
		ψ = sol[:,end]
		ω = ψ[1:n]
		θ = ψ[n+1:end]
		x = fsteadystate(θ, A, P, κ)
		if norm(x,2) > 1e-3
			println("Warning: no steady state found.")
		end

		#Compute the flow and normalise
		flow_swing = edgepower(θ, E, κ)
		flowmax_swing = maximum(abs.(flow_swing))
		flow_swing = flow_swing./flowmax_swing

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

		#Knock out most loaded edge
		d = argmax(abs.(flow_swing))
		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)

		totedgessurviving = 0
		timetakenvec = []
		for i=1:components
			chunk = table[i]
			n_subset = length(chunk)
			nodesubset = zeros(Int64, length(chunk))
			ω_subset = zeros(Float64, length(chunk))
			θ_subset = zeros(Float64, length(chunk))
			ψ_subset = zeros(Float64, 2*n_subset)
			for j=1:length(chunk)
				nodesubset[j] = nodeset[chunk[j]]
				ω_subset[j] = ω[chunk[j]]
				θ_subset[j] = θ[chunk[j]]
			end
			ψ_subset[1:n_subset] = ω_subset
			ψ_subset[n_subset+1:end] = θ_subset
			edgessurviving, timetaken = swingfracturewithtime!(net, ψ_subset, nodesubset, P, 3.0, alpha, κ, flowmax_swing, D)
			totedgessurviving += edgessurviving
			push!(timetakenvec, timetaken)
		end
		vedges[z] = totedgessurviving/m
		vtimes[z] = maximum(timetakenvec)
	end
	return vedges, vtimes
end

function swingcascadebisection(ensemblesize::Int64, tol::Float64, n::Int64, ns::Int64, nd::Int64,
								q::Float64, k::Int64, κ::Float64, I::Float64, D::Float64)

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

		# Construct each network
		net, E = smallworldnetwork(n, k, q)
		# net, E = austriannetwork()
		sources, sinks = sourcesinklocsmod(net, ns, nd, n)
		P = sourcesinkvector(sources, sinks, n, 1.0)
		m = numedges(net)

		#Make the initial state
		A = adjacencymatrix(E)
		ψ = rand(2*n)
		f(u, p, t) = fswing(u, A, P, I, D, κ)
		tspan = (0.0, 250.0)
		prob = ODEProblem(f, ψ, tspan)
		sol = solve(prob,reltol=1e-8,abstol=1e-8,save_everystep=false,dense=false,save_end=true)

		#Check for convergence
		ψ = sol[:,end]
		ω = ψ[1:n]
		θ = ψ[n+1:end]
		x = fsteadystate(θ, A, P, κ)
		if norm(x,2) > 1e-3
			println("Warning: no steady state found.")
		end

		#Compute the flow pattern and normalise
		flow_swing = edgepower(θ, E, κ)
		flowmax_swing = maximum(abs.(flow_swing))
		fmaxes[i] = flowmax_swing
		flow_swing = flow_swing./flowmax_swing

		#Store the ensemble member
		ensemblemember = CascadeNetworkSwing(net, E, sources, sinks, flow_swing, ψ)
		ensemble[i] = ensemblemember
	end

	alpha = 0.01
	stepsize = 0.3
	beneath0 = true
	beneath1 = true
	iterations = 1

	while abs(stepsize) > tol

		# println(alpha)
		runningsurvivors = Array{Float64, 1}(undef, ensemblesize)
		for z=1:ensemblesize

			#Unpack the network
			net = ensemble[z].network
			E = ensemble[z].E
			m = numedges(net)
			sources = ensemble[z].sources
			sinks = ensemble[z].sinks
			F = ensemble[z].initialflow
			ψ = ensemble[z].initialnodestate
			ω = ψ[1:n]
			θ = ψ[n+1:end]
			nodeset = collect(1:n)
			edgeset = collect(1:m)
			P = sourcesinkvector(sources, sinks, n, 1.0)

			#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
			totedgessurviving = 0
			for i=1:components
				chunk = table[i]
				n_subset = length(chunk)
				nodesubset = zeros(Int64, length(chunk))
				ω_subset = zeros(Float64, length(chunk))
				θ_subset = zeros(Float64, length(chunk))
				ψ_subset = zeros(Float64, 2*n_subset)
				for j=1:length(chunk)
					nodesubset[j] = nodeset[chunk[j]]
					ω_subset[j] = ω[chunk[j]]
					θ_subset[j] = θ[chunk[j]]
				end
				ψ_subset[1:n_subset] = ω_subset
				ψ_subset[n_subset+1:end] = θ_subset
				edgessurviving = swingfracture!(net, ψ_subset, nodesubset, P, 3.0, alpha, κ, fmaxes[z])
				totedgessurviving += edgessurviving
			end

			#Record fraction of edges remaining
			runningsurvivors[z] = totedgessurviving/Float64(m)
		end

		#If fraction has crossed the target then reverse and half the step size
		global avremaining = mean(runningsurvivors)
		# println("avremaining", avremaining)
		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
end

function criticalcouplingbisection(ensemblesize::Int64, tol::Float64, n::Int64, ns::Int64, nd::Int64,
								q::Float64, k::Int64, κ::Float64, I::Float64, D::Float64)

	#Iterate over ensemble
	crit = zeros(Float64, ensemblesize)
	κ0 = κ
	for z=1:ensemblesize
		println(z)
		#Make the network
		net, E = austriannetwork()
		sources, sinks = sourcesinklocsmod(net, ns, nd, n)
		P = sourcesinkvector(sources, sinks, n, 1.0)
		m = numedges(net)

		#Make the initial state
		A = adjacencymatrix(E)
		ψlast = rand(2*n)
		ψnew = rand(2*n)
		f(u, p, t) = fswing(u, A, P, I, D, κ)
		tspan = (0.0, 200.0)
		prob = ODEProblem(f, ψlast, tspan)
		sol = solve(prob,reltol=1e-8,abstol=1e-8,save_everystep=false,
			dense=false,save_end=true)
		ψlast = sol[:,end]

		obj1 = NetStore2(A, P, E, κ, 0.0, 0.0, 0.0, true, false)

		tol = 1e-3

		#Iterate over coupling
		stepsize = 0.01
		bisecting = true
		κ = κ0
		κnew = κ
		κold = κ
		while bisecting
			println(κ)

			f(u, p, t) = fswing(u, A, P, I, D, κ)
			obj1.kappa = κ
			κnew = κ
			tspan = (0.0, 100.0)
			prob = ODEProblem(f, ψlast, tspan)

			#Set up callback function
			convergedflag = false
			function conditiontemp(u1, t1, integrator, ob::NetStore2)
				θ = u1[n+1:end]
				if (t1 > 10.0)
					resid = fsteadystate(θ, ob.adjmat, ob.pvec, ob.kappa)
					if norm(resid,2) < 1e-5
						convergedflag = true
						return true
					end
				end
				return false
			end
			condition(u,t,integrator) = conditiontemp(u,t,integrator,obj1)
			affect!(integrator) = terminate!(integrator)
			cb = DiscreteCallback(condition, affect!)
			sol = solve(prob,reltol=1e-8,abstol=1e-8,save_everystep=false,
				dense=false,save_end=true, callback=cb)
			if convergedflag == true
				ψlast = sol[:,end]

				#If bisection done, then stop
				if stepsize < tol
					crit[z] = κ
					bisecting = false
				end

				#Else, continue
				κold = κ
				κ = κ - stepsize
			else
				stepsize = stepsize/2.0
				κ = κold - stepsize
			end
		end
	end
	return crit

end

function swingcascadeouterbisection(ensemblesize::Int64, tol::Float64, n::Int64, ns::Int64, nd::Int64,
								q::Float64, k::Int64, κ::Float64, I::Float64, D::Float64)

	rhos = zeros(Float64, ensemblesize)
	for z=1:ensemblesize
		println("Ensemble ", z)
		net, E = smallworldnetwork(n, k, q)
		sources, sinks = sourcesinklocsmod(net, ns, nd, n)
		P = sourcesinkvector(sources, sinks, n, 1.0)
		m = numedges(net)

		#Make the initial state
		A = adjacencymatrix(E)
		ψ = rand(2*n)
		f(u, p, t) = fswing(u, A, P, I, D, κ)
		tspan = (0.0, 250.0)
		prob = ODEProblem(f, ψ, tspan)
		sol = solve(prob,reltol=1e-8,abstol=1e-8,save_everystep=false,dense=false,save_end=true)

		#Check for convergence
		alpha = 0.01
		ψ = sol[:,end]
		ω = ψ[1:n]
		θ = ψ[n+1:end]
		x = fsteadystate(θ, A, P, κ)
		if norm(x,2) > 1e-3
			println("Warning: no steady state found.")
			alpha = NaN
		else

			#Compute the flow pattern and normalise
			F = edgepower(θ, E, κ)
			fmax = maximum(abs.(F))
			F = F./fmax

			#Now do bisection
			stepsize = 0.3
			beneath0 = true
			beneath1 = true
			iterations = 1
			while abs(stepsize) > tol

				#Reactivate all edges
				nodeset = collect(1:n)
				edgeset = collect(1:m)
				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
				totedgessurviving = 0
				for i=1:components
					chunk = table[i]
					n_subset = length(chunk)
					nodesubset = zeros(Int64, length(chunk))
					ω_subset = zeros(Float64, length(chunk))
					θ_subset = zeros(Float64, length(chunk))
					ψ_subset = zeros(Float64, 2*n_subset)
					for j=1:length(chunk)
						nodesubset[j] = nodeset[chunk[j]]
						ω_subset[j] = ω[chunk[j]]
						θ_subset[j] = θ[chunk[j]]
					end
					ψ_subset[1:n_subset] = ω_subset
					ψ_subset[n_subset+1:end] = θ_subset
					edgessurviving = swingfracture!(net, ψ_subset, nodesubset, P, 2.0, alpha, κ, fmax)
					totedgessurviving += edgessurviving
				end
				S = totedgessurviving/Float64(m)

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

			end
		end
		rhos[z] = alpha
	end
	return rhos

end

function swingcascadebreakdownvsalpha(ensemblesize::Int64, n::Int64, ns::Int64, nd::Int64, q::Float64,
								k::Int64, alphares::Int64, alphamin::Float64, alphamax::Float64,
								I::Float64, D::Float64, κ::Float64)

	alphas = collect(range(alphamin, length=alphares, stop=alphamax))
	edgemat = zeros(Float64, ensemblesize, alphares)
	overloadmat = zeros(Float64, ensemblesize, alphares)
	desyncmat = zeros(Float64, ensemblesize, alphares)
	unbalancemat = zeros(Float64, ensemblesize, alphares)
	fmaxes = zeros(Float64, ensemblesize)
	critvalues = zeros(Float64, ensemblesize)

	println("Starting")
	for z=1:ensemblesize
		println("Ensemble = ", z)

		net, E = smallworldnetwork(n, k, q)
		# net, E = scalefreenetwork(n, 2*n, 3)
		sources, sinks = sourcesinklocs(ns, nd, n)
		P = sourcesinkvector(sources, sinks, n, 1.0)
		m = numedges(net)
		nodeset = collect(1:n)
		edgeset = collect(1:m)
		A = adjacencymatrix(E)
		ψ = rand(2*n)

		#Time step to find steady state
		f(u, p, t) = fswing(u, A, P, I, D, κ)
		tspan = (0.0, 100.0)
		prob = ODEProblem(f, ψ, tspan)
		sol = solve(prob,reltol=1e-8,abstol=1e-8,save_everystep=false,dense=false,save_end=true)

		#Check for convergence
		ψ = sol[:,end]
		ω = ψ[1:n]
		θ = ψ[n+1:end]
		x = fsteadystate(θ, A, P, κ)
		if norm(x,2) > 1e-3
			println("Warning: no steady state found.")
		end

		#Compute the flow and normalise
		flow_swing = edgepower(θ, E, κ)
		flowmax_swing = maximum(abs.(flow_swing))
		fmaxes[z] = flowmax_swing
		flow_swing = flow_swing./flowmax_swing

		redges = zeros(Float64, alphares)
		roverloads = zeros(Float64, alphares)
		rdesyncs = zeros(Float64, alphares)
		runbalances = zeros(Float64, alphares)

		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.(flow_swing))
			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)

			totedgessurviving = 0
			totaloverloads = 0
			totaldesyncs = 0
			totalunbalances = 0
			for i=1:components
				chunk = table[i]
				n_subset = length(chunk)
				nodesubset = zeros(Int64, length(chunk))
				ω_subset = zeros(Float64, length(chunk))
				θ_subset = zeros(Float64, length(chunk))
				ψ_subset = zeros(Float64, 2*n_subset)
				for j=1:length(chunk)
					nodesubset[j] = nodeset[chunk[j]]
					ω_subset[j] = ω[chunk[j]]
					θ_subset[j] = θ[chunk[j]]
				end
				ψ_subset[1:n_subset] = ω_subset
				ψ_subset[n_subset+1:end] = θ_subset
				edgessurviving, overloads, desyncs, unbalances = swingfracturewithbreakdown!(net, ψ_subset, nodesubset, P, 3.0, alpha, κ, flowmax_swing)
				totedgessurviving += edgessurviving
				totaloverloads += overloads
				totaldesyncs += desyncs
				totalunbalances += unbalances
			end
			redges[ind] = totedgessurviving/m
			roverloads[ind] = (m - totedgessurviving - totaldesyncs - totalunbalances)/m
			rdesyncs[ind] = totaldesyncs/m
			runbalances[ind] = totalunbalances/m
		end

		#Find the critical value
		aftertransition = findall(x->x>0.5, redges)
		if length(aftertransition) < 1
			critvalues[z] = NaN
		else
			ac = alphas[aftertransition[1]]
			critvalues[z] = ac
		end

		#Store S
		edgemat[z,:] = redges
		overloadmat[z,:] = roverloads
		desyncmat[z,:] = rdesyncs
		unbalancemat[z,:] = runbalances
	end
	return alphas, mean(edgemat, dims=1)[1,:], mean(overloadmat, dims=1)[1,:], mean(desyncmat, dims=1)[1,:], mean(unbalancemat, dims=1)[1,:], mean(fmaxes), critvalues
end

max_stars_repo_name: o-smith/GridResilience

max_issues_repo_name: o-smith/GridResilience

max_forks_repo_name: o-smith/GridResilience


text:
# include("dependencies_for_runtests.jl")

using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBoundary
using Oceananigans.ImmersedBoundaries: conditional_length
using Statistics: mean, mean!, norm
using CUDA: @allowscalar

@testset "Conditional Reductions" begin
    for arch in archs
        @info "    Testing Reductions on Immersed fields"

        grid = RectilinearGrid(arch, size = (6, 1, 1), extent = (1, 1, 1))

        ibg  = ImmersedBoundaryGrid(grid, GridFittedBoundary((x, y, z) -> (x < 0.5)))

        fful = Field{Center, Center, Center}(grid)
        fimm = Field{Center, Center, Center}(ibg)

        @test conditional_length(fimm) == length(fimm) / 2

        fful .= 2
        fimm .= 2

        @allowscalar fimm[1, :, :] .= 1e6
        @allowscalar fimm[2, :, :] .= -1e4
        @allowscalar fimm[3, :, :] .= -12.5

        @test norm(fful) ≈ √2 * norm(fimm)
        @test mean(fful) ≈      mean(fimm)

        for reduc in (mean, maximum, minimum)
            @test reduc(fful) == reduc(fimm)
            @test all(Array(interior(reduc(fful, dims=1)) .== interior(reduc(fimm, dims=1))))
        end
        @test sum(fful) == sum(fimm) * 2
        @test all(Array(interior(sum(fful, dims=1)) .== interior(sum(fimm, dims=1)) .* 2))

        @test prod(fful) == prod(fimm) * 8
        @test all(Array(interior(prod(fful, dims=1)) .== interior(prod(fimm, dims=1)) .* 8))

        @info "    Testing Reductions in Standard fields"

        fcon = Field{Center, Center, Center}(grid)

        fcon .= 2

        @allowscalar fcon[1, :, :] .= 1e6
        @allowscalar fcon[2, :, :] .= -1e4
        @allowscalar fcon[3, :, :] .= -12.5

        @test norm(fful) ≈ √2 * norm(fcon, condition = (i, j, k, x, y) -> i > 3)

        for reduc in (mean, maximum, minimum)
            @test reduc(fful) == reduc(fcon, condition = (i, j, k, x, y) -> i > 3)
            @test all(Array(interior(reduc(fful, dims=1)) .== interior(reduc(fcon, condition = (i, j, k, x, y) -> i > 3, dims=1))))
        end
        @test sum(fful) == sum(fcon, condition = (i, j, k, x, y) -> i > 3) * 2
        @test all(Array(interior(sum(fful, dims=1)) .== interior(sum(fcon, condition = (i, j, k, x, y) -> i > 3, dims=1)) .* 2))

        @test prod(fful) == prod(fcon, condition = (i, j, k, x, y) -> i > 3) * 8
        @test all(Array(interior(prod(fful, dims=1)) .== interior(prod(fcon, condition = (i, j, k, x, y) -> i > 3, dims=1)) .* 8))

        @info "    Testing in-place conditional reductions"

        redimm = Field{Nothing, Center, Center}(ibg)
        for (reduc, reduc!) in zip((mean, maximum, minimum, sum, prod), (mean!, maximum!, minimum!, sum!, prod!))
            @test reduc!(redimm, fimm)[1, 1 , 1] == reduc(fcon, condition = (i, j, k, x, y) -> i > 3, dims = 1)[1, 1, 1]
        end
    end
end


text:
# Todo, check all packages for what they check isdefined on
using Tokenize.Tokens: COLON

# Lower commit introduced, upper commit removed
const DEFINES = Dict{Symbol, NTuple{2, VersionNumber}}()

# When a symbol changed name
function add_deprecation(from::Symbol, to::Symbol, v::VersionNumber)
    push!(DEFINES, from => (v"0.0", v))
    push!(DEFINES, to   => (v, typemax(VersionNumber)))
end

# When a symbol was removed
add_removal(sym::Symbol, v::VersionNumber) = push!(DEFINES, sym => (v"0.0", v))

# When a symbol was added
add_addition(sym::Symbol, v::VersionNumber) = push!(DEFINES, sym => (v, typemax(VersionNumber)))

let
    # STDLIB
    add_removal(:Test, v"0.7.0-DEV.2004")

    add_addition(:error_color, v"0.6.0-dev.1545") # 18628
    add_addition(:unwrap_unionall, v"0.6.0-dev.2123") # 18457
    add_addition(:xor, v"0.6.0-dev.1236") # 18977

    add_addition(Symbol("@nospecialize"), v"0.7.0-DEV.969")

    add_deprecation(:isleaftype, :isconcrete, v"0.7.0-DEV.1775") # 23666
    add_addition(:__precompile__, v"0.4.0-dev+6521")
end

begin
    struct IsDefinedBaseCheck; vers; end
    register(IsDefinedBaseCheck, Deprecation(
        "This isdefined check for Base is always true or false for the versions supported by this package",
        "julia",
        typemin(VersionNumber), typemin(VersionNumber), typemax(VersionNumber)))
        IsDefinedBaseCheck() = IsDefinedBaseCheck(Pkg.Reqs.parse(IOBuffer("julia 0.6"))["julia"])

    function dep_for_vers(::Type{IsDefinedBaseCheck}, vers)
        IsDefinedBaseCheck(vers["julia"])
    end

  get_symbol(x) = nothing
  function get_symbol(x::EXPR) #, id)
        c = children(x)

        if isexpr(x, CSTParser.Quotenode) && length(c) == 2 &&
           isexpr(c[1], OPERATOR, Tokens.COLON) && isexpr(c[2], CSTParser.IDENTIFIER)
           return c[2].val
        end

        if isexpr(x, CSTParser.Call) && isexpr(c[1], CSTParser.IDENTIFIER) &&
                c[1].val == "Symbol" && isexpr(c[3], CSTParser.LITERAL)
            return c[3].val
        end

        return nothing
    end

    match(IsDefinedBaseCheck, CSTParser.If) do x
        dep, expr, resolutions, context = x
        resolve_isdefined_call(children(expr)[2], dep, expr, resolutions, context)
    end

    match(IsDefinedBaseCheck, CSTParser.BinarySyntaxOpCall) do x
        dep, expr, resolutions, context = x
        !isexpr(children(expr)[2], OPERATOR, Tokens.LAZY_AND) && return
        resolve_isdefined_call(children(expr)[1], dep, expr, resolutions, context)
    end

    function resolve_isdefined_call(call, dep, expr, resolutions, context)
        call = call.expr

        isnegated = false
        if isexpr(call, CSTParser.UnaryOpCall)
            if isexpr(children(call)[1], OPERATOR, Tokens.NOT)
                call = children(call)[2]
                isnegated = true
            else
                return
            end
        end

        n_args = (length(children(call)) - 1  # function name
                                         - 3) # punctuation ( , )
        !(n_args == 2) && return
        function_identifier = children(call)[1]
        !is_identifier(function_identifier, "isdefined") && return
        first_arg, second_arg = children(call)[3], children(call)[5]
        symval = get_symbol(second_arg)
        symval == nothing && return
        !is_identifier(first_arg, "Base") && return
        sym = Symbol(symval)

        # Also applies in @static context, but not necessarily in other macro contexts
        replace_expr = expr
        !haskey(DEFINES, sym) && return
        if context.in_macrocall
            context.top_macrocall == parent(expr) || return
            is_macroname(context.top_macrocall, "static") || return
            replace_expr = context.top_macrocall
        end

        baselow, basehigh = DEFINES[sym]

        # bl ----------------- bh
        #     il ------- ih
        alwaysdefined = all(interval -> (interval.lower >= baselow && (interval.upper < basehigh || basehigh == typemax(VersionNumber))), dep.vers.intervals)

        # bl ---- bh                   bl ---- bh
        #             il ------- ih
        neverdefined = all(interval -> (interval.lower >= basehigh || interval.upper < baselow), dep.vers.intervals)
        if isnegated
            alwaysdefined, neverdefined = neverdefined, alwaysdefined
        end
        @assert !(alwaysdefined && neverdefined)
        if alwaysdefined
            resolve_inline_body(dep, resolutions, expr, replace_expr)
        elseif neverdefined
            resolve_delete_expr(dep, resolutions, expr, replace_expr)
        end
    end
end


text:
using Test
using InteractiveUtils

const_int() = 1

const lno = LineNumberNode(1, :none)

let ci = @code_lowered const_int()
    @eval function oc_trivial()
        $(Expr(:new_opaque_closure, Tuple{}, false, Any, Any,
            Expr(:opaque_closure_method, nothing, 0, lno, ci)))
    end
end
@test isa(oc_trivial(), Core.OpaqueClosure{Tuple{}, Any})
@test oc_trivial()() == 1

let ci = @code_lowered const_int()
    @eval function oc_simple_inf()
        $(Expr(:new_opaque_closure, Tuple{}, false, Union{}, Any,
            Expr(:opaque_closure_method, nothing, 0, lno, ci)))
    end
end
@test isa(oc_simple_inf(), Core.OpaqueClosure{Tuple{}, Int})
@test oc_simple_inf()() == 1

struct OcClos2Int
    a::Int
    b::Int
end
(a::OcClos2Int)() = getfield(a, 1) + getfield(a, 2)
let ci = @code_lowered OcClos2Int(1, 2)();
    @eval function oc_trivial_clos()
        $(Expr(:new_opaque_closure, Tuple{}, false, Int, Int,
            Expr(:opaque_closure_method, nothing, 0, lno, ci),
            1, 2))
    end
end
@test oc_trivial_clos()() == 3

let ci = @code_lowered OcClos2Int(1, 2)();
    @eval function oc_self_call_clos()
        $(Expr(:new_opaque_closure, Tuple{}, false, Int, Int,
            Expr(:opaque_closure_method, nothing, 0, lno, ci),
            1, 2))()
    end
end
@test @inferred(oc_self_call_clos()) == 3
let opt = @code_typed oc_self_call_clos()
    @test length(opt[1].code) == 1
    @test isa(opt[1].code[1], Core.ReturnNode)
end

struct OcClos1Any
    a
end
(a::OcClos1Any)() = getfield(a, 1)
let ci = @code_lowered OcClos1Any(1)()
    @eval function oc_pass_clos(x)
        $(Expr(:new_opaque_closure, Tuple{}, false, Any, Any,
            Expr(:opaque_closure_method, nothing, 0, lno, ci),
            :x))
    end
end
@test oc_pass_clos(1)() == 1
@test oc_pass_clos("a")() == "a"

let ci = @code_lowered OcClos1Any(1)()
    @eval function oc_infer_pass_clos(x)
        $(Expr(:new_opaque_closure, Tuple{}, false, Union{}, Any,
            Expr(:opaque_closure_method, nothing, 0, lno, ci),
            :x))
    end
end
@test isa(oc_infer_pass_clos(1), Core.OpaqueClosure{Tuple{}, typeof(1)})
@test isa(oc_infer_pass_clos("a"), Core.OpaqueClosure{Tuple{}, typeof("a")})
@test oc_infer_pass_clos(1)() == 1
@test oc_infer_pass_clos("a")() == "a"

let ci = @code_lowered identity(1)
    @eval function oc_infer_pass_id()
        $(Expr(:new_opaque_closure, Tuple{Any}, false, Any, Any,
            Expr(:opaque_closure_method, nothing, 1, lno, ci)))
    end
end
function complicated_identity(x)
    oc_infer_pass_id()(x)
end
@test @inferred(complicated_identity(1)) == 1
@test @inferred(complicated_identity("a")) == "a"
let ci = (@code_typed complicated_identity(1))[1]
    @test length(ci.code) == 1
    @test isa(ci.code[1], Core.ReturnNode)
end

struct OcOpt
    A
end

(A::OcOpt)() = ndims(getfield(A, 1))

let ci = @code_lowered OcOpt([1 2])()
    @eval function oc_opt_ndims(A)
        $(Expr(:new_opaque_closure, Tuple{}, false, Union{}, Any,
            Expr(:opaque_closure_method, nothing, 0, lno, ci),
            :A))
    end
end
oc_opt_ndims([1 2])

let A = [1 2]
    let Oc = oc_opt_ndims(A)
        @test_broken sizeof(Oc.env) == 0
        @test Oc() == 2
    end
end

using Base.Experimental: @opaque

@test @opaque(x->2x)(8) == 16
let f = @opaque (x::Int, y::Float64)->(2x, 3y)
    @test_throws TypeError f(1, 1)
    @test f(2, 3.0) === (4, 9.0)
end
function uses_frontend_opaque(x)
    @opaque y->x+y
end
@test uses_frontend_opaque(10)(8) == 18

# World age mechanism
function test_oc_world_age end
mk_oc_world_age() = @opaque ()->test_oc_world_age()
g_world_age = @opaque ()->test_oc_world_age()
h_world_age = mk_oc_world_age()
@test isa(h_world_age, Core.OpaqueClosure{Tuple{}, Union{}})
test_oc_world_age() = 1
@test_throws MethodError g_world_age()
@test_throws MethodError h_world_age()
@test mk_oc_world_age()() == 1
g_world_age = @opaque ()->test_oc_world_age()
@test g_world_age() == 1
@test isa(mk_oc_world_age(), Core.OpaqueClosure{Tuple{}, Int})

# Evil, dynamic Vararg stuff (don't do this - made to work for consistency)
function maybe_opaque(isva::Bool)
    T = isva ? Vararg{Int, 1} : Int
    @opaque (x::T)->x
end
@test maybe_opaque(false)(1) == 1
@test maybe_opaque(true)(1) == (1,)

# Vargarg in complied mode
mk_va_opaque() = @opaque (x...)->x
@test mk_va_opaque()(1) == (1,)
@test mk_va_opaque()(1,2) == (1,2)

# OpaqueClosure show method
@test repr(@opaque x->1) == "(::Any)::Any->◌"

# Opaque closure in CodeInfo returned from generated functions
function mk_ocg(args...)
    ci = @code_lowered const_int()
    cig = Meta.lower(@__MODULE__, Expr(:new_opaque_closure, Tuple{}, false, Any, Any,
        Expr(:opaque_closure_method, nothing, 0, lno, ci))).args[1]
    cig.slotnames = Symbol[Symbol("#self#")]
    cig.slottypes = Any[Any]
    cig.slotflags = UInt8[0x00]
    cig
end

@eval function oc_trivial_generated()
    $(Expr(:meta, :generated_only))
    $(Expr(:meta,
            :generated,
            Expr(:new,
                Core.GeneratedFunctionStub,
                :mk_ocg,
                Any[:oc_trivial_generated],
                Any[],
                @__LINE__,
                QuoteNode(Symbol(@__FILE__)),
                true)))
end
@test isa(oc_trivial_generated(), Core.OpaqueClosure{Tuple{}, Any})
@test oc_trivial_generated()() == 1

# Constprop through varargs OpaqueClosure
function oc_varargs_constprop()
    oc = @opaque (args...)->args[1]+args[2]+args[3]
    return Val{oc(1,2,3)}()
end
@test Base.return_types(oc_varargs_constprop, Tuple{}) == Any[Val{6}]

# OpaqueClosure ABI
f_oc_noinline(x) = @opaque function (y)
    @Base._noinline_meta
    x + y
end

let oc = Base.inferencebarrier(f_oc_noinline(1))
    @test oc(2) == 3
end

function f_oc_noinline_call(x, y)
    return f_oc_noinline(x)(y)
end
@test f_oc_noinline_call(1, 2) == 3


text:
function buildnet()
	net = SpikingLayersNet()

	# create layers
	P_L = add_layer!(net, SpikingLayer(SimpleSpikingNeuron,"P_L",100))
	P_M = add_layer!(net, SpikingLayer(SimpleSpikingNeuron,"P_M",100))
	P_S = add_layer!(net, SpikingLayer(SimpleSpikingNeuron,"P_S",100))

	LGN_L = add_layer!(net, SpikingLayer(HHNeuron,"LGN_L",900))
	LGN_C1 = add_layer!(net, SpikingLayer(HHNeuron,"LGN_C1",900))
	LGN_C2 = add_layer!(net, SpikingLayer(HHNeuron,"LGN_C2",900))

	V1_4 = add_layer!(net, SpikingLayer(HHNeuron,"V1_4",900))
	V1_23 = add_layer!(net, SpikingLayer(HHNeuron,"V1_23",900))
	V1_5 = add_layer!(net, SpikingLayer(HHNeuron,"V1_5",900))

	# wire it up!
	# L Photoreceptor -> LGN
	w = 0.005
	#w = 3. # ridiculous weights - just to make sure it will spike

	# L Photoreceptor -> LGN
	for target in ["LGN_L", "LGN_C1", "LGN_C2"]
		c = p_lgn_connectivity(P_L.N, LGN_L.N, EXC, w, false)
		connect!(net, "P_L", target, c)
		println(length(c.synapses))
	end

	# M Photoreceptor -> LGN
	for target in ["LGN_L", "LGN_C2"]
		c = p_lgn_connectivity(P_M.N, LGN_L.N, EXC, w, false) # be careful about this kind of code. assumes LGN are same size!
		connect!(net, "P_M", target, c)
		println(length(c.synapses))
	end

	# M->LGN is inhibitory
	c = p_lgn_connectivity(P_M.N, LGN_C1.N, INH, w, false)
	connect!(net, "P_M", "LGN_C1", c)
	println(length(c.synapses))

	# S Photoreceptor -> LGN
	c = p_lgn_connectivity(P_S.N, LGN_C1.N, INH, w, false)
	connect!(net, "P_S", "LGN_C1", c)
	println(length(c.synapses))

	# LGN->V1 cortex
	c = forward_connectivity(LGN_L.N, V1_4.N, EXC, w, true)
	connect!(net, "LGN_L",  "V1_4", c)
	println(length(c.synapses))

	c = forward_connectivity(LGN_C1.N, V1_4.N, EXC, w, true)
	connect!(net, "LGN_C1", "V1_4", c)
	println(length(c.synapses))

	c = forward_connectivity(LGN_C2.N, V1_23.N, EXC, w*1.5, true)
	connect!(net, "LGN_C2", "V1_23", c)
	println(length(c.synapses))

	# V1->V1
	c = forward_connectivity(V1_4.N, V1_23.N, EXC, w*1.5, true)
	connect!(net, "V1_4", "V1_23", c)
	println(length(c.synapses))

	c = forward_connectivity(V1_23.N, V1_5.N, EXC, w*3, true)
	connect!(net, "V1_23", "V1_5", c)
	println(length(c.synapses))

	# lateral connectivities
	w_pot = w * 0.9; w_dep = w * 1.5 * 1.7;
	c = lateral_connectivity(V1_4.N, w_pot, w_dep)
	connect!(net, "V1_4", "V1_4", c)
	println(length(c.synapses))

	w_pot = w * 0.9; w_dep = w * 1.5 * 1.7; # weights are wrong
	c = lateral_connectivity(V1_4.N, w_pot, w_dep) # replace with V1_23!
	connect!(net, "V1_23", "V1_23", c)
	println(length(c.synapses))

	w_pot = w * 0.9; w_dep = w * 1.5 * 1.7;
	c = lateral_connectivity(V1_4.N, w_pot, w_dep)
	connect!(net, "V1_5", "V1_5", c)
	println(length(c.synapses))

	net
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using LinearAlgebra, Random, StableRNGs
using DensityInterface, ValueShapes
using ArraysOfArrays, Distributions, PDMats, StatsBase

struct _TestDensityStruct{T} <: AbstractDensity
    mvn::T
end
@inline DensityInterface.DensityKind(::_TestDensityStruct) = IsDensity()
DensityInterface.logdensityof(density::_TestDensityStruct, v::Any) = Distributions.logpdf(density.mvn, v)
ValueShapes.totalndof(td::_TestDensityStruct) = Int(3)
BAT.sampler(td::_TestDensityStruct) = BAT.sampler(td.mvn)

struct _UniformDensityStruct{T} <: AbstractDensity
    mvu::T
end
@inline DensityInterface.DensityKind(::_UniformDensityStruct) = IsDensity()
DensityInterface.logdensityof(ud::_UniformDensityStruct, v::Any) = logpdf(ud.mvu, v)
ValueShapes.varshape(ud::_UniformDensityStruct) = varshape(ud.mvu)
ValueShapes.totalndof(ud::_UniformDensityStruct) = Int(3)
BAT.var_bounds(ud::_UniformDensityStruct) = BAT.HyperRectBounds(BAT.HyperRectVolume(zeros(3), ones(3)))

struct _DeltaDensityStruct{T} <: AbstractDensity
    value::T
end
@inline DensityInterface.DensityKind(::_DeltaDensityStruct) = IsDensity()
DensityInterface.logdensityof(dd::_DeltaDensityStruct, v::Any) = (v == dd.value) ? Inf : -Inf
ValueShapes.totalndof(dd::_DeltaDensityStruct) = Int(1)

struct _ShapeDensityStruct{T} <: AbstractDensity
    ntdist::T
end
@inline DensityInterface.DensityKind(::_ShapeDensityStruct) = IsDensity()
DensityInterface.logdensityof(sd::_ShapeDensityStruct, v) = logpdf(sd.ntdist, v)
ValueShapes.varshape(sd::_ShapeDensityStruct) = varshape(sd.ntdist)

struct _NonBATDensity end
@inline DensityInterface.DensityKind(::_NonBATDensity) = IsDensity()
DensityInterface.logdensityof(d::_NonBATDensity, v) = log(norm(v)^2)
ValueShapes.varshape(d::_NonBATDensity) = ArrayShape{Real}(2)

@testset "abstract_density" begin
    mvn = MvNormal(ones(3), PDMat(Matrix{Float64}(I,3,3)))
    td = _TestDensityStruct(mvn)
    tds = BAT.DensityWithShape(td, ArrayShape{Real}(3))

    @test @inferred(isequal(@inferred(BAT.var_bounds(td)), missing))
    @test @inferred(isequal(@inferred(varshape(td)), missing))

    x = rand(3)
    DensityInterface.test_density_interface(tds, x, logpdf(mvn, x))
    @test_throws ArgumentError BAT.checked_logdensityof(td, [Inf, Inf, Inf])
    @test_throws BAT.DensityEvalException BAT.checked_logdensityof(tds, [Inf, Inf, Inf])
    @test_throws BAT.DensityEvalException BAT.checked_logdensityof(tds, [Inf, Inf, Inf])
    @test_throws BAT.DensityEvalException BAT.checked_logdensityof(tds, [NaN, NaN, NaN])
    @test_throws ArgumentError BAT.checked_logdensityof(tds, rand(length(mvn)+1))
    @test_throws ArgumentError BAT.checked_logdensityof(tds, rand(length(mvn)-1))

    @test @inferred(BAT.checked_logdensityof(tds, x)) == @inferred(logpdf(mvn, x))

    x = [-Inf, 0, Inf]
    ud_shape_1 = NamedTupleShape(a=ArrayShape{Real}(1), b=ArrayShape{Real}(1), c=ArrayShape{Real}(1))
    ud_shape_2 = NamedTupleShape(a=ArrayShape{Real}(3))
    mvu = product_distribution([Uniform() for i in 1:3])
    ud = _UniformDensityStruct(mvu)

    DensityInterface.test_density_interface(ud, x, -Inf)
    @test @inferred(BAT.checked_logdensityof(ud, x)) == -Inf
    DensityInterface.test_density_interface(ud_shape_1(ud), ud_shape_1(x), logpdf(mvu, x))
    @test @inferred(BAT.checked_logdensityof(ud_shape_1(ud), ud_shape_1(x))) == @inferred(logpdf(mvu, x))
    @test @inferred(BAT.checked_logdensityof(ud_shape_2(ud), ud_shape_2(x))) == @inferred(logpdf(mvu, x))

    @test_throws ArgumentError BAT.checked_logdensityof(ud, vcat(x,x))

    @test BAT.checked_logdensityof(ud, x .- eps(1.0)) == -Inf

    @test_throws ArgumentError @inferred(BAT.checked_logdensityof(ud, [0 0 0]))

    ntshape = NamedTupleShape(a=ScalarShape{Real}(), b=ScalarShape{Real}(), c=ScalarShape{Real}())
    shapedasnt = ShapedAsNT(x, ntshape)

    dd = _DeltaDensityStruct(0)
    dds = BAT.DensityWithShape(dd, ScalarShape{Real}())
    DensityInterface.test_density_interface(dd, 0, Inf)
    @test_throws BAT.DensityEvalException BAT.checked_logdensityof(dds, 0)
    DensityInterface.test_density_interface(dds, 0, Inf)

    ntdist = NamedTupleDist(a=mvn, b=mvu)

    x1_for_sd = rand(length(ntdist.a))
    x2_for_sd = rand(length(ntdist.b))
    x_for_sd = vcat(x1_for_sd, x2_for_sd)

    correct_shape_of_sd = NamedTupleShape(a=ArrayShape{Real}(length(ntdist.a)), b=ArrayShape{Real}(length(ntdist.b)))
    incorrect_shape_of_sd = NamedTupleShape(a=ArrayShape{Real}(length(ntdist.a)-1), b=ArrayShape{Real}(length(ntdist.b)+1))

    x_for_sd_good_shape = correct_shape_of_sd(x_for_sd)
    x_for_sd_bad_shape = incorrect_shape_of_sd(x_for_sd)

    sd = _ShapeDensityStruct(ntdist)

    DensityInterface.test_density_interface(sd, x_for_sd_good_shape, logpdf(mvn, x1_for_sd) + logpdf(mvu, x2_for_sd))
    @test BAT.checked_logdensityof(sd, x_for_sd_good_shape) == logpdf(mvn, x1_for_sd) + logpdf(mvu, x2_for_sd)
    @test_throws BAT.DensityEvalException BAT.checked_logdensityof(sd, x_for_sd_bad_shape)

    @testset "rand" begin
        td = _TestDensityStruct(mvn)
        @test rand(StableRNG(7002), BAT.bat_sampler(td)) ≈ [2.386799038, 1.072161895, 0.791486531]
    end

    @testset "non-BAT densities" begin
        d = _NonBATDensity()
        x = randn(3)
        @test @inferred(convert(AbstractDensity, d)) isa BAT.WrappedNonBATDensity
        bd = convert(AbstractDensity, d)
        DensityInterface.test_density_interface(bd, x, logdensityof(d, x))
        @test @inferred(logdensityof(bd, x)) == logdensityof(d, x)
        @test @inferred(logdensityof(bd)) == logdensityof(d)
        @test @inferred(varshape(bd)) == varshape(d)
    end
end


text:
using Meganet
using Base.Test

# network
nf = 18
nc = 3
nex = 4
for TYPE=[Float64,Float32]
K = getDenseKernel(TYPE,[nf,nf])
Bout = randn(TYPE,nf,3)
nLayer = getBatchNormLayer(TYPE,[div(nf,nc),nc],isTrainable=true)
L     = getSingleLayer(TYPE,K,nLayer,Bout=Bout)
net = getResNN(TYPE,L,4,.1)

# regularizers
pRegTh = getTikhonovReg(TYPE)
pRegW = getTikhonovReg(TYPE)
# loss
pLoss = getSoftMaxLoss(TYPE)
# data
Y   = randn(TYPE,nf,nex)/10
nw2 = (pLoss.addBias) ? nf+1 : nf
C   = full(sparse(full(rand(1:nc,nex)),collect(1:nex),ones(nex),nc,nex))
C = convert.(TYPE,C)
#initialize
theta = initTheta(net)
W     = randn(TYPE,nc,nw2)/10

objFun = dnnObjFctn(net,pLoss,pRegTh,pRegW)

@testset "dThLoss $TYPE" begin
function testdThLoss(x,v=nothing)
        F,his,dF, = getMisfit(objFun,x,W,Y,C,[],true)
        if v!==nothing
            return F,dot(dF,v)
        else
            return F
        end
end
ckDer, = checkDerivative(testdThLoss,theta,out=false)
@test ckDer
end

@testset "dWLoss $TYPE" begin
function testdWLoss(x,v=nothing)
        F,his,dFth,dF = getMisfit(objFun,theta,x,Y,C,[],true)
        if v!==nothing
            return F,dot(dF,v)
        else
            return F
        end
end
chkDer, = checkDerivative(testdWLoss,vec(W),out=false)
@test chkDer
end

@testset "dJ $TYPE" begin
function testdJ(x,v=nothing)
        F,his,dF = evalObjFctn(objFun,x,Y,C,[],true)
        if v!==nothing
            return F,dot(dF,v)
        else
            return F
        end
end
chkDer, = checkDerivative(testdJ,[vec(theta);vec(W)],out=false)
@test chkDer
end
end


text:
"""
The barrel vault (Scordelis-Lo) roof is one of the benchmarks for linear elastic
analysis of shells.

The candidate element's usefulness in irregular geometries (and most practical
cases involve a high degree of geometric irregularity) is tested. As would be
expected,the irregular mesh results are not as good as those provided by a
regular meshwith the same number of variables.

Problem description

The physical basis of the problem is a deeply arched roof supported only
bydiaphragms at its curved edges (an aircraft hanger), deforming under its own
weight. It is interesting to observe that the geometry is such that the
centerpoint of the roof moves upward under the self-weight(downwardly directed)
load. Perhaps this is one reason why the problem is not straightforward
numerically.
"""
module scordelis_lo_examples

using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FESetShellQ4Module: FESetShellQ4
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json

function _execute_model(n = 8, visualize = true)
    # analytical solution for the vertical deflection and the midpoint of the
    # free edge
    analyt_sol=-0.3024;
    # Parameters:
    E=4.32e8;
    nu=0.0;
    thickness = 0.25; # geometrical dimensions are in feet
    R = 25.0;
    L = 50.0;
    formul = FEMMShellT3FFModule

    tolerance = R/n/1000
    fens, fes = T3block(40/360*2*pi,L/2,n,n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a), y, R*(cos(a)-1))
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    femm.transv_shear_formulation = formul.__TRANSV_SHEAR_FORMULATION_AVERAGE_B
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # rigid diaphragm
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    for i in [1,3,5]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Y
    l1 = selectnode(fens; box = Float64[-Inf Inf L/2 L/2 -Inf Inf], inflate = tolerance)
    for i in [2,4,6]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to X
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,5,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Midpoint of the free edge
    nl = selectnode(fens; box = Float64[sin(40/360*2*pi)*25 sin(40/360*2*pi)*25 L/2 L/2 -Inf Inf], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt[0, 0, -90, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    resultpercent =   dchi.values[nl, 3][1]/analyt_sol*100
    @info "Solution for $(count(fens)*6) dofs: $(round(resultpercent, digits = 4))%"

    # Visualization
    if visualize
        scattersysvec!(dchi, (L/8)/maximum(abs.(U)).*U)
        update_rotation_field!(Rfield0, dchi)
        plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
            #plot_nodes(fens),
            plot_midsurface(fens, fes; x = geom0.values, facecolor = "rgb(12, 12, 123)"),
            plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
            dims = 1)
        pl = render(plots)
    end

    return true
end

function _execute_model_w_units(n = 8, visualize = true)
    # analytical solution for the vertical deflection and the midpoint of the
    # free edge
    analyt_sol=-0.3024*phun("ft");
    # Parameters:
    E=4.32e8*phun("lbf/ft^2");
    nu=0.0;
    thickness = 0.25*phun("ft"); # geometrical dimensions are in feet
    R = 25.0*phun("ft");
    L = 50.0*phun("ft");
    q = -90*phun("lbf/ft^2")
    formul = FEMMShellT3FFModule

    tolerance = R/n/1000
    fens, fes = T3block(40/360*2*pi,L/2,n,n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a), y, R*(cos(a)-1))
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    femm.transv_shear_formulation = formul.__TRANSV_SHEAR_FORMULATION_AVERAGE_B
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # rigid diaphragm
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    for i in [1,3,5]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Y
    l1 = selectnode(fens; box = Float64[-Inf Inf L/2 L/2 -Inf Inf], inflate = tolerance)
    for i in [2,4,6]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to X
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,5,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Midpoint of the free edge
    nl = selectnode(fens; box = Float64[sin(40/360*2*pi)*R sin(40/360*2*pi)*R L/2 L/2 -Inf Inf], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt[0, 0, q, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    resultpercent =   dchi.values[nl, 3][1]/analyt_sol*100
    @info "Solution for $(count(fens)*6) dofs: $(round(resultpercent, digits = 4))%"

    # Visualization
    if visualize
        scattersysvec!(dchi, (L/8)/maximum(abs.(U)).*U)
        update_rotation_field!(Rfield0, dchi)
        plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
            #plot_nodes(fens),
            plot_midsurface(fens, fes; x = geom0.values, facecolor = "rgb(12, 12, 123)"),
            plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
            dims = 1)
        pl = render(plots)
    end

    return true
end

function test_convergence(ns = [32, 64, 128])
    @info "Scordelis-Lo shell"
    for n in ns
        _execute_model(n, false)
    end
    return true
end

function test_convergence_w_units(ns = [32, 64, 128])
    @info "Scordelis-Lo shell"
    for n in ns
        _execute_model_w_units(n, false)
    end
    return true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test_convergence()
    test_convergence_w_units()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing

# ns = [2, 4, 8]
# [ Info: Solution for 54 dofs: 73.982%
# [ Info: Solution for 150 dofs: 66.6193%
# [ Info: Solution for 486 dofs: 85.5858%
# [ Info: Solution for 1734 dofs: 95.4205%


text:
include("euler/euler.jl")
using .Calculus: fibonacci_index, get_fibonacci

compute(n::Int)::Int = trunc(Int, get_fibonacci(fibonacci_index(n) + 2) / 2)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# reset_default_graph()
# using SymPy

# x, y = @vars x y
# u = x*(1-x)*y*(1-y)
# v = x*(1-x)*y*(1-y)^2

# p = x*(1-x)*y*(1-y)

# ux = diff(u,x)
# uy = diff(u,y)
# vx = diff(v,x)
# vy = diff(v,y)
# px = diff(p,x)
# py = diff(p,y)
# f = u*ux + v*uy - (diff(ux,x)+diff(uy,y)) + px
# g = u*vx + v*vy - (diff(vx,x)+diff(vy,y)) + py
# h = diff(u, x) + diff(v, y)
# println(replace(replace(sympy.julia_code(simplify(f)), ".^"=>"^"), ".*"=>"*"))
# println(replace(replace(sympy.julia_code(simplify(g)), ".^"=>"^"), ".*"=>"*"))
# println(replace(replace(sympy.julia_code(simplify(h)), ".^"=>"^"), ".*"=>"*"))

function u_exact(x,y)
    if y == 0.0
        return 1.0
    else
        return 0.0
    end
end

function v_exact(x,y)
    0.0
end

function p_exact(x,y)
    0.0
end

function ffunc(x, y)
    #-x^2*y*(x - 1)^2*(y - 1)^2*(2*y - 1) + x*y^2*(x - 1)*(2*x - 1)*(y - 1)^2 + x*y*(y - 1) - 2*x*(x - 1) + y*(x - 1)*(y - 1) - 2*y*(y - 1)
    0.0
end
function gfunc(x, y)
    #x^2*y*(x - 1)^2*(y - 1)^3*(3*y - 1) - x*y^2*(x - 1)*(2*x - 1)*(y - 1)^3 + 3*x*y*(x - 1) + 5*x*(x - 1)*(y - 1) + 2*y*(y - 1)^2
    0.0
end

function hfunc(x,y)
    # (y - 1)*(-2*x*y*(x - 1) + x*y - x*(x - 1)*(y - 1) + y*(x - 1))
    0.0
end

# nn to estimate space varying viscosity
function nu_exact(x, y)
    1 + 6 * x^2 + x / (1 + 2 * y^2)
    # 1 + x^2 + x / (1 + 2 * y^2)
    # 1 + x^2
end

using LinearAlgebra
using MAT
using AdFem
using PyPlot
using SparseArrays

m = 20
n = 20
h = 1/n
# nu = 0.01

F1 = compute_fem_source_term1(eval_f_on_gauss_pts(ffunc, m, n, h), m, n, h)
F2 = compute_fem_source_term1(eval_f_on_gauss_pts(gfunc, m, n, h), m, n, h)
H = h^2*eval_f_on_fvm_pts(hfunc, m, n, h)
B = constant(compute_interaction_matrix(m, n, h))

# compute F
xy = fem_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
u0 = @. u_exact(x,y)
v0 = @. v_exact(x,y)

xy = fvm_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
p0 = @. p_exact(x,y)
# nu_gauss = Variable(ones(4*m*n))
nu_fvm = Variable(ones(m*n))
nu_gauss = reshape(repeat(nu_fvm, 1, 4), (-1,))

Laplace = compute_fem_laplace_matrix1(nu_gauss, m, n, h)

# Laplace = constant(nu*compute_fem_laplace_matrix1(m, n, h))
function compute_residual(S)
    u, v, p = S[1:(m+1)*(n+1)], S[(m+1)*(n+1)+1:2(m+1)*(n+1)], S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    interaction = compute_interaction_term(p, m, n, h) # julia kernel needed
    f1 = compute_fem_source_term1(ugauss.*ux, m, n, h)
    f2 = compute_fem_source_term1(vgauss.*uy, m, n, h)
    f3 = -interaction[1:(m+1)*(n+1)]
    f4 = Laplace*u
    f5 = -F1
    F = f1 + f2 + f3 + f4 + f5

    g1 = compute_fem_source_term1(ugauss.*vx, m, n, h)
    g2 = compute_fem_source_term1(vgauss.*vy, m, n, h)
    g3 = -interaction[(m+1)*(n+1)+1:end]
    g4 = Laplace*v
    g5 = -F2
    G = g1 + g2 + g3 + g4 + g5

    H0 = -B * [u;v] + H

    R = [F;G;H0]
    return R
end

function compute_jacobian(S)
    u, v, p = S[1:(m+1)*(n+1)], S[(m+1)*(n+1)+1:2(m+1)*(n+1)], S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    M1 = constant(compute_fem_mass_matrix1(ux, m, n, h))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, m, n, h))

    N1 = constant(compute_fem_mass_matrix1(vy, m, n, h))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, m, n, h))

    J0 = [Fu Fv
          Gu Gv]
    J = [J0 -B'
        -B spdiag(zeros(size(B,1)))]
end

NT = 5
# S = zeros(m*n+2(m+1)*(n+1), NT+1)
# S[1:m+1, 1] = ones(m+1,1)

bd = bcnode("all", m, n, h)
# bd = [bd; bd .+ (m+1)*(n+1); ((1:m) .+ 2(m+1)*(n+1))]
bd = [bd; bd .+ (m+1)*(n+1); 2*(m+1)*(n+1)+(n-1)*m+1:2*(m+1)*(n+1)+(n-1)*m+2] # only apply Dirichlet to velocity

function solve_steady_cavityflow_one_step(S)
    residual = compute_residual(S)
    J = compute_jacobian(S)

    J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, m, n, h)
    # residual[bd] .= 0.0
    residual = scatter_update(residual, bd, zeros(length(bd)))

    d = J\residual
    S_new = S - d
    return S_new
end

function condition(i, S_arr)
    i <= NT + 1
end

function body(i, S_arr)
    S = read(S_arr, i-1)
    op = tf.print("i=",i)
    i = bind(i, op)
    S_new = solve_steady_cavityflow_one_step(S)
    S_arr = write(S_arr, i, S_new)
    return i+1, S_arr
end

# for i = 1:NT
#     residual = compute_residual(S[:,i])
#     J = compute_jacobian(S[:,i])

#     J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, m, n, h)
#     residual[bd] .= 0.0

#     d = J\residual
#     S[:,i+1] = S[:,i] - d
#     @info i, norm(residual)
# end

function plot_velocity_pressure_viscosity(k)

    figure(figsize=(14,4));
    subplot(131)
    visualize_scalar_on_fem_points(S_true[end, 1:(m+1)*(n+1)], m, n, h); title("velocity x data")
    subplot(132)
    visualize_scalar_on_fem_points(run(sess, S[end,1:(m+1)*(n+1)]), m, n, h); title("velocity x prediction")
    subplot(133)
    visualize_scalar_on_fem_points(S_true[end, 1:(m+1)*(n+1)] - run(sess, S[end,1:(m+1)*(n+1)]), m, n, h); title("velocity x error")
    savefig("figures_ones6_fvm/steady_xy_nn_velox$k.png")
    # savefig("steady_xy_ones_velox.png")

    figure(figsize=(14,4));
    subplot(131)
    visualize_scalar_on_fem_points(S_true[end, (m+1)*(n+1)+1:2*(m+1)*(n+1)], m, n, h); title("velocity y data")
    subplot(132)
    visualize_scalar_on_fem_points(run(sess, S[end,(m+1)*(n+1)+1:2*(m+1)*(n+1)]), m, n, h); title("velocity y prediction")
    subplot(133)
    visualize_scalar_on_fem_points(S_true[end, (m+1)*(n+1)+1:2*(m+1)*(n+1)] - run(sess, S[end,(m+1)*(n+1)+1:2*(m+1)*(n+1)]), m, n, h); title("velocity y error")
    savefig("figures_ones6_fvm/steady_xy_nn_veloy$k.png")
    # savefig("steady_xy_ones_veloy.png")

    figure(figsize=(14,4));
    subplot(131)
    visualize_scalar_on_fvm_points(S_true[end, 2*(m+1)*(n+1)+1:end], m, n, h); title("pressure data")
    subplot(132)
    visualize_scalar_on_fvm_points(run(sess, S[end,2*(m+1)*(n+1)+1:end]), m, n, h); title("pressure prediction")
    subplot(133)
    visualize_scalar_on_fvm_points(S_true[end, 2*(m+1)*(n+1)+1:end] - run(sess, S[end,2*(m+1)*(n+1)+1:end]), m, n, h); title("pressure difference")
    savefig("figures_ones6_fvm/steady_xy_nn_pres$k.png")
    # savefig("steady_xy_ones_pres.png")

    figure(figsize=(14,4));
    subplot(131)
    visualize_scalar_on_fvm_points(nu_exact.(x,y), m, n, h); title("viscosity exact")
    subplot(132)
    visualize_scalar_on_gauss_points(run(sess, nu_gauss), m, n, h); title("viscosity prediction");gca().invert_yaxis()
    subplot(133)
    visualize_scalar_on_fvm_points(nu_exact.(x,y).-run(sess, nu_gauss)[1:4:end], m, n, h); title("viscosity difference")
    savefig("figures_ones6_fvm/steady_xy_nn_visc$k.png")
    # savefig("steady_xy_ones_visc.png")

end

S_arr = TensorArray(NT+1)
S_arr = write(S_arr, 1, [u0; v0; p0])

i = constant(2, dtype=Int32)

_, S = while_loop(condition, body, [i, S_arr])
S = set_shape(stack(S), (NT+1, 2*(m+1)*(n+1)+m*n))

S_true = matread("steady_cavity_data.mat")["V"]
loss = mean((S[end,1:2*(m+1)*(n+1)] - S_true[end,1:2*(m+1)*(n+1)])^2)

# loss = mean((S[end,:] - S_true[end,:])^2)
loss = loss * 1e10

sess = Session(); init(sess)
@info run(sess, loss)

max_iter = 1000

# ADCME.load(sess, "nn$k.mat")

for k = 1:100
    loss_ = BFGS!(sess, loss, max_iter)
    matwrite("figures_ones6_fvm/loss$k.mat", Dict("L"=>loss_))
    close("all"); semilogy(loss_); title("loss vs. iteration")
    savefig("figures_ones6_fvm/steady_xy_loss$k.png")
    plot_velocity_pressure_viscosity(k)
    ADCME.save(sess, "figures_ones6_fvm/nn$k.mat")
end
# savefig("steady_xy_ones_loss.png")


text:
"""
	rates(param::parameters,iterations::Int64,divergence::Array,sfs::Array)

Function to solve randomly *N* scenarios. The function will create *N* models, defined by ```Analytical.parameters()```, to estimate analytically fixation and polymorphic rates for each model. The rates will be used to compute summary statistics required at ABC. The function output a HDF5 file containing the solved models, the selected DAC and the analytical rates.

If rho and/or theta are set to ```nothing```, the function will input random values given the range 0.0005:0.0005:0.01. Otherwise you can fix the values.

If gL is set to ```nothing```, the function will not account the role of the weakly selected alleles in the estimation.

# Arguments
 - `param::parameters`: mutable structure containing the model
 - `convolutedSamples::binomialDict` : structure containing the binomial convolution
 - `gH::Array{Int64,1}` : Range of strong selection coefficients
 - `gL::Union{Array{Int64,1},Nothing}`: Range of weak selection coefficients
 - `gamNeg::Array{Int64,1}` : Range of deleterious selection coefficients
 - `theta::Union{Float64,Nothing}` : Population-scaled mutation rate on coding region
 - `rho::Union{Float64,Nothing}` : Population-scaled recombination rate
 - `shape::Float64=0.184` : DFE shape parameter
 - `iterations::Int64` : Number of solutions
 - `output::String` : File to output HDF5 file
# Returns
 - `Array`: summary statistics
 - `Output`: HDF5 file containing models solved and rates.
"""
function rates(;param::parameters,
				convolutedSamples::binomialDict,
				gH::S,
				gL::S,
				gamNeg::S,
				theta::Union{Float64,Nothing}=0.001,
				rho::Union{Float64,Nothing}=0.001,
				shape::Float64=0.184,
				iterations::Int64,
				output::String,
				scheduler::String="local") where S <: Union{Array{Int64,1},UnitRange{Int64},Nothing}

	# Iterations = models to solve
	# Factor to modify input Γ(shape) parameter. Flexible Γ distribution over negative alleles
	fac     = rand(-2:0.05:2,iterations)
	afac    = @. param.al*(2^fac)

	# Deleting shape > 1. Negative alpha_x values
	idx = findall(afac .> 1)
	if !isempty(idx)
		afac[idx] = rand(afac[afac .< 1],size(idx,1))
	end

	# Random α values
	nTot    = rand(0.1:0.01:0.9,iterations)

	# Defining αW. It is possible to solve non-accounting for weak fixations
	if isnothing(gL)
		# Setting αW to 0 for all estimations
		nLow    = fill(0.0,iterations)
		# Random strong selection coefficients
		ngl     = rand(repeat([1],iterations),iterations);
	else
		# Setting αW as proportion of α
		lfac    = rand(0.0:0.05:0.9,iterations)
		nLow    = @. nTot * lfac
		# Random weak selection coefficients
		ngl     = rand(repeat(gL,iterations),iterations);
	end

	# Creating N models to iter in threads. Set N models (paramerters) and sampling probabilites (binomialDict)
	nParam  = [param for i in 1:iterations];
	nBinom  = [convolutedSamples for i in 1:iterations];

	# Random strong selection coefficients
	ngh     = rand(repeat(gH,iterations),iterations);
	# Random negative selection coefficients
	ngamNeg = rand(repeat(gamNeg,iterations),iterations);

	# Random θ on coding regions
    #=obsNeut        = sfs.p0./sum(sfs.p0)
    n              = 1 ./collect(1:(param.nn-1))
    expNeut        = n ./ sum(n)=#

	if !isnothing(theta)
		θ = fill(theta,iterations)
		#=θᵣ= fill(theta .* obsNeut ./ expNeut,iterations)=#
	else
		θ = rand(0.0005:0.0005:0.01,iterations)
		#=θᵣ= map(x -> x .* obsNeut ./ expNeut,θ)=#
	end

	# Random ρ on coding regions
	if !isnothing(rho)
		ρ = fill(rho,iterations)
	else
		ρ = rand(0.0005:0.0005:0.05,iterations)
	end

	# Estimations to distributed workers
	if scheduler == "local"
		println(scheduler)
		@time out = Distributed.pmap(iterRates,nParam, nBinom, nTot, nLow, ngh, ngl, ngamNeg, afac, θ, ρ);
		#=@time out = Distributed.pmap(iterRates,nParam, nBinom, nTot, nLow, ngh, ngl, ngamNeg, afac, θ, θᵣ, ρ);=#
	else
		println(scheduler)
		@time out = ParallelUtilities.pmapbatch( iterRates,nParam, nBinom, nTot, nLow, ngh, ngl, ngamNeg, afac, θ, ρ);
		#=@time out = ParallelUtilities.pmapbatch( iterRates,nParam, nBinom, nTot, nLow, ngh, ngl, ngamNeg, afac, θ, θᵣ, ρ);=#
	end

	# Remove the workers to free memory resources
	#=for i in workers()
		rmprocs(i)
	end=#

	# Reducing output array
	df = vcat(out...)

	# Saving models and rates
	models = DataFrame(df[:,1:8],[:B,:alLow,:alTot,:gamNeg,:gL,:gH,:al,:ρ])
	neut   = df[:,9:(8+size(param.dac,1))]
	sel    = df[:,(9+size(param.dac,1)):(8+size(param.dac,1)*2)]
	dsdn   = Array(df[:,(end-3):end])

	# Saving multiple summary statistics
	n = OrderedDict{Int,Array}()
	s = OrderedDict{Int,Array}()
	for i in eachindex(param.dac)
		n[param.dac[i]] = neut[:,i]
		s[param.dac[i]] = sel[:,i]
	end

	# Writting HDF5 file
	JLD2.jldopen(output, "a+") do file
		file[string(param.N)* "/" * string(param.n) * "/models"] = models;
		file[string(param.N)* "/" * string(param.n) * "/neut"]   = n;
		file[string(param.N)* "/" * string(param.n) * "/sel"]    = s;
		file[string(param.N)* "/" * string(param.n) * "/dsdn"]   = dsdn;
		file[string(param.N)* "/" * string(param.n) * "/dac"]    = param.dac;
	end;
end

"""
	iterRates(param::parameters,afac::Float64,bfac::Float64,alTot::Float64,alLow::Float64,divergence::Array,sfs::Array)

Estimating rates given a model for all B range.

# Arguments
 - `param::parameters`
 - `convolutedSamples::binomialDict`
 - `alTot::Float64`
 - `alLow::Float64`
 - `gH::Int64`
 - `gL::Int64`
 - `gamNeg::Int64`
 - `afac::Float64`
 - `ρ::Float64`
 - `θ::Float64`
# Output
 - `Array{Float64,2}`
"""
function iterRates(param::parameters,convolutedSamples::binomialDict,alTot::Float64,alLow::Float64,gH::Int64,gL::Int64,gamNeg::Int64,afac::Float64,θ::Float64,ρ::Float64)

	# Creating model to solve
	# Γ distribution
	param.al    = afac; param.be = abs(afac/gamNeg); param.gamNeg = gamNeg
	# α, αW
	param.alLow = alLow; param.alTot = alTot;
	# Positive selection coefficients
	param.gH    = gH;param.gL = gL
	# Mutation rate and recomb
	param.thetaMidNeutral = θ; param.rho = ρ
	#=param.thetaMidNeutral = θ; param.θᵣ .= θᵣ; param.rho = ρ=#
	# Solving θ on non-coding region and probabilites to get α value without BGS
	param.B = 0.999
	setThetaF!(param)
	setPpos!(param)

	# Allocate array to solve the model for all B values
	r = zeros(size(param.B_bins,1),(size(param.dac,1) * 2) + 12)
	for j in eachindex(param.B_bins)
		# Set B value
		param.B = param.B_bins[j]
		# Solve θ non-coding for the B value.
		setThetaF!(param)
		# Solve model for the B value
		tmp = try
			gettingRates(param,convolutedSamples.bn[param.B])
		catch e
			zeros(size(param.dac,1) *2+ 12)'
		end
		@inbounds r[j,:] = tmp
	end

	return r
end

"""
	gettingRates(gammaL,gammaH,pposL,pposH,observedData,nopos)

Estimating analytical rates of fixation and polymorphism to approach α value accouting for background selection, weakly and strong positive selection. Output values will be used to sample from a Poisson distribution the total counts of polymorphism and divergence using observed data.

# Arguments
 - `param::parameters`
 - `cnvBinom::SparseMatrixCSC{Float64,Int64}`
# Returns
 - `Array{Float64,2}` containing solved model, fixation and polymorphic rates
"""
function gettingRates(param::parameters,cnvBinom::SparseMatrixCSC{Float64,Int64})

	################################################
	# Subset rates accounting for positive alleles #
	################################################

	# Fixation
	fN       = param.B*fixNeut(param)
	fNeg     = param.B*fixNegB(param,0.5*param.pposH+0.5*param.pposL)
	fPosL    = fixPosSim(param,param.gL,0.5*param.pposL)
	fPosH    = fixPosSim(param,param.gH,0.5*param.pposH)

	ds       = fN
	dn       = fNeg + fPosL + fPosH

	# Polymorphism
	neut::Array{Float64,1} = DiscSFSNeutDown(param,cnvBinom)
	selH::Array{Float64,1} = if isinf(exp(param.gH * 2))
		DiscSFSSelPosDownArb(param,param.gH,param.pposH,cnvBinom)
	else
		DiscSFSSelPosDown(param,param.gH,param.pposH,cnvBinom)
	end
	selL::Array{Float64,1} = DiscSFSSelPosDown(param,param.gL,param.pposL,cnvBinom)
	selN::Array{Float64,1} = DiscSFSSelNegDown(param,param.pposH+param.pposL,cnvBinom)
	# Cumulative rates
	tmp = cumulativeSfs(hcat(neut,selH,selL,selN),false)
	splitColumns(matrix::Array{Float64,2}) = (view(matrix, :, i) for i in 1:size(matrix, 2));
	neut, selH, selL, selN = splitColumns(tmp)
	sel = (selH+selL)+selN

	##########
	# Output #
	##########
	analyticalValues::Array{Float64,2} = vcat(param.B,param.alLow,param.alTot,param.gamNeg,param.gL,param.gH,param.al,param.thetaMidNeutral,neut[param.dac],sel[param.dac],ds,dn,fPosL,fPosH)'

	return (analyticalValues)
end


text:
module M3GTools

  using DelimitedFiles
  using FortranFiles
  using PDBTools # from https://github.com/m3g/PDBTools

  include("./Atom.jl")
  export Atom

  include("./Bond.jl")
  export Bond

  include("./Angle.jl")
  export Angle

  include("./Dihedral.jl")
  export Dihedral

  include("./Improper.jl")
  export Improper

  include("./readpsf.jl")
  export readpsf

  include("./readprm.jl")
  export readprm!

  include("./getpdbcoords.jl")
  export getpdbcoords!

  include("./coulomb.jl")
  export coulomb

  include("./vdw.jl")
  export vdw

  include("./nhbonds.jl")
  export nhbonds

  include("./Simulation.jl")

  include("./init.jl")

  include("./select.jl")
  include("./dcdio.jl")

  include("./cm.jl")
  include("./pbc.jl")
  include("./distance.jl")
  include("./tcf.jl")
  include("./procrustes.jl")

  include("./gmdget.jl")

  include("./density.jl")
  export density
  include("./movingaverage.jl")
  export movingaverage

  include("./version.jl")

  # For gmd3D
  include("./GMD3DGrid.jl")
  include("./gmd3D_setgrid.jl")
  include("./gmd3D.jl")
  include("./gmd3D_write.jl")
  export gmd3D, gmd3D_write

  # time correlation functions
  include("./tcorr.jl")

  # compute interpolated data
  include("./sample_interpolation.jl")

end


text:
module ThreeBodyDecay

using QuadGK
using StaticArrays
using Cuba
using PartialWaveFunctions
using Parameters
using RecipesBase

export Kallen, KallenFact, Kibble,
       inphrange,
       lims,lims1,lims2,lims3
export phase
export ijk, ij_from_k
#
export cosθij
export cosθ12, cosθ23, cosθ31
#
export σkofi, σkofj
export σ1of2, σ2of3, σ3of1,
       σ1of3, σ2of1, σ3of2
include("kinematics.jl")

export x_str, minusone
include("utils.jl")

export ispositive
export TriavialWignerRotation
export wr
export cosζ
#
export cosζ12_for0, cosζ23_for0, cosζ31_for0
#
export cosζ21_for1, cosζ21_for2
export cosζ13_for1, cosζ13_for3
export cosζ32_for3, cosζ32_for2
#
export cosζ12_for3, cosζ23_for1,cosζ31_for2
#
export cosθhatk1, cosθhatk2, cosθhatk3
export cosζk1_for1, cosζk2_for2, cosζk3_for3
#
include("wignerrotations.jl")

export ThreeBodyMasses, ThreeBodySpins, ThreeBodyParities
export ThreeBodySystem
export DalitzPlotPoint
export Invariants
export randomPoint
export nt, x2, over2
export possible_helicities
export border, border31, border12, border23
export flatDalitzPlotSample
include("tbs_struct.jl")

# density
export getbinned1dDensity,
       getbinned2dDensity
export gridded_density_function
include("rand_corr.jl")

export rotz!, roty!, roty_cos!, roty_cos_inv!, boostz!
export invmasssq
export four_vectors_in_binary_decay
include("rand_gen.jl")

#lineshape
export pole, BW, BWdw
export Lineshape, BreitWigner, amp
export ScattLenApproximation, ScattLen
export Rho, iRho, ChewMandestam
export RhoQTB, iRhoQTB
#
export three_body_phase_space_integral
include("lineshape.jl")

export change_basis_3from1,
        change_basis_1from2,
        change_basis_2from3
include("cross_channel_relations.jl")

#
export jp, str2jp
export @jp_str
export ⊗
export possible_ls, possible_lsLS
export possible_coupling_schemes
include("coupling_scheme.jl")

export Zksτ
include("general_case.jl")

export NoRecoupling, ParityRecoupling, RecouplingLS
export DecayChain, DecayChainLS, DecayChainsLS
export amplitude
export itr, summed_over_polarization
include("decay_channel.jl")

include("dalitzplotsrecipe.jl")

export squaredalitz,
       squaredalitz1,
       squaredalitz2,
       squaredalitz3
export invsquaredalitz,
       invsquaredalitz1,
       invsquaredalitz2,
       invsquaredalitz3
export jacobean_squaredalitz,
       jacobean_squaredalitz1,
       jacobean_squaredalitz2,
       jacobean_squaredalitz3
include("squaredalitz.jl")

end  # module ThreeBodyDecay

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
"""
    resizekern(source_size::Tuple{Int,Int}, dest_size::Tuple{Int,Int})

Create an image resize blur kernel, for use before reducing image size to avoid aliasing.
"""
function resizekern(source_size::Tuple{Int,Int}, dest_size::Tuple{Int,Int})
    # Blur before resizing to prevent aliasing (kernel size dependent on both source and target image size)
    σ = map((o, n) -> 0.75 * o / n, source_size, dest_size)
    if first(σ) < 1
        return ImageFiltering.KernelFactors.gaussian(σ)
    else
        return ImageFiltering.KernelFactors.IIRGaussian(σ)
    end
end

"""
    sizethatfits(original_size::Tuple{Int,Int},target_shape::Tuple{Int,Int})

Takes an original image size, and fits it to a target shape for image resizing. Maintains aspect ratio.
"""
function sizethatfits(
    original_size::Tuple{Int,Int},
    target_shape::Tuple{Int,Int},
)
    if original_size[1] > original_size[2]
        target_img_size = (
            target_shape[1],
            floor(Int, target_shape[2] * (original_size[2] / original_size[1])),
        )
    else
        target_img_size = (
            floor(Int, target_shape[1] * (original_size[1] / original_size[2])),
            target_shape[2],
        )
    end
    return target_img_size
end

"""
    loadResizePadImageToFit(img_path::String, sets::Settings)

Load image and reesize it to fit inside the target image shap, while maintaining
aspect ratio and preventing aliasing.
"""
loadResizePadImageToFit(img_path::String, sets::Settings) =
    loadResizePadImageToFit(FileIO.load(img_path), sets)
function loadResizePadImageToFit(img::T, sets::Settings) where {T<:AbstractArray}
    img_size = size(img)
    target_img_size = sizethatfits(img_size, sets.image_shape)
    kern = resizekern(img_size, target_img_size)
    return resizePadImageToFit(img, target_img_size, sets, kern)
end

"""
    resizePadImageToFit(img_path::String, sets::Settings, kern::Tuple{ImageFiltering.KernelFactors.ReshapedOneD,ImageFiltering.KernelFactors.ReshapedOneD})

Loads and prepares (resizes + pads) an image to fit within a given shape.
Returns the image and the padding.
"""
function resizePadImageToFit(
    img::Array{T},
    target_img_size::Tuple{Int,Int},
    settings::Settings,
    kern::Tuple{
        ImageFiltering.KernelFactors.ReshapedOneD,
        ImageFiltering.KernelFactors.ReshapedOneD,
    },
) where {T<:ColorTypes.Color}

    if size(img,1) > target_img_size[1] #Apply blur first if reducing size, to avoid aliasing
        imgr = ImageTransformations.imresize(
            ImageFiltering.imfilter(img, kern, NA()),
            target_img_size,
        )
    else
        imgr = ImageTransformations.imresize(img, target_img_size)
    end

    # Determine top and left padding
    vpad_top = floor(Int, (settings.image_shape[1] - target_img_size[1]) / 2)
    hpad_left = floor(Int, (settings.image_shape[2] - target_img_size[2]) / 2)

    # Determine bottom and right padding accounting for rounding of top and left (to ensure accuate result image size if source has odd dimensions)
    vpad_bottom = settings.image_shape[1] - (vpad_top + target_img_size[1])
    hpad_right = settings.image_shape[2] - (hpad_left + target_img_size[2])

    padding = [hpad_left, vpad_top, hpad_right, vpad_bottom]

    # Pad image
    return padarray(
            imgr,
            Fill(
                zero(eltype(img)),
                (vpad_top, hpad_left),
                (vpad_bottom, hpad_right),
            ),
        ),
        padding
end

"""
    load(ds::LabelledImageDataset, settings::Settings; limitfirst::Int = -1)

Load images from a populated `LabelledImageDataset` into memory.
"""
function load(
    ds::LabelledImageDataset,
    settings::Settings;
    indexes::Union{Array{Int}} = [],
)
    if length(indexes) == 1
        numimages = 1
    elseif length(indexes) > 0 && length(indexes) < length(ds.image_paths)
        numimages = length(indexes)
        @info "Loading $(numimages) images from $(ds.name) dataset into memory"
    else
        numimages = length(ds.image_paths)
        indexes = 1:numimages
        @info "Loading all ($(numimages)) images from $(ds.name) dataset into memory"
    end
    kern = resizekern(ds.image_size_lims, settings.image_shape)

    firstimg, padding = loadResizePadImageToFit(
        ds.image_paths[1],
        settings,
    )
    imgsize = size(firstimg)
    lds = LoadedDataset(
            imstack_mat = xtype(undef,imgsize[1],imgsize[2],settings.image_channels,numimages),
            paddings = Vector{Vector{Int}}(undef, 0),
            labels = ds.labels[indexes],
        )

    j = 1
    ProgressMeter.@showprogress 2 "Loading images..." for i in indexes
        img = FileIO.load(ds.image_paths[i])
        img_size = size(img)
        target_img_size = sizethatfits(img_size, settings.image_shape)
        img_resized, padding = resizePadImageToFit(
            img,
            target_img_size,
            settings,
            kern,
        )
        lds.imstack_mat[:, :, :, j] = collect(permutedims(
            channelview(img_resized)[1:settings.image_channels, :, :],
            [2, 3, 1],
        ))
        push!(lds.paddings, padding)
        j += 1
    end
    return lds
end