text:
#TODO: convert two function calls into Union{COOTen,ThirdOrderSymTensor}
#=------------------------------------------------------------------------------
              Routines for searching over alpha/beta parameters
------------------------------------------------------------------------------=#
"""------------------------------------------------------------------------------
  This function aligns graphs using their tensor representations. These routines
  call the param_search functions for the associated method used.

  Inputs
  ------
  * A, B - (ThirdOrderSymTensor):
	Two third order tensors representing the presence of triangles within the
	network. A must be larger than B, else the routines will be called with the
	parameters swapped.
  * method - (String):
	The choice of method used to align the methods. Options include 'LambdaTAME',
	'LowRankTAME', and 'TAME'.

  Outputs
  -------
	* 'best_TAME_PP_tris' - (Int):
		The largest number of triangles matched over all iterations.
	* 'max_triangle_match' - (Int):
	  The maximum number of triangles matchable. This is the minimum between the
	  number of triangles in graphs A and B.
	* 'best_TAME_PP_x' - (Array{Float,2}):
	  The best iterate found over all the alphas and betas specified by the user.
	  when 'LambdaTAME' and 'LowRankTAME' are called, this is replaced by the U
	  and V components of the best iterate.

-----------------------------------------------------------------------------"""
function align_tensors(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor;
					   method::String="LambdaTAME",no_matching=false,kwargs...)

	#put larger tensor on the left
	if B.n > A.n
		results = align_tensors(B,A;method = method, no_matching=no_matching,kwargs...)
		#flip the matchings if A and B were swapped
		if method == "LambdaTAME" ||  method == "LowRankTAME"
			best_TAME_PP_tris, max_triangle_match, U_best, V_best, best_matching = results
			return best_TAME_PP_tris, max_triangle_match, U_best, V_best, Dict((j,i) for (i,j) in best_matching)
		elseif method == "TAME"
			best_TAME_PP_tris, max_triangle_match, best_TAME_PP_x, best_matching = results
			return best_TAME_PP_tris, max_triangle_match, best_TAME_PP_x, Dict((j,i) for (i,j) in best_matching)
		end
	end

	if method == "LambdaTAME"
		return ΛTAME_param_search(A,B;kwargs...)
	elseif method == "LowRankTAME"
		return LowRankTAME_param_search(A,B;no_matching = no_matching,kwargs...)
	elseif method == "TAME"
		return TAME_param_search(A,B;no_matching = no_matching,kwargs...)
	else
		throw(ArgumentError("method must be one of 'LambdaTAME', 'LowRankTAME',or 'TAME'."))
	end

end

"""------------------------------------------------------------------------------
  This function aligns graphs using their tensor representations and returns the
  profiled version of the algorithms. These routines call the 'param_search'
  functions for the associated method used.

  Inputs
  ------
  * A, B - (ThirdOrderSymTensor):
	Two third order tensors representing the presence of triangles within the
	network. A must be larger than B, else the routines will be called with the
	parameters swapped.
  * method - (String):
	The choice of method used to align the methods. Options include 'LambdaTAME',
	'LowRankTAME', and 'TAME'.
  * 'no_matching' - (Bool):
	Will not run the matching routines when True. This is useful when studying
	the ranks of the iterates. Any counts which may rely on the matchings are
	replaced by -1.
  Outputs
  -------
  * 'best_TAME_PP_tris' - (Int):
    The largest number of triangles matched over all iterations.
  * 'max_triangle_match' - (Int):
	The maximum number of triangles matchable. This is the minimum between the
	number of triangles in graphs A and B.
  * 'best_TAME_PP_x' - (Array{Float,2}):
	The best iterate found over all the alphas and betas specified by the user.
	when 'LambdaTAME' and 'LowRankTAME' are called, this is replaced by the U
	and V components of the best iterate.
  * 'best_matching' - (Dict{Int,Int}):
	The best matching found, maps from A to B. When the method flips A and B,
	the dictionary is also flipped when returned.
  * profile - (Dict):
	A dictionary storing the profiling results from each of the methods. Please
	see the '_profiled' versions of the code to see what is returned by each
	function.
------------------------------------------------------------------------------"""
function align_tensors_profiled(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor;
					            method::String="LambdaTAME",no_matching=false,kwargs...)

	#put larger tensor on the left
	if B.n > A.n
		results =  align_tensors_profiled(B,A;method = method, no_matching=no_matching,kwargs...)
		#flip the matchings if A and B were swapped
		if method == "LambdaTAME" ||  method == "LowRankTAME"
			best_TAME_PP_tris, max_triangle_match, U_best, V_best, best_matching,profile = results
			return best_TAME_PP_tris, max_triangle_match, U_best, V_best, Dict((j,i) for (i,j) in best_matching), profile
		elseif method == "TAME"
			best_TAME_PP_tris, max_triangle_match, best_TAME_PP_x, best_matching,profile = results
			return best_TAME_PP_tris, max_triangle_match, best_TAME_PP_x, Dict((j,i) for (i,j) in best_matching), profile
		end

	end

	if method == "LambdaTAME"
		return ΛTAME_param_search_profiled(A,B;kwargs...)
	elseif method == "LowRankTAME"
		return LowRankTAME_param_search_profiled(A,B;no_matching = no_matching,kwargs...)
	elseif method == "TAME"
		return TAME_param_search_profiled(A,B;no_matching = no_matching,kwargs...)
	else
		throw(ArgumentError("method must be one of 'LambdaTAME', 'LowRankTAME', or 'TAME'."))
	end
end

function ΛTAME_param_search(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor;
                            iter::Int = 15,tol::Float64=1e-6,
							alphas::Array{F,1}=[.5,1.0],
							betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001]) where {F <: AbstractFloat}

    max_triangle_match = min(size(A.indices,1),size(B.indices,1))
    total_triangles = size(A.indices,1) + size(B.indices,1)

    U = Array{Float64,2}(undef,A.n,iter)
    V = Array{Float64,2}(undef,B.n,iter)

    best_TAME_PP_tris = -1
    best_i  = -1
	best_j = -1
	best_matching = Dict{Int,Int}()

    for α in alphas
        for beta in betas

			U,V = ΛTAME(A,B,beta,iter,tol,α)
			search_tris, i, j, matching = search_Krylov_space(A,B,U,V)
			println("α:$(α) -- β:$(beta) finished -- tri_match:$search_tris -- max_tris $(max_triangle_match) -- best tri_match:$best_TAME_PP_tris")
            if search_tris > best_TAME_PP_tris
                best_TAME_PP_tris = search_tris
                best_i = i
				best_j = j
				best_matching = matching
            end
        end
    end
	println("best i:$best_i -- best j:$best_j")
	return best_TAME_PP_tris, max_triangle_match, U[best_i,:], V[best_j,:], best_matching
end

function ΛTAME_param_search_profiled(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor;
	                                 iter::Int = 15,tol::Float64=1e-6, alphas::Array{F,1}=[.5,1.0],
							         betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001]) where {F <: AbstractFloat}

    max_triangle_match = min(size(A.indices,1),size(B.indices,1))
    best_TAME_PP_tris = -1
    best_i  = -1
	best_j = -1
	best_matching = Dict{Int,Int}()

	m = A.n
	n = B.n

	results = Dict(
		"TAME_timings" => Array{Float64,1}(undef,length(alphas)*length(betas)),
		"Krylov Timings"=> Array{Float64,1}(undef,length(alphas)*length(betas))
	)
	exp_index = 1

    U = Array{Float64,2}(undef,m,iter)
    V = Array{Float64,2}(undef,n,iter)

    for α in alphas
        for beta in betas

			((U,V),runtime) = @timed ΛTAME(A,B,beta,iter,tol,α)
			results["TAME_timings"][exp_index] = runtime

			#search the Krylov Subspace
			((search_tris, i, j, matching),runtime) = @timed search_Krylov_space(A,B,U,V)
			results["Krylov Timings"][exp_index] = runtime
			exp_index += 1

			if search_tris > best_TAME_PP_tris
				best_matching = matching
                best_TAME_PP_tris = search_tris
                best_i = i
                best_j = j
            end

			println("α:$(α) -- β:$(beta) finished -- tri_match:$search_tris -- max_tris $(max_triangle_match) -- best tri_match: $best_TAME_PP_tris")
        end
    end

	println("best i:$best_i -- best j:$best_j")
	return best_TAME_PP_tris, max_triangle_match, U[:,best_i], V[:,best_j], best_matching, results

end

#add in SparseSymmetricTensors.jl function definitions
function TAME_param_search(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor;
                           iter::Int = 15,tol::Float64=1e-6, alphas::Array{F,1}=[.5,1.0],
						   betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
						   kwargs...) where {F <: AbstractFloat}

    max_triangle_match = min(size(A.indices,1),size(B.indices,1))
    total_triangles = size(A.indices,1) + size(B.indices,1)
	best_TAME_PP_tris::Int = -1
	best_matching = Dict{Int,Int}()

	m = A.n
	n = B.n

	best_TAME_PP_x = Array{Float64,2}(undef,m,n)

    for α in alphas
        for β in betas

			x, triangle_count, matching = TAME(A,B,β,iter,tol,α;W = ones(A.n,B.n),kwargs...)

			if triangle_count > best_TAME_PP_tris
				best_matching = matching
                best_TAME_PP_tris = triangle_count
				best_TAME_PP_x = copy(x)
            end
			println("α:$(α) -- β:$β finished -- tri_match:$(triangle_count) -- max_tris $(max_triangle_match) -- best tri_match:$(best_TAME_PP_tris)")
        end

    end

	return best_TAME_PP_tris, max_triangle_match, best_TAME_PP_x, best_matching
end

function TAME_param_search_profiled(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor;
                                    iter::Int = 15,tol::Float64=1e-6, alphas::Array{F,1}=[.5,1.0],
						            betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
						            profile::Bool=false,profile_aggregation="all",
						            kwargs...) where {F <: AbstractFloat}
    max_triangle_match = min(size(A.indices,1),size(B.indices,1))
    total_triangles = size(A.indices,1) + size(B.indices,1)
	best_TAME_PP_tris = -1
	best_matching = Dict{Int,Int}()

	m = A.n
	n = B.n

	best_TAME_PP_x = Array{Float64,2}(undef,m,n)
	experiment_profiles = Array{Tuple{String,Dict{String,Union{Array{Float64,1},Array{Array{Float64,1},1}}}},1}(undef,0)

    for α in alphas
        for β in betas

			x, triangle_count, matching, experiment_profile = TAME_profiled(A,B,β,iter,tol,α;W = ones(m,n),kwargs...)
			push!(experiment_profiles,("α:$(α)_β:$(β)",experiment_profile))

			if triangle_count > best_TAME_PP_tris
				best_matching = matching
                best_TAME_PP_tris = triangle_count
				best_TAME_PP_x = copy(x)
            end
			println("α:$(α) -- β:$β finished -- tri_match:$(best_TAME_PP_tris) -- max_tris $(max_triangle_match)")
        end

    end

	return best_TAME_PP_tris, max_triangle_match, best_TAME_PP_x, best_matching, experiment_profiles

end

function LowRankTAME_param_search(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor;
						          iter::Int = 15,tol::Float64=1e-6,
						          U_0::Array{Float64,2} = ones(A.n,1),
						          V_0::Array{Float64,2} = ones(B.n,1),
						          alphas::Array{F,1}=[.5,1.0],
						          betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
						          kwargs...) where {F <: AbstractFloat}
    max_triangle_match = min(size(A.indices,1),size(B.indices,1))
    total_triangles = size(A.indices,1) + size(B.indices,1)
    best_TAME_PP_tris = -1
	best_matching = Dict{Int,Int}()

	m = A.n
	n = B.n

	best_TAME_PP_U = ones(m,1)
	best_TAME_PP_V = ones(n,1)

    for α in alphas
        for β in betas

			U, V, triangle_count,matching =
				LowRankTAME(A,B,U_0,V_0,β,iter,tol,α;kwargs...)

            if triangle_count > best_TAME_PP_tris
                best_TAME_PP_tris = triangle_count
				best_matching = matching
				best_TAME_PP_U = copy(U)
				best_TAME_PP_V = copy(V)

            end
			println("α:$(α) -- β:$(β) -- tri_match:$(triangle_count) -- max_tris:$(max_triangle_match) -- best tri match:$best_TAME_PP_tris")
        end

    end

	return best_TAME_PP_tris, max_triangle_match, best_TAME_PP_U, best_TAME_PP_V, best_matching

end

function LowRankTAME_param_search_profiled(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor;
                           iter::Int = 15,tol::Float64=1e-6,
						   alphas::Array{F,1}=[.5,1.0],
						   betas::Array{F,1} =[1000.0,100.0,10.0,1.0,0.0,0.1,0.01,0.001],
						   kwargs...) where {F <: AbstractFloat}
    max_triangle_match = min(size(A.indices,1),size(B.indices,1))
    total_triangles = size(A.indices,1) + size(B.indices,1)
	best_TAME_PP_tris = -1
	best_matching = Dict{Int,Int}()

	best_TAME_PP_U = ones(A.n,1)
	best_TAME_PP_V = ones(B.n,1)

	experiment_profiles = []

    for α in alphas
        for β in betas

			U, V, triangle_count,matching, experiment_profile =
				LowRankTAME_profiled(A,B,ones(A.n,1),ones(B.n,1), β,iter,tol,α;kwargs...)

			push!(experiment_profiles,("α:$(α)_β:$(β)",experiment_profile))

            if triangle_count > best_TAME_PP_tris
				best_TAME_PP_tris = triangle_count
				best_matching = matching
				best_TAME_PP_U = copy(U)
				best_TAME_PP_V = copy(V)

            end
			println("α:$(α) -- β:$(β) -- tri_match:$(best_TAME_PP_tris) -- max_tris $(max_triangle_match)")
        end

    end

	return best_TAME_PP_tris, max_triangle_match, best_TAME_PP_U, best_TAME_PP_V,best_matching, experiment_profiles

end
#=------------------------------------------------------------------------------
             		    Spectral Relaxation Routines
------------------------------------------------------------------------------=#

"""------------------------------------------------------------------------------
  The LambdaTAME method. Method starts with a uniform prior and computes the
  power method for each of the tensors passed it. These iterates are stored and
  the best rank one matching is picked out of all the quadratic pairs of vectors
  computed. This portion of the algorith computes the contractions, the full
  procedure can be found in the associated '_param_search' algorithms.

  Inputs
  ------
  * A,B - (ThirdOrderSymTensor):
	The tensors to align against one another.
  * β - (Float):
	The shift to use on the iterations.
  * α -(float):
	The mixing parameter for combining in the starting iterates.
  * 'max_iter' - (Int):
	The maximum number of iterations to run.
  * tol - (Float):
	The tolerence to solve the algorithm to. Computes the tolerance by measuring
	the absolute value of the difference between the computed eigenvalues.
  * 'update_user' - (Int):
	Specifies how frequently output messages should be printed to the user.
	Default is -1 which means no output messages are printed, else if the value
	is k, then the kth iterate will print out an update.
  Output
  ------
  * U, V - (Array{Float,2})
    Returns the computed contractions started with a uniform prior.
------------------------------------------------------------------------------"""
function ΛTAME(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor, β::Float64,
               max_iter::Int,tol::Float64,α::Float64;update_user=-1)

    U = zeros(A.n,max_iter+1)
    V = zeros(B.n,max_iter+1) #store initial in first column

    U[:,1] = ones(A.n)
    U[:,1] /=norm(U[:,1])

    V[:,1] = ones(B.n)
    V[:,1] /=norm(U[:,1])

    sqrt_β = β^(.5)

    lambda = Inf
    i = 1

    while true

        U[:,i+1] = tensor_vector_contraction(A,U[:,i])
        V[:,i+1] = tensor_vector_contraction(B,V[:,i])

        lambda_A = (U[:,i+1]'*U[:,i])
        lambda_B = (V[:,i+1]'*V[:,i])
        new_lambda = lambda_A*lambda_B

        if β != 0.0
            U[:,i+1] .+= sqrt_β*U[:,i+1]
            V[:,i+1] .+= sqrt_β*V[:,i+1]
        end

        if α != 1.0
            U[:,i+1] = α*U[:,i+1] + (1 -α)*U[:,1]
            V[:,i+1] = α*V[:,i+1] + (1 -α)*V[:,1]
        end

        U[:,i+1] ./= norm(U[:,i+1])
        V[:,i+1] ./= norm(V[:,i+1])

		if update_user != -1 && i % update_user == 0
            println("iteration $(i)    λ_A: $(lambda_A) -- λ_B: $(lambda_B) -- newλ: $(new_lambda)")
		end

        if abs(new_lambda - lambda) < tol || i >= max_iter
            return U[:,1:i], V[:,1:i]
        else
            lambda = new_lambda
            i += 1
        end

    end

end

#runs TAME, but reduces down to lowest rank form first
function LowRankTAME(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,W::Array{F,2},
                       β::F, max_iter::Int,tol::F,α::F;
					   max_rank::Int = minimum((A.n,B.n)),kwargs...) where {F <:AbstractFloat}

    dimension = minimum((A.n,B.n))

	(U_k,S,VT),t = @timed svd(W)
	singular_indexes = [i for i in 1:minimum((max_rank,length(C_S))) if S[i] > S[1]*eps(Float64)*dimension]

	U = U_k[:,singular_indexes]
	V = VT[:,singular_indexes]*diagm(S[singular_indexes])

	return LowRankTAME(A,B,U,V,β,max_iter,tol,α;kwargs...)
end

"""------------------------------------------------------------------------------
  The low rank implementation of TAME, computes the terms use the mixed property
  generalized to the tensor case.

  Inputs
  ------
  * A,B - (ThirdOrderSymTensor):
	The tensors to align against one another.
  * 'U_0', 'V_0' - (Array{Float,2}):
	The low rank components of the starting iteration X = UV'. Iterates are
	normalized before the iterations begin.
  * β - (Float):
	The shift to use on the iterations.
  * α -(float):
	The mixing parameter for combining in the starting iterates.
  * 'max_iter' - (Int):
	The maximum number of iterations to run.
  * tol - (Float):
	The tolerence to solve the algorithm to. Computes the tolerance by measuring
	the absolute value of the difference between the computed eigenvalues.
  * 'max_rank' - (Int):
	Specify the maximum rank of each of the iterates. Default makes it so that
	only singular values small enough to be considered zero are truncated.
  * 'update_user' - (Int):
	Specifies how frequently output messages should be printed to the user.
	Default is -1 which means no output messages are printed, else if the value
	is k, then the kth iterate will print out an update.
  * 'no_matching' - (Bool):
	Specifies whether or not to run the matching and scoring portions of the
	algorithm. Useful if only the iterates are desired.
  * 'low_rank_matching' - (Bool):
	Specifies whether or not to run the low rank matching procedure from [1].
	This is useful when speed is needed, but may lead to regressions in the
	matching performance.

  Output
  ------
  * 'best_U', 'best_V' - (Array{Float,2})
    Returns the components to the iteration which matched the most triangles.
  * 'best_triangle_count' - (Int)
    The maximum number of triangles matched.
  * 'best_matching' - (Dict{Int,Int})
	The matching computed between the two graphs, maps the vertices of A to the
	vertices of B.

  Citation
  --------
  [1] - H. Nassar, N. Veldt, S. Mohammadi, A. Grama, and D. F. Gleich,
		“Low rank spectral network alignment,” in Proceedings of the 2018
		World Wide Web Conference, 2018, pp. 619–628
------------------------------------------------------------------------------"""
function LowRankTAME(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,
                          U_0::Array{F,2},V_0::Array{F,2}, β::F, max_iter::Int,tol::F,α::F;
						  max_rank::Int = minimum((A.n,B.n)),update_user::Int=-1,
						  no_matching=false,low_rank_matching=false) where {F <:AbstractFloat}

	@assert size(U_0,2) == size(V_0,2)

	dimension = minimum((A.n,B.n))

	best_triangle_count::Int = -1
	best_matching = Dict{Int,Int}()
    best_x = zeros(size(U_0,1),size(V_0,1))
    best_index = -1
	triangles = -1
	gaped_triangles = -1

	normalization_factor = sqrt(tr((V_0'*V_0)*(U_0'*U_0)))
	U_0 ./= sqrt(normalization_factor)
	V_0 ./= sqrt(normalization_factor)

	U_k = copy(U_0)
	V_k = copy(V_0)

	best_U::Array{F,2} = copy(U_k)
	best_V::Array{F,2} = copy(U_k)

    lambda = Inf

    for i in 1:max_iter

		A_comps, B_comps = get_kron_contract_comps(A,B,U_k,V_k)
		lam = tr((B_comps'*V_k)*(U_k'*A_comps))

		if α != 1.0 && β != 0.0
			U_temp = hcat(sqrt(α) * A_comps, sqrt(α * β) * U_k, sqrt(1-α) * U_0)
			V_temp = hcat(sqrt(α) * B_comps, sqrt(α * β) * V_k, sqrt(1-α) * V_0)
		elseif α != 1.0
			U_temp = hcat(sqrt(α)*A_comps, sqrt(1-α)*U_0)
			V_temp = hcat(sqrt(α)*B_comps, sqrt(1-α)*V_0)
		elseif β != 0.0
			U_temp = hcat(A_comps, sqrt(β) * U_k)
			V_temp = hcat(B_comps, sqrt(β) * V_k)
		else
			U_temp = A_comps
			V_temp = B_comps
		end

		A_Q,A_R = qr(U_temp)
		B_Q,B_R = qr(V_temp)

		core = A_R*B_R'
		C_U,C_S,C_Vt = svd(core)
		singular_indexes= [i for i in 1:1:minimum((max_rank,length(C_S))) if C_S[i] > C_S[1]*eps(Float64)*dimension]

		U_k_1 = A_Q*C_U[:,singular_indexes]
		V_k_1 = B_Q*(C_Vt[:,singular_indexes]*diagm(C_S[singular_indexes]))

		normalization_factor = sqrt(tr((V_k_1'*V_k_1)*(U_k_1'*U_k_1)))

		U_k_1 ./= sqrt(normalization_factor)
		V_k_1 ./= sqrt(normalization_factor)

		#TODO: need to redo or remove this
		Y, Z = get_kron_contract_comps(A,B,U_k_1,V_k_1)

		lam = tr((Y'*U_k_1)*(V_k_1'*Z))

		if !no_matching
			#evaluate the matchings
			if low_rank_matching
				triangles, gaped_triangles,matching = TAME_score(A,B,U_k_1,V_k_1)
			else
				triangles, gaped_triangles,matching = TAME_score(A,B,U_k_1*V_k_1')
			end

			if triangles > best_triangle_count
				best_matching = matching
				best_triangle_count  = triangles
				best_U = copy(U_k_1)
				best_V = copy(V_k_1)
			end
		end

		if update_user != -1 && i % update_user == 0
			println("λ_$i: $(lam) -- rank:$(length(singular_indexes)) -- tris:$(triangles) -- gaped_t:$(gaped_triangles)")
		end

        if abs(lam - lambda) < tol
    		return best_U, best_V, best_triangle_count, best_matching
        end

		#get the low rank factorization for the next one
		U_k = copy(U_k_1)
		V_k = copy(V_k_1)

		lambda = lam

    end

	return best_U, best_V, best_triangle_count, best_matching
end

"""------------------------------------------------------------------------------
  The low rank implementation of TAME, computes the terms use the mixed property
  generalized to the tensor case.

  Inputs
  ------
  * A,B - (ThirdOrderSymTensor):
	The tensors to align against one another.
  * 'U_0', 'V_0' - (Array{Float,2}):
	The low rank components of the starting iteration X = UV'. Iterates are
	normalized before the iterations begin.
  * β - (Float):
	The shift to use on the iterations.
  * α -(float):
	The mixing parameter for combining in the starting iterates.
  * 'max_iter' - (Int):
	The maximum number of iterations to run.
  * tol - (Float):
	The tolerence to solve the algorithm to. Computes the tolerance by measuring
	the absolute value of the difference between the computed eigenvalues.
  * 'max_rank' - (Int):
	Specify the maximum rank of each of the iterates. Default makes it so that
	only singular values small enough to be considered zero are truncated.
  * 'update_user' - (Int):
	Specifies how frequently output messages should be printed to the user.
	Default is -1 which means no output messages are printed, else if the value
	is k, then the kth iterate will print out an update.
  * 'no_matching' - (Bool):
	Specifies whether or not to run the matching and scoring portions of the
	algorithm. Useful if only the iterates are desired.
  * 'low_rank_matching' - (Bool):
	Specifies whether or not to run the low rank matching procedure from [1].
	This is useful when speed is needed, but may lead to regressions in the
	matching performance.

  Output
  ------
	* 'best_U', 'best_V' - (Array{Float,2})
	  Returns the components to the iteration which matched the most triangles.
	* 'best_triangle_count' - (Int)
	  The maximum number of triangles matched.
	* 'best_matching' - (Dict{Int,Int})
	  The matching computed between the two graphs, maps the vertices of A to the
	  vertices of B.
	* 'experiment_profile' - (Dict{String,Union{Array{F,1},Array{Array{F,1},1}}}):
	  The experiment profile computed, keys for the experiment data collected are
	  as follows.
	  +'ranks' - The ranks of each iterate X_k.
	  +'contraction_timings' - Time taken to compute each contraction.
	  +'svd_timings' - Time taken to compute the svd.
	  +'qr_timings' - Time taken to compute the QR factorations.
	  +'matched_tris' - The number of triangles matched by each iterate.
	  +'sing_vals' - The singular values of each iterate X_k.
	  +'matching_times' - Time taken to solve the matchings.
	  +'soring_timings' - Time taken to score each matching.

  Citation
  --------
	[1] - H. Nassar, N. Veldt, S. Mohammadi, A. Grama, and D. F. Gleich,
		  “Low rank spectral network alignment,” in Proceedings of the 2018
		  World Wide Web Conference, 2018, pp. 619–628
------------------------------------------------------------------------------"""
function LowRankTAME_profiled(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,
                          U_0::Array{F,2},V_0::Array{F,2}, β::F, max_iter::Int,tol::F,α::F;
						  max_rank::Int = minimum((A.n,B.n)),update_user::Int=-1,
						  no_matching::Bool=false,low_rank_matching::Bool=false) where {F <:AbstractFloat}

	@assert size(U_0,2) == size(V_0,2)

	dimension = minimum((A.n,B.n))

	best_triangle_count::Int = -1
	best_matching = Dict{Int,Int}()
    best_x = zeros(size(U_0,1),size(V_0,1))
    best_index = -1
	triangles = -1
	gaped_triangles = -1

	normalization_factor = sqrt(tr((V_0'*V_0)*(U_0'*U_0)))
	U_0 ./= sqrt(normalization_factor)
	V_0 ./= sqrt(normalization_factor)

	experiment_profile = Dict{String,Union{Array{F,1},Array{Array{F,1},1}}}(
		"ranks"=>Array{Float64,1}(undef,0),
		"contraction_timings"=>Array{Float64,1}(undef,0),
		"svd_timings"=>Array{Float64,1}(undef,0),
		"qr_timings"=>Array{Float64,1}(undef,0),
		"matched_tris"=>Array{Float64,1}(undef,0),
		"sing_vals"=>Array{Array{Float64,1},1}(undef,0),
		"matching_timings"=>Array{Float64,1}(undef,0),
		"scoring_timings"=>Array{Float64,1}(undef,0)
	)

	U_k = copy(U_0)
	V_k = copy(V_0)

	best_U::Array{F,2} = copy(U_k)
	best_V::Array{F,2} = copy(U_k)

    lambda = Inf

    for i in 1:max_iter

		(A_comps, B_comps),t = @timed get_kron_contract_comps(A,B,U_k,V_k)
		push!(experiment_profile["contraction_timings"],t)

		lambda_k_1 = tr((B_comps'*V_k)*(U_k'*A_comps))

		if α != 1.0 && β != 0.0
			U_temp = hcat(sqrt(α) * A_comps, sqrt(α * β) * U_k, sqrt(1-α) * U_0)
			V_temp = hcat(sqrt(α) * B_comps, sqrt(α * β) * V_k, sqrt(1-α) * V_0)
		elseif α != 1.0
			U_temp = hcat(sqrt(α)*A_comps, sqrt(1-α)*U_0)
			V_temp = hcat(sqrt(α)*B_comps, sqrt(1-α)*V_0)
		elseif β != 0.0
			U_temp = hcat(A_comps, sqrt(β) * U_k)
			V_temp = hcat(B_comps, sqrt(β) * V_k)
		else
			U_temp = A_comps
			V_temp = B_comps
		end

		(A_Q,A_R),t_A = @timed qr(U_temp)
		(B_Q,B_R),t_B = @timed qr(V_temp)
		push!(experiment_profile["qr_timings"],t_A + t_B)

		core = A_R*B_R'
		(C_U,C_S::Array{Float64,1},C_Vt),t = @timed svd(core)
		push!(experiment_profile["svd_timings"],t)

		singular_indexes= [i for i in 1:1:minimum((max_rank,length(C_S))) if C_S[i] > C_S[1]*eps(Float64)*dimension]
		push!(experiment_profile["sing_vals"],C_S)
		push!(experiment_profile["ranks"],float(length(singular_indexes)))

		U_k_1 = A_Q*C_U[:,singular_indexes]
		V_k_1 = B_Q*(C_Vt[:,singular_indexes]*diagm(C_S[singular_indexes]))

		normalization_factor = sqrt(tr((V_k_1'*V_k_1)*(U_k_1'*U_k_1)))

		U_k_1 ./= sqrt(normalization_factor)
		V_k_1 ./= sqrt(normalization_factor)

		if !no_matching
			#evaluate the matchings
			if low_rank_matching
				triangles, gaped_tris, matching, matching_time, scoring_time = TAME_score(A,B,U_k_1,V_k_1;return_timings=true)
			else
				triangles, gaped_tris, matching, matching_time, scoring_time = TAME_score(A,B,U_k_1*V_k_1';return_timings=true)
			end

			push!(experiment_profile["matched_tris"],float(triangles))
			push!(experiment_profile["matching_timings"],float(matching_time))
			push!(experiment_profile["scoring_timings"], float(scoring_time))

			if triangles > best_triangle_count
				best_matching = matching
				best_triangle_count  = triangles
				best_U = copy(U_k_1)
				best_V = copy(V_k_1)
			end
		end

		if update_user != -1 && i % update_user == 0
			println("λ_$i: $(lambda_k_1) -- rank:$(length(singular_indexes)) -- tris:$(triangles) -- gaped_t:$(gaped_tris)")
		end

		if abs(lambda_k_1 - lambda) < tol || i >= max_iter
			#=
			triangles,_= TAME_score(A,B,sparse(best_U*best_V');return_timings=false)
			if triangles > best_triangle_count
				best_triangle_count = triangles
			end
			=#
			return best_U, best_V, best_triangle_count,best_matching, experiment_profile
		end
		#get the low rank factorization for the next one
		U_k = copy(U_k_1)
		V_k = copy(V_k_1)

		lambda = lambda_k_1
	end

	return best_U, best_V, best_triangle_count,best_matching, experiment_profile
end

function setup_tame_data(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor)
	return _index_triangles_nodesym(A.n,A.indices), _index_triangles_nodesym(B.n,B.indices)
end

function _index_triangles_nodesym(n,Tris::Array{Int,2})

	Ti = [ Vector{Tuple{Int,Int}}(undef, 0) for i in 1:n ]
	for (ti,tj,tk) in eachrow(Tris)
		push!(Ti[ti], (tj,tk))
		push!(Ti[tj], (ti,tk))
		push!(Ti[tk], (ti,tj))
	end
	sort!.(Ti)
	return Ti
end

"""------------------------------------------------------------------------------
  An implementation of TAME, for our experiments we use the original C++ code,
  but we include these routines for any additional experimentation desired.

  Inputs
  ------
  * A,B - (ThirdOrderSymTensor):
	The tensors to align against one another.
  * W - (Array{Float,2}):
	The starting iteration, Iterate is normalized before the iterations begin.
  * β - (Float):
	The shift to use on the iterations.
  * α -(float):
	The mixing parameter for combining in the starting iterates.
  * 'max_iter' - (Int):
	The maximum number of iterations to run.
  * tol - (Float):
	The tolerence to solve the algorithm to. Computes the tolerance by measuring
	the absolute value of the difference between the computed eigenvalues.
  * 'update_user' - (Int):
	Specifies how frequently output messages should be printed to the user.
	Default is -1 which means no output messages are printed, else if the value
	is k, then the kth iterate will print out an update.
  * 'no_matching' - (Bool):
	Specifies whether or not to run the matching and scoring portions of the
	algorithm. Useful if only the iterates are desired.

  Output
  ------
	* 'best_x'- (Array{Float,2})
	  Returns the components to the iteration which matched the most triangles.
	  Reshapes the iterate x into a matrix.
	* 'best_triangle_count' - (Int)
	  The maximum number of triangles matched.
	* 'best_matching' - (Dict{Int,Int})
	  The matching computed between the two graphs, maps the vertices of A to the
	  vertices of B.
------------------------------------------------------------------------------"""
function TAME(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor,β::F, max_iter::Int,
			  tol::F,α::F;update_user::Int=-1,W::Array{F,2}=ones(A.n,B.n),
			  no_matching=false,) where {F <:AbstractFloat}

    dimension = minimum((A.n,B.n))

	A_Ti, B_Ti = setup_tame_data(A,B)

	best_x = Array{Float64,2}(undef,A.n,B.n)
    best_triangle_count = -1
	best_index = -1
	best_matching = Dict{Int,Int}()

    x0 = reshape(W,A.n*B.n)
	x0 ./=norm(x0)
    x_k = copy(x0)

    i = 1
    lambda = Inf

    while true

	    x_k_1 = impTTVnodesym(A.n, B.n, x_k, A_Ti, B_Ti)
		#x_k_1 = implicit_contraction(A,B,x_k)

		#println("norm diff is:",norm(x_k_1_test - x_k_1)/norm(x_k_1_test))
        new_lambda = dot(x_k_1,x_k)

        if β != 0.0
            x_k_1 .+= β * x_k
        end

        if α != 1.0
            x_k_1 = α * x_k_1 + (1 - α) * x0
        end

		x_k_1 ./= norm(x_k_1)

		if !no_matching

			X = reshape(x_k_1,A.n,B.n)
			triangles, gaped_triangles, matching =  TAME_score(A,B,X)

			if update_user != -1 && i % update_user == 0
				println("finished iterate $(i):tris:$triangles -- gaped_t:$gaped_triangles")
			end

			if triangles > best_triangle_count
				best_matching = matching
				best_x = copy(x_k_1)
				best_triangle_count = triangles
				best_iterate = i
			end

		end

		if update_user != -1 && i % update_user == 0
			println("λ: $(new_lambda)")
		end

        if abs(new_lambda - lambda) < tol || i >= max_iter
  			return reshape(best_x,A.n,B.n), best_triangle_count, best_matching
        else
            x_k = copy(x_k_1)
            lambda = new_lambda
            i += 1
        end

    end

end

"""------------------------------------------------------------------------------
  The low rank implementation of TAME, computes the terms use the mixed property
  generalized to the tensor case.

  Inputs
  ------
  * A,B - (ThirdOrderSymTensor):
	The tensors to align against one another.
  * W - (Array{Float,2}):
	The starting iteration, Iterate is normalized before the iterations begin.
  * β - (Float):
	The shift to use on the iterations.
  * α -(float):
	The mixing parameter for combining in the starting iterates.
  * 'max_iter' - (Int):
	The maximum number of iterations to run.
  * tol - (Float):
	The tolerence to solve the algorithm to. Computes the tolerance by measuring
	the absolute value of the difference between the computed eigenvalues.
  * 'update_user' - (Int):
	Specifies how frequently output messages should be printed to the user.
	Default is -1 which means no output messages are printed, else if the value
	is k, then the kth iterate will print out an update.
  * 'no_matching' - (Bool):
	Specifies whether or not to run the matching and scoring portions of the
	algorithm. Useful if only the iterates or profiling is desired.

  Output
  ------
    * 'best_x'- (Array{Float,2})
      Returns the components to the iteration which matched the most triangles.
      Reshapes the iterate x into a matrix.
	* 'best_triangle_count' - (Int)
	  The maximum number of triangles matched.
	* 'best_matching' - (Dict{Int,Int})
	  The matching computed between the two graphs, maps the vertices of A to the
	  vertices of B.
	* 'experiment_profile' - (Dict{String,Union{Array{F,1},Array{Array{F,1}}}}):
	  The experiment profile computed, keys for the experiment data collected are
	  as follows.
	  +'contraction_timings' - Time taken to compute each contraction.
	  +'matched_tris' - The number of triangles matched by each iterate.
	  +'gaped triangles' - The number of unmatched triangles, equal to the
	                       maximum - matched triangles.
	  +'sing_vals' - The singular values of each iterate X_k.
	  +'ranks' - The ranks of each iterate X_k.
	  +'matching_timings' - Time taken to solve the matchings.
	  +'soring_timings' - Time taken to score each matching.

------------------------------------------------------------------------------"""
function TAME_profiled(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor, β::F, max_iter::Int,
	                   tol::F,α::F;update_user::Int=-1,W::Array{F,2} = ones(m,n),
					   no_matching::Bool=false) where {F <:AbstractFloat}

    dimension = minimum((A.n,B.n))

    experiment_profile = Dict{String,Union{Array{F,1},Array{Array{F,1}}}}(
		"contraction_timings"=>Array{F,1}(undef,0),
		"matching_timings"=>Array{F,1}(undef,0),
		"scoring_timings"=>Array{F,1}(undef,0),
		"matched triangles"=>Array{F,1}(undef,0),
		"gaped triangles"=>Array{F,1}(undef,0),
		"sing_vals"=>Array{Array{F,1},1}(undef,0),
		"ranks"=>Array{F,1}(undef,0)
	)

	A_Ti, B_Ti = setup_tame_data(A,B)

	best_x = Array{Float64,2}(undef,A.n,B.n)
    best_triangle_count::Int = -1
	best_index = -1
	best_matching = Dict{Int,Int}()

    x0 = reshape(W,A.n*B.n)
	x0 ./=norm(x0)
    x_k = copy(x0)

    i = 1
    lambda = Inf

    while true

		x_k_1,t = @timed impTTVnodesym(A.n, B.n, x_k, A_Ti, B_Ti)
		push!(experiment_profile["contraction_timings"],t)

        new_lambda = dot(x_k_1,x_k)

        if β != 0.0
            x_k_1 .+= β * x_k
        end

        if α != 1.0
            x_k_1 = α * x_k_1 + (1 - α) * x0
        end

        x_k_1 ./= norm(x_k_1)

		S = svdvals(reshape(x_k_1,(A.n,B.n)))

		rank = 0.0
		for i in 1:length(S)
			if S[i] > S[1]*eps(Float64)*dimension
				rank = rank + 1
			end
		end

		push!(experiment_profile["sing_vals"],S)
		push!(experiment_profile["ranks"],rank)

		if !no_matching

			triangles, gaped_triangles, matching, matching_time, scoring_time = TAME_score(A,B,reshape(x_k_1,A.n,B.n);return_timings=true)

			push!(experiment_profile["matching_timings"],matching_time)
			push!(experiment_profile["scoring_timings"],scoring_time)
			push!(experiment_profile["matched triangles"],float(triangles))
			push!(experiment_profile["gaped triangles"],float(gaped_triangles))

			if update_user != -1 && i % update_user == 0
				println("finished iterate $(i):tris:$(triangles) -- gaped_t:$(gaped_triangles)")
			end

			if triangles > best_triangle_count
				best_matching = matching
				best_x = copy(x_k_1)
				best_triangle_count = triangles
				best_iterate = i
			end
		end

		if update_user != -1 && i % update_user == 0
			println("λ_$i: $(new_lambda)")
		end

        if abs(new_lambda - lambda) < tol || i >= max_iter
   			return reshape(best_x,A.n,B.n), best_triangle_count, best_matching, experiment_profile
        else
            x_k = copy(x_k_1)
            lambda = new_lambda
            i += 1
        end

    end

end


text:
# 8bit Posits
one(::Type{T}) where {T<:PositAll8} = T(0x40)
minusone(::Type{T}) where {T<:PositAll8} = T(0xc0)
zero(::Type{T}) where {T<:PositAll8} = T(0x00)
floatmax(::Type{T}) where {T<:PositAll8} = T(0x7f)
floatmin(::Type{T}) where {T<:PositAll8} = T(0x01)

# 16bit Posits
one(::Type{T}) where {T<:PositAll16} = T(0x4000)
minusone(::Type{T}) where {T<:PositAll16} = T(0xc000)
zero(::Type{T}) where {T<:PositAll16} = T(0x0000)
floatmax(::Type{T}) where {T<:PositAll16} = T(0x7fff)
floatmin(::Type{T}) where {T<:PositAll16} = T(0x0001)

# 24bit Posits
one(::Type{T}) where {T<:PositAll24} = T(0x4000_0000)
minusone(::Type{T}) where {T<:PositAll24} = T(0xc000_0000)
zero(::Type{T}) where {T<:PositAll24} = T(0x0000_0000)
floatmax(::Type{T}) where {T<:PositAll24} = T(0x7fff_ff00)
floatmin(::Type{T}) where {T<:PositAll24} = T(0x0000_0100)

# 32bit Posits
one(::Type{Posit32}) = Posit32(0x4000_0000)
minusone(::Type{Posit32}) = Posit32(0xc000_0000)
zero(::Type{Posit32}) = Posit32(0x0000_0000)
floatmax(::Type{Posit32}) = Posit32(0x7fff_ffff)
floatmin(::Type{Posit32}) = Posit32(0x0000_0001)

# -(x::T) where {T<:PositAll8} = x*minusone(T)
# -(x::T) where {T<:PositAll16} = x*minusone(T)
# -(x::T) where {T<:PositAll24} = x*minusone(T)
# -(x::Posit32) = x*minusone(Posit32)

function -(x::T) where {T<:PositAll8}
    if x == zero(T) || x == notareal(T) # don't change sign for 0 and NaR
        return x
    else    # subtracting from 0x00 (two's complement def for neg)
        return T(0x00 - UInt8(x))
    end
end

function -(x::T) where {T<:PositAll16}
    if x == zero(T) || x == notareal(T) # don't change sign for 0 and NaR
        return x
    else    # subtracting from 0x0000 (two's complement def for neg)
        return T(0x0000 - UInt16(x))
    end
end

function -(x::T) where {T<:Union{Posit32,PositAll24}}
    if x == zero(T) || x == notareal(T) # don't change sign for 0 and NaR
        return x
    else    # subtracting from 0x0000 (two's complement def for neg)
        return T(0x0000_0000 - UInt32(x))
    end
end

# generalize also for objects of the type AbstractPosit
minusone(p::AbstractPosit) = minusone(typeof(p))

notareal(::Type{T}) where {T<:PositAll8} = T(0x80)
notareal(::Type{T}) where {T<:PositAll16} = T(0x8000)
notareal(::Type{T}) where {T<:PositAll24} = T(0x8000_0000)
notareal(::Type{Posit32}) = Posit32(0x8000_0000)
notareal(p::AbstractPosit) = notareal(typeof(p))

signbit(p::Posit8) = signbit(reinterpret(Int8,p))
signbit(p::Posit16) = signbit(reinterpret(Int16,p))
signbit(p::Posit32) = signbit(reinterpret(Int32,p))
signbit(p::T) where {T<:Union{PositX1,PositX2}} = signbit(reinterpret(Int32,p))

isfinite(p::T) where {T<:AbstractPosit} = p != notareal(T)

function sign(p::T) where {T <: AbstractPosit}
    if signbit(p)       # negative and infinity case
        if isfinite(p)  # negative
            return minusone(T)
        else            # infinity
            return zero(T)
        end
    else                # positive and zero case
        if iszero(p)    # zero
            return zero(T)
        else            # positive
            return one(T)
        end
    end
end

ext: jl


text:
import Libdl: dlext

DIR = @__DIR__

if DIR != nothing
    const DEP_DIR = joinpath(DIR, "..", "deps", "usr")
else
    const DEP_DIR = Pkg.dir("RDKafka", "deps", "usr")
end

function find_librdkafka(dep_dir)
    base = "librdkafka"
    for fname in [
        "$base.$dlext",
        "$base.$dlext.1",
        "$base.1.$dlext",
    ]
        path = abspath(joinpath(dep_dir, fname))
        if isfile(path)
            return path
        end
    end
    return nothing
end

const LIBRDKAFKA = find_librdkafka(DEP_DIR)


text:
module Utils
using jInv.Mesh
using Distributed
using SparseArrays
using LinearAlgebra
using Statistics
using Random
using DelimitedFiles

include("SourceParallelism.jl");
include("modelUtils.jl");
include("SrcRcvUtils.jl");
include("BoundModel.jl");
include("readModelAndGenerateMeshMref.jl");
include("readModelAndGenerateNodalMeshMref.jl");
include("getHelmholtzOperator.jl")
include("remoteChangePmis.jl")

end


text:
#Is this usable data?
function isconcrete{T<:Number}(x::T)
    !isna(x) && isfinite(x)
end

function isconcrete(x::MathConst)
    return true
end

function isconcrete(x)
    !isna(x)
end

function isallconcrete(xs)
    ans = true
    for x in xs
        if !isconcrete(x)
            ans = false
            break
        end
    end
    return ans
end

function concretize(xss::AbstractVector...)
    if all(map(isallconcrete, xss))
        return xss
    end

    count = 0
    for j in 1:length(xss[1])
        for xs in xss
            if !isconcrete(xs[j])
                @goto next_j1
            end
        end

        count += 1
        @label next_j1
    end

    yss = Array(AbstractVector, length(xss))
    for (i, xs) in enumerate(xss)
        yss[i] = Array(eltype(xs), count)
    end

    k = 1
    for j in 1:length(xss[1])
        for xs in xss
            if !isconcrete(xs[j])
                @goto next_j2
            end
        end

        for (i, xs) in enumerate(xss)
            yss[i][k] = xs[j]
        end
        k += 1

        @label next_j2
    end

    return tuple(yss...)
end

# How many concrete elements in an iterable
function concrete_length(xs)
    n = 0
    for x in xs
        if isconcrete(x)
            n += 1
        end
    end
    n
end

function concrete_length{T}(xs::DataArray{T})
    n = 0
    for i = 1:length(xs)
        if !xs.na[i] && isconcrete(xs.data[i]::T)
            n += 1
        end
    end
    n
end

function concrete_length(xs::Iterators.Chain)
    n = 0
    for obj in xs.xss
        n += concrete_length(obj)
    end
    n
end

function concrete_minimum(xs)
    if done(xs, start(xs))
        error("argument must not be empty")
    end

    x_min = first(xs)
    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x)
            x_min = x
            break
        end
    end

    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x) && x < x_min
            x_min = x
        end
    end
    return x_min
end

function concrete_maximum(xs)
    if done(xs, start(xs))
        error("argument must not be empty")
    end

    x_max = first(xs)
    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x)
            x_max = x
            break
        end
    end

    for x in xs
        if Gadfly.isconcrete(x) && isfinite(x) && x > x_max
            x_max = x
        end
    end
    return x_max
end

function concrete_minmax{T<:Real}(xs, xmin::T, xmax::T)
    if eltype(xs) <: Base.Callable
        return xmin, xmax
    end

    for x in xs
        if isconcrete(x)
            xT = convert(T, x)
            if isnan(xmin) || xT < xmin
                xmin = xT
            end
            if isnan(xmax) || xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax{T<:Real, TA}(xs::DataArray{TA}, xmin::T, xmax::T)
    for i = 1:length(xs)
        if !xs.na[i]
            x = xs.data[i]::TA
            xT = convert(T, x)
            if isnan(xmin) || xT < xmin
                xmin = xT
            end
            if isnan(xmax) || xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax{T}(xs, xmin::T, xmax::T)
    for x in xs
        if isconcrete(x)
            xT = convert(T, x)
            if xT < xmin
                xmin = xT
            end
            if xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax{T, TA}(xs::DataArray{TA}, xmin::T, xmax::T)
    for i = 1:length(xs)
        if !xs.na[i]
            x = xs.data[i]::TA
            xT = convert(T, x)
            if xT < xmin
                xmin = xT
            end
            if xT > xmax
                xmax = xT
            end
        end
    end
    xmin, xmax
end

function concrete_minmax{T<:Real}(xs::Iterators.Chain, xmin::T, xmax::T)
    for obj in xs.xss
        xmin, xmax = concrete_minmax(obj, xmin, xmax)
    end
    xmin, xmax
end

function nonzero_length(xs)
    n = 0
    for x in xs
        if x != 0
            n += 1
        end
    end
    n
end

# Create a new object of type T from a with missing values (i.e., those set to
# nothing) inherited from b.
function inherit{T}(a::T, b::T)
    c = copy(a)
    inherit!(c, b)
    c
end

function inherit!{T}(a::T, b::T)
    for field in fieldnames(T)
        aval = getfield(a, field)
        bval = getfield(b, field)
        # TODO: this is a hack to let non-default labelers overide the defaults
        if aval === nothing || aval === string || aval == showoff
            setfield!(a, field, bval)
        elseif typeof(aval) <: Dict && typeof(bval) <: Dict
            merge!(aval, getfield(b, field))
        end
    end
    nothing
end

isnothing(u) = is(u, nothing)
issomething(u) = !isnothing(u)

negate(f) = x -> !f(x)

function has{T,N}(xs::AbstractArray{T,N}, y::T)
    for x in xs
        if x == y
            return true
        end
    end
    return false
end

Maybe(T) = Union(T, Nothing)

function lerp(x::Float64, a::Float64, b::Float64)
    a + (b - a) * max(min(x, 1.0), 0.0)
end

# Generate a unique id, primarily for assigning IDs to SVG elements.
let next_unique_svg_id_num = 0
    global unique_svg_id
    function unique_svg_id()
        uid = @sprintf("id%d", next_unique_svg_id_num)
        next_unique_svg_id_num += 1
        uid
    end
end

# Remove any markup or whitespace from a string.
function escape_id(s::String)
    s = replace(s, r"<[^>]*>", "")
    s = replace(s, r"\s", "_")
    s
end

# Can a numerical value be treated as an integer
is_int_compatable(::Integer) = true
is_int_compatable{T <: FloatingPoint}(x::T) = abs(x) < maxintfloat(T) && float(int(x)) == x
is_int_compatable(::Any) = false

## Return a DataFrame with x, y column suitable for plotting a function.
#
# Args:
#  f: Function/Expression to be evaluated.
#  a: Lower bound.
#  b: Upper bound.
#  n: Number of points to evaluate the function at.
#
# Returns:
#  A data frame with "x" and "f(x)" columns.
#
function evalfunc(f::Function, a, b, n)
    @assert n > 1

    step = (b - a) / (n - 1)
    xs = Array(typeof(a + step), n)
    for i in 1:n
        xs[i] = a + (i-1) * step
    end

    df = DataFrame(xs, map(f, xs))
    # NOTE: 'colnames!' is the older deprecated name. 'names!' was also defined
    # but threw an error.
    try
        names!(df, [:x, :f_x])
    catch
        colnames!(df, ["x", "f_x"])
    end
    df
end

# Construct a jscall to store arbitrary data in the element
function jsdata(key::String, value::String, arg::Vector{Measure}=Measure[])
    return jscall(
        """
        data("$(escape_string(key))", $(value))
        """, arg)
end

# Construct jscall properties to store arbitrary data in plotroot elements.
function jsplotdata(key::String, value::String, arg::Vector{Measure}=Measure[])
    return jscall(
        """
        plotroot().data("$(escape_string(key))", $(value))
        """, arg)
end

function svg_color_class_from_label(label::String)
    return @sprintf("color_%s", escape_id(label))
end

# TODO: Delete when 0.3 compatibility is dropped
if VERSION < v"0.4-dev"
    using Dates

    function Showoff.showoff{T <: Union(Date, DateTime)}(ds::AbstractArray{T}, style=:none)
        years = Set()
        months = Set()
        days = Set()
        hours = Set()
        minutes = Set()
        seconds = Set()
        for d in ds
            push!(years, Dates.year(d))
            push!(months, Dates.month(d))
            push!(days, Dates.day(d))
            push!(hours, Dates.hour(d))
            push!(minutes, Dates.minute(d))
            push!(seconds, Dates.second(d))
        end
        all_same_year         = length(years)   == 1
        all_one_month         = length(months)  == 1 && 1 in months
        all_one_day           = length(days)    == 1 && 1 in days
        all_zero_hour         = length(hours)   == 1 && 0 in hours
        all_zero_minute       = length(minutes) == 1 && 0 in minutes
        all_zero_seconds      = length(minutes) == 1 && 0 in minutes
        all_zero_milliseconds = length(minutes) == 1 && 0 in minutes

        # first label format
        label_months = false
        label_days = false
        f1 = "u d, yyyy"
        f2 = ""
        if !all_zero_seconds
            f2 = "HH:MM:SS.sss"
        elseif !all_zero_seconds
            f2 = "HH:MM:SS"
        elseif !all_zero_hour || !all_zero_minute
            f2 = "HH:MM"
        else
            if !all_one_day
                first_label_format = "u d yyyy"
            elseif !all_one_month
                first_label_format = "u yyyy"
            elseif !all_one_day
                first_label_format = "yyyy"
            end
        end
        if f2 != ""
            first_label_format = string(f1, " ", f2)
        else
            first_label_format = f1
        end

        labels = Array(String, length(ds))
        labels[1] = Dates.format(ds[1], first_label_format)
        d_last = ds[1]
        for (i, d) in enumerate(ds[2:end])
            if Dates.year(d) != Dates.year(d_last)
                if all_one_day && all_one_month
                    f1 = "yyyy"
                elseif all_one_day && !all_one_month
                    f1 = "u yyyy"
                else
                    f1 = "u d, yyyy"
                end
            elseif Dates.month(d) != Dates.month(d_last)
                f1 = all_one_day ? "u" : "u d"
            elseif Dates.day(d) != Dates.day(d_last)
                f1 = "d"
            else
                f1 = ""
            end

            if f2 != ""
                f = string(f1, " ", f2)
            elseif f1 != ""
                f = f1
            else
                f = first_label_format
            end

            labels[i+1] = Dates.format(d, f)
            d_last = d
        end

        return labels
    end
else
    using Base.Dates
end

# TODO: This is a clusterfuck. I should really just wrap Date types to force
# them to behave how I want.

if !method_exists(/, (Dates.Day, Dates.Day))
    /(a::Dates.Day, b::Dates.Day) = a.value / b.value
end

if VERSION < v"0.4.0-dev"
    Base.convert(::Type{Dates.Millisecond}, d::Dates.Day) =
        Dates.Millisecond(24 * 60 * 60 * 1000 * Dates.value(d))
end

if !method_exists(/, (Dates.Day, Real))
    /(a::Dates.Day, b::Real) = Dates.Day(round(Integer, (a.value / b)))
end
/(a::Dates.Day, b::FloatingPoint) = convert(Dates.Millisecond, a) / b

if !method_exists(/, (Dates.Millisecond, Dates.Millisecond))
    /(a::Dates.Millisecond, b::Dates.Millisecond) = a.value / b.value
end

if !method_exists(/, (Dates.Millisecond, Real))
    /(a::Dates.Millisecond, b::Real) = Dates.Millisecond(round(Integer, (a.value / b)))
end
/(a::Dates.Millisecond, b::FloatingPoint) = Dates.Millisecond(round(Integer, (a.value / b)))

if !method_exists(-, (Dates.Date, Dates.DateTime))
    -(a::Dates.Date, b::Dates.DateTime) = convert(Dates.DateTime, a) - b
end

+(a::Dates.Date, b::Dates.Millisecond) = convert(Dates.DateTime, a) + b

if !method_exists(-, (Dates.DateTime, Dates.Date))
    -(a::Dates.DateTime, b::Dates.Date) = a - convert(Dates.DateTime, b)
end

if !method_exists(/, (Dates.Day, Dates.Millisecond))
    /(a::Dates.Day, b::Dates.Millisecond) = convert(Dates.Millisecond, a) / b
end

for T in [Dates.Hour, Dates.Minute, Dates.Second, Dates.Millisecond]
    if !method_exists(-, (Dates.Date, T))
        @eval begin
            -(a::Dates.Date, b::$(T)) = convert(Dates.DateTime, a) - b
        end
    end
end

#if !method_exists(*, (FloatingPoint, Dates.Day))
    *(a::FloatingPoint, b::Dates.Day) = Dates.Day(round(Integer, (a * b.value)))
    *(a::Dates.Day, b::FloatingPoint) = b * a
    *(a::FloatingPoint, b::Dates.Millisecond) = Dates.Millisecond(round(Integer, (a * b.value)))
    *(a::Dates.Millisecond, b::FloatingPoint) = b * a
#end

# Arbitrarily order colors
function color_isless(a::ColorValue, b::ColorValue)
    return color_isless(convert(RGB{Float32}, a), convert(RGB{Float32}, b))
end

function color_isless(a::AlphaColorValue, b::AlphaColorValue)
    return color_isless(convert(RGBA{Float32}, a), convert(RGBA{Float32}, b))
end

function color_isless(a::RGB{Float32}, b::RGB{Float32})
    if a.r < b.r
        return true
    elseif a.r == b.r
        if a.g < b.g
            return true
        elseif a.g == b.g
            return a.b < b.b
        else
            return false
        end
    else
        return false
    end
end

function color_isless(a::RGBA{Float32}, b::RGBA{Float32})
    if a.c < b.c
        return true
    elseif a.c == b.c
        return a.alpha < b.alpha
    else
        return false
    end
end

function group_color_isless{S, T <: ColorValue}(a::(@compat Tuple{S, T}),
                                                b::(@compat Tuple{S, T}))
    if a[1] < b[1]
        return true
    elseif a[1] == b[1]
        return color_isless(a[2], b[2])
    else
        return false
    end
end


text:
function remove_node!(integrator::DiffEqBase.DEIntegrator, I...)
    idxs = getindices(integrator.u, I...)
    for c in full_cache(integrator)
        remove_node!(c, I...)
    end
    if DiffEqBase.is_diagonal_noise(integrator.sol.prob)
      for c in DiffEqBase.ratenoise_cache(integrator)
          remove_node!(c, I...)
      end
    end
    remove_node_non_user_cache!(integrator, I) # required to do noise correctly
end

function add_node!(integrator::DiffEqBase.DEIntegrator, x, I...)
    cur_len = length(integrator.u)
    add_len = length(x)
    last_idx = length(integrator.u[I...].nodes)
    idx_start = getindices(integrator.u, last_idx)[end] + 1
    idxs = idx_start:idx_start+add_len-1
    for c in full_cache(integrator)
        add_node!(c, fill!(similar(x, eltype(c)),0), I...)
    end
    if DiffEqBase.is_diagonal_noise(integrator.sol.prob)
      for c in DiffEqBase.ratenoise_cache(integrator)
          add_node!(c, fill!(similar(x, eltype(c)),0), I...)
      end
    end
    add_node_non_user_cache!(integrator, idxs, x, I...) # required to do noise correctly
end

function add_node!(integrator::DiffEqBase.DEIntegrator, x)
    cur_len = length(integrator.u)
    add_len = length(x)
    last_idx = length(integrator.u.nodes)
    idx_start = getindices(integrator.u, last_idx)[end] + 1
    idxs = idx_start:idx_start+add_len-1
    for c in full_cache(integrator)
        add_node!(c, fill!(similar(x, eltype(c)),0))
    end
    if DiffEqBase.is_diagonal_noise(integrator.sol.prob)
      for c in DiffEqBase.ratenoise_cache(integrator)
          add_node!(c, fill!(similar(x, eltype(c)),0))
      end
    end
    add_node_non_user_cache!(integrator, idxs, fill!(similar(x, eltype(x)),0)) # required to do noise correctly
end

reshape(m::AbstractMultiScaleArray, i::Int...) = m

function remove_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,node)
  i = length(integrator.u)
  resize_non_user_cache!(integrator,integrator.cache,i)
end
function remove_node_non_user_cache!(integrator::DiffEqBase.AbstractSDEIntegrator,node)
  if DiffEqBase.is_diagonal_noise(integrator.sol.prob)
    remove_node_noise!(integrator,node)
    for c in rand_cache(integrator)
      remove_node!(c,node...)
    end
  end
end

function remove_node_noise!(integrator,node)
  for c in integrator.W.S₁
    remove_node!(c[2],node...)
    if DiffEqBase.alg_needs_extra_process(integrator.alg)
      remove_node!(c[3],node...)
    end
  end
  for c in integrator.W.S₂
    remove_node!(c[2],node...)
    if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
      remove_node!(c[3],node...)
    end
  end
  remove_node!(integrator.W.dW,node...)
  remove_node!(integrator.W.dWtilde,node...)
  remove_node!(integrator.W.dWtmp,node...)
  remove_node!(integrator.W.curW,node...)

  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    remove_node!(integrator.W.curZ,node...)
    remove_node!(integrator.W.dZtmp,node...)
    remove_node!(integrator.W.dZtilde,node...)
    remove_node!(integrator.W.dZ,node...)
  end
end

function add_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,idxs,x)
  i = length(integrator.u)
  resize_non_user_cache!(integrator,integrator.cache,i)
end
function add_node_non_user_cache!(integrator::DiffEqBase.AbstractODEIntegrator,idxs,x,node...)
  i = length(integrator.u)
  resize_non_user_cache!(integrator,integrator.cache,i)
end
function add_node_non_user_cache!(integrator::DiffEqBase.AbstractSDEIntegrator,idxs,x,node...)
  if DiffEqBase.is_diagonal_noise(integrator.sol.prob)
    add_node_noise!(integrator,idxs,x,node...)
    for c in rand_cache(integrator)
      add_node!(c,copy(x),node...)
    end
  end
end
function add_node_non_user_cache!(integrator::DiffEqBase.AbstractSDEIntegrator,idxs,x)
  if DiffEqBase.is_diagonal_noise(integrator.sol.prob)
    add_node_noise!(integrator,idxs,x)
    for c in rand_cache(integrator)
      add_node!(c,copy(x))
    end
  end
end

function add_node_noise!(integrator,idxs,x,node...)
  for c in integrator.W.S₁
    add_node!(c[2],copy(x),node...)
    if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
      add_node!(c[3],copy(x),node...)
    end
    StochasticDiffEq.fill_new_noise_caches!(integrator,c,c[1],idxs)
  end
  for c in integrator.W.S₂
    add_node!(c[2],copy(x),node...)
    if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
      add_node!(c[3],copy(x),node...)
    end
    StochasticDiffEq.fill_new_noise_caches!(integrator,c,c[1],idxs)
  end

  add_node!(integrator.W.dW,copy(x),node...)
  integrator.W.dW[idxs] .= zero(eltype(integrator.u))
  add_node!(integrator.W.curW,copy(x),node...)
  integrator.W.curW[idxs] .= zero(eltype(integrator.u))
  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    add_node!(integrator.W.dZ,copy(x),node...)
    integrator.W.dZ[idxs] .= zero(eltype(integrator.u))
    add_node!(integrator.W.curZ,copy(x),node...)
    integrator.W.curZ[idxs] .= zero(eltype(integrator.u))
  end

  i = length(integrator.u)
  add_node!(integrator.W.dWtilde,copy(x),node...)
  add_node!(integrator.W.dWtmp,copy(x),node...)
  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    add_node!(integrator.W.dZtmp,copy(x),node...)
    add_node!(integrator.W.dZtilde,copy(x),node...)
  end

  # fill in rands
  fill!(@view(integrator.W.curW[idxs]),zero(eltype(integrator.u)))
  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    fill!(@view(integrator.W.curZ[idxs]),zero(eltype(integrator.u)))
  end
end

function add_node_noise!(integrator,idxs,x)
  for c in integrator.W.S₁
    add_node!(c[2],copy(x))
    if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
      add_node!(c[3],copy(x))
    end
    StochasticDiffEq.fill_new_noise_caches!(integrator,c,c[1],idxs)
  end
  for c in integrator.W.S₂
    add_node!(c[2],copy(x))
    if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
      add_node!(c[3],copy(x))
    end
    StochasticDiffEq.fill_new_noise_caches!(integrator,c,c[1],idxs)
  end

  add_node!(integrator.W.dW,copy(x))
  integrator.W.dW[idxs] .= zero(eltype(integrator.u))
  add_node!(integrator.W.curW,copy(x))
  integrator.W.curW[idxs] .= zero(eltype(integrator.u))
  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    add_node!(integrator.W.dZ,copy(x))
    integrator.W.dZ[idxs] .= zero(eltype(integrator.u))
    add_node!(integrator.W.curZ,copy(x))
    integrator.W.curZ[idxs] .= zero(eltype(integrator.u))
  end

  i = length(integrator.u)
  add_node!(integrator.W.dWtilde,copy(x))
  add_node!(integrator.W.dWtmp,copy(x))
  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    add_node!(integrator.W.dZtmp,copy(x))
    add_node!(integrator.W.dZtilde,copy(x))
  end

  # fill in rands
  fill!(@view(integrator.W.curW[idxs]),zero(eltype(integrator.u)))
  if StochasticDiffEq.alg_needs_extra_process(integrator.alg)
    fill!(@view(integrator.W.curZ[idxs]),zero(eltype(integrator.u)))
  end
end


text:
mutable struct DLNode{T}
    value::T
    pred::Union{DLNode{T}, Nothing}
    succ::Union{DLNode{T}, Nothing}
    DLNode(v) = new{typeof(v)}(v, nothing, nothing)
end

function insertpost(prevnode, node)
    succ = prevnode.succ
    prevnode.succ = node
    node.pred = prevnode
    node.succ = succ
    if succ != nothing
        succ.pred =  node
    end
    node
end

function printfromroot(root)
    print(root.value)
    while root.succ != nothing
        root = root.succ
        print(" -> $(root.value)")
    end
    println()
end

node1 = DLNode(1)
printfromroot(node1)
node2 = DLNode(2)
node3 = DLNode(3)
insertpost(node1, node2)
printfromroot(node1)
insertpost(node1, node3)
printfromroot(node1)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module PharmaceuticalClassification

import LightGraphs
import MetaGraphs
import ProgressMeter

export PharmClass
export PharmGraph
export all_neighbors
export available_systems
export build_graph
export build_graph!
export children
export equivalent
export inneighbors
export normalize
export outneighbors
export parents
export system_matches
export vertices

const all_neighbors = MetaGraphs.all_neighbors
const inneighbors = MetaGraphs.inneighbors
const outneighbors = MetaGraphs.outneighbors

include("types.jl")

include("atc_fix.jl")
include("build.jl")
include("countlines.jl")
include("edge.jl")
include("iterate.jl")
include("normalize.jl")
include("system.jl")
include("traversal.jl")
include("vertex.jl")

end


text:
using Chordal
using Documenter
using Literate

# For Plots.jl
# https://discourse.julialang.org/t/plotting-errors-when-building-documentation-using-plots-jl-and-documenter-jl/67849
ENV["GKSwstype"]="100"

EXCLUDED_EXAMPLES = ["chordal_graph.jl", "sdp_standard.jl"]

# utility function from https://github.com/JuliaOpt/Convex.jl/blob/master/docs/make.jl
fix_math_md(content) = replace(content, r"\$\$(.*?)\$\$"s => s"```math\1```")

# utility functions from https://github.com/oxfordcontrol/COSMO.jl/blob/master/docs/make.jl
fix_suffix(filename) = replace(filename, ".jl" => ".md")
function postprocess(content)
      """
      The source files for all examples can be found in [/examples](https://github.com/tjdiamandis/Chordal.jl/tree/main/examples).
      """ * content
end

examples_path = joinpath(@__DIR__, "../examples/")
examples = filter(x -> endswith(x, ".jl") && !in(x, EXCLUDED_EXAMPLES), readdir(examples_path))
build_path =  joinpath(@__DIR__, "src", "examples/")

for example in examples
      Literate.markdown(
        examples_path * example, build_path;
        preprocess = fix_math_md,
        postprocess = postprocess,
        flavor = Literate.DocumenterFlavor(),
        credit = true
    )
end

examples_nav = fix_suffix.(joinpath.("examples", examples))

makedocs(;
    modules=[Chordal],
    authors="Theo Diamandis",
    repo="https://github.com/tjdiamandis/Chordal.jl/blob/{commit}{path}#L{line}",
    sitename="Chordal.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://tjdiamandis.github.io/Chordal.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => examples_nav,
        "API Reference" => "api.md"
    ],
)

deploydocs(;
    repo="github.com/tjdiamandis/Chordal.jl",
    devbranch = "main"
)


text:
module ClimateTools

# External modules
using Reexport
@reexport using ClimateBase
using NetCDF
@reexport using NCDatasets
using Shapefile
using AxisArrays
using NaNMath
const axes = Base.axes
using ArgCheck
using DataFrames
using Interpolations
using ProgressMeter
using Polynomials
using IterTools
using Statistics
using Random
using Dates
using GeoStats
# using InverseDistanceWeighting
using Extremes
using Distances
using PyCall
using RollingFunctions
import Base.vcat
import Base.getindex
import Base.show
import Base.size
import Base.setindex!
import Base.similar
import Base.write
import Statistics.minimum
import Statistics.maximum
import Statistics.std
import Statistics.var
import Statistics.mean
import Statistics.quantile
import Base: +
import Base: -
import Base: *
import Base: /
import Base.findmax

const scipy = PyNULL()

function __init__()
    copy!(scipy, pyimport_conda("scipy.interpolate", "scipy"))
end

# Included files
include("functions.jl")
include("indices.jl")
include("indicators.jl")
include("extract.jl")
include("interface.jl")
include("cf_conventions.jl")
include("biascorrect.jl")
include("export.jl")
include("time.jl")
include("spatial.jl")
include("analysis.jl")

export inpoly, inpolygrid, meshgrid, inpolyvec, ndgrid
export findmax, findmin
export frostdays, summerdays, icingdays, tropicalnights
export daysabove10 #, daysbelow0, degdaysabove, degdaysbelow
export RXday
export vaporpressure, approx_surfacepressure
export customthresover, customthresunder, annualmax, annualmin
export annualmean, annualsum, prcp1
export wbgt, diurnaltemperature, meantemperature
export drought_dc
export ensemble_mean, ensemble_std, ensemble_max, ensemble_min
export load, load2D
export regrid, applymask
export griddata
export shapefile_coords, shapefile_coords_poly
export resample, spatialsubset, timestep
export qqmap, qqmaptf
export biascorrect_extremes
export permute_west_east
export getdim_lat, getdim_lon, getdim_tim, isdefined, extractpoly
export get_dimname, get_calendar
export polyfit, polyval
export @isdefined
export merge
export minimum, maximum, std, var, mean
export quantile
export daymean, daysum
export monthmean, monthsum, temporalmean
export yearmonthdayhour
export write, findmindist

end #module


text:
## Iterate ###################################################################

"""
    Iterate(iter) :: Distribution{eltype(iter)}

Create a distribution yielding one by one elements from iterator `iter`,
in order.

# Examples
```jldoctest
julia> rand(Fill(Iterate(Iterators.countfrom(1)), 3), 2)
2-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [1, 2, 3]

julia> rand(Iterate(1:9), 10)
ERROR: ArgumentError: iterator exhausted
[...]
```

!!! note
    `Iterate(iter)` is semantically equivalent to the following definition:
    ```julia
    function Iterate(iter)
        local st
        Thunk() do
            it = @isdefined(st) ? iterate(iter, st) :
                                  iterate(iter)
            it === nothing && throw(ArgumentError("iterator exhausted"))
            st = it[2]
            Pure(it[1])
        end
    end
    ```
"""
struct Iterate{T,I} <: Distribution{T}
    iter::I

    Iterate(iter) = new{eltype(iter),typeof(iter)}(iter)
end

mutable struct SamplerIterate{T,I,E,S} <: SamplerReset{T}
    iter::I
    elem::Union{Some{E},Nothing} # Some only on first iteration
    state::S
end

function Sampler(::Type{RNG}, it::Iterate{T,I},
                 n::Repetition) where {RNG<:AbstractRNG,T,I}
    r = iterate(it.iter)
    if r === nothing
        SamplerIterate{T,Nothing,Nothing,Nothing}(nothing, nothing, nothing)
    else
        x, st = r
        SamplerIterate{T,I,typeof(x),typeof(st)}(it.iter,Some(x),st)
    end
end

function reset!(sp::SamplerIterate, _...)
    x, st = iterate(sp.iter) # we know it's not nothing
    sp.elem = Some(x)
    sp.state = st
    sp
end

reset!(sp::SamplerIterate{T,Nothing}, _...) where {T} = sp # nothing to do

function rand(rng::AbstractRNG, sp::SamplerIterate)
    if sp.elem !== nothing
        elem = sp.elem
        sp.elem = nothing
        return something(elem)
    end
    r = iterate(sp.iter, sp.state)
    r === nothing &&
        throw(ArgumentError("iterator exhausted"))
    sp.state = r[2]
    r[1]
end

rand(rng::AbstractRNG, sp::SamplerIterate{T,Nothing}) where {T} =
    throw(ArgumentError("iterator exhausted"))

## SubSeq ####################################################################

"""
    SubSeq(A::AbstractArray{T}, p::Real) where {T} :: Distribution{Vector{T}}

Create a distribution yielding vectors consisting of a random subsequence of
the given array `A`, where each element of `A` is included (in order) with
independent probability `p`. (Complexity is linear in `p*length(A)`, so this
is efficient even if `p` is small and `A` is large.) Technically,
this process is known as "Bernoulli sampling" of `A`.

# Examples
```julia-repl
julia> rand(SubSeq(1:8, 0.3))
3-element Array{Int64,1}:
 3
 4
 7

julia> rand!(Int[], SubSeq(1:8, 0.3), Val(1))
2-element Array{Int64,1}:
 4
 5
```

!!! note
    `SubSeq` is implemented as a thin wrapper over `Random.randsubseq!`, but
    is semantically equivalent to the following definition:
    ```julia
    SubSeq(A, p) = Lift(filter,
                        Bind(Fill(Bernoulli(p), length(A))) do selected
                            Pure(_ -> pop!(selected))
                        end,
                        Pure(A))
    ```
"""
struct SubSeq{T,A} <: Distribution{T}
    a::A
    p::Float64

    SubSeq(a::AbstractArray{T}, p::Real) where {T} =
        new{Vector{T},typeof(a)}(a, p)
end

rand!(rng::AbstractRNG, s::AbstractVector,
      sp::SamplerTrivial{<:SubSeq}, ::Val{1}) =
    randsubseq!(rng, s, sp[].a, sp[].p)

rand(rng::AbstractRNG, sp::SamplerTrivial{<:SubSeq{T}}) where {T} =
    rand!(rng, T(), sp, Val(1))

## SubIter ###################################################################

"""
    SubIter(iter, p::Real) :: Distribution{eltype(iter)}

Create a distribution yielding one by one some elements
from iterator `iter`, in order. Each element is included with
independent probability `p`. Note that
+ `SubIter(iter, 1)` is equivalent to `Iterate(iter)`
+ `rand(SubIter(array, p), n)` is somewhat equivalent to
  `rand(SubSeq(array, p))[1:n]`, assuming `n` is small enough so that
  not all the values of `array` are consumed.

# Examples
```julia-repl
julia> rand(Fill(SubIter(Iterators.countfrom(1), .2), 10), 2)
2-element Array{Array{Int64,1},1}:
 [4, 5, 6, 13, 17, 21, 23, 41, 42, 48]
 [2, 6, 12, 28, 49, 51, 58, 60, 65, 66]
```

!!! note
    `SubIter(iter, p)` is semantically equivalent (roughly) to:
    ```
    SubIter(iter, p) =
        Keep(x -> x !== nothing,
             Bind(Zip(Bernoulli(p), Iterate(iter))) do (succ, x)
                 Pure(succ ? x : nothing)
             end)
    ```
    This is assuming that `iter` doesn't contain `nothing`.
"""
struct SubIter{T,I} <: Distribution{T}
    iter::I
    p::Float64

    SubIter(iter, p::Real) = new{eltype(iter),typeof(iter)}(iter, p)
end

Sampler(::Type{RNG}, sub::SubIter, n::Repetition) where {RNG<:AbstractRNG} =
    SamplerTag{typeof(sub)}((sp = Sampler(RNG, Iterate(sub.iter), n),
                             bernouilli = Sampler(RNG, Bernoulli(sub.p), n)))

reset!(sp::SamplerTag{<:SubIter}, n...) = (reset!(sp.data.sp, n...); sp)

function rand(rng::AbstractRNG, sp::SamplerTag{<:SubIter{T}})::T where {T}
    it = sp.data.sp
    bernouilli = sp.data.bernouilli
    while true
        x = rand(rng, it)
        rand(rng, bernouilli) && return x
    end
end

## Shuffling #################################################################

### Fisher-Yates

struct FisherYates{T,A} <: Distribution{T}
    a::A

    FisherYates(a::AbstractArray{T}) where {T} =
        new{T,typeof(a)}(a)
end

struct SamplerFisherYates{T,A} <: SamplerReset{T}
    a::A
    inds::Vector{Int}
end

fy_inds(a) = Vector{Int}(undef, length(a) + 1)

Sampler(::Type{RNG}, fy::FisherYates{T,A}, ::Repetition) where {RNG<:AbstractRNG} where {T,A} =
    reset!(SamplerFisherYates{T,A}(fy.a, fy_inds(fy.a)))

fy_reset!(inds, n) = @inbounds inds[end] = -n # < 0 means not yet initialized

reset!(sp::SamplerFisherYates, n=length(sp.a)) = (fy_reset!(sp.inds, n); sp)

@noinline function fy_initialize!(rng, inds, k)
    k == 0 &&
        throw(ArgumentError("FisherYates: all elements have been consumed"))
    n = length(inds) - 1
    kn = n
    copyto!(inds, 1:n)
    m = n + k
    mask = nextpow(2, n) - 1
    while n != m
        (mask >> 1) == n && (mask >>= 1)
        i = 1 + rand(rng, Random.ltm52(n, mask))
        #^^^ faster equivalent to i = rand(rng, 1:n) (cf. Base.shuffle!)
        @inbounds inds[i], inds[n] = inds[n], inds[i]
        n -= 1
    end
    kn
end

function fy_rand(rng::AbstractRNG, inds, a)
    @inbounds begin
        k = inds[end] # contains the index in inds where the index in a is located
        if k <= 0
            k = fy_initialize!(rng, inds, k)
        end
        inds[end] = k - 1
        a[inds[k]]
    end
end

rand(rng::AbstractRNG, sp::SamplerFisherYates) = fy_rand(rng, sp.inds, sp.a)

### SelfAvoid

# cf. `self_avoid_sample!` in StatsBase.jl

struct SelfAvoid{T,A} <: Distribution{T}
    a::A

    SelfAvoid(a::AbstractArray{T}) where {T} =
        new{T,typeof(a)}(a)
end

struct SamplerSelfAvoid{T,A,S} <: SamplerReset{T}
    a::A
    seen::Set{Int}
    idx::S
end

sa_idx(::Type{RNG}, a) where {RNG} = Sampler(RNG, Base.OneTo(length(a)), Val(Inf))

function Sampler(RNG::Type{<:AbstractRNG}, sa::SelfAvoid{T,A}, ::Repetition) where {T,A}
    idx  = sa_idx(RNG, sa.a)
    SamplerSelfAvoid{T,A,typeof(idx)}(sa.a, Set{Int}(), idx)
end

sa_reset!(sp) = empty!(sp.seen)

reset!(sp::SamplerSelfAvoid, _=0) = (sa_reset!(sp); sp)

function sa_rand(rng::AbstractRNG, sp)
    seen = sp.seen
    idx = sp.idx
    while true
        i = rand(rng, idx)
        if !(i in seen)
            push!(seen, i)
            return @inbounds sp.a[i]
        end
    end
end

rand(rng::AbstractRNG, sp::SamplerSelfAvoid) = sa_rand(rng, sp)

## Shuffle ###################################################################

struct Shuffle{T,A} <: Distribution{T}
    a::A
    n::Int

    # univariate constructor
    Shuffle(a::AbstractArray{T}) where {T} = new{T,typeof(a)}(a)

    # multivariate constructor
    Shuffle(a::AbstractArray{T}, n::Integer) where {T} =
        new{Vector{T},typeof(a)}(a, n)
end

function variate_size(sh::Shuffle)
    if isunivariate(sh)
        T = gentype(sh)
        T <: Number ? () : throw(ArgumentError("variate_size undefined"))
    else
        (sh.n,)
    end
end

mutable struct SamplerShuffle{T,A,S<:Sampler} <: SamplerReset{T}
    a::A
    n::Int
    alg::Int
    idx::S
    seen::Set{Int}
    inds::Vector{Int}

    SamplerShuffle{TT}(a::AbstractArray{T}, n::Int, idx::S) where {TT,T,S} =
        new{TT,typeof(a),S}(a, n, 0, idx)
end

Sampler(::Type{RNG}, sh::Shuffle{T}, n::Repetition) where {RNG<:AbstractRNG,T} =
    SamplerShuffle{T}(sh.a, sh.n, sa_idx(RNG, sh.a))

isunivariate(::Union{Shuffle{T,A},SamplerShuffle{T,A}}) where {T,A} = T === eltype(A)

function _reset!(sp::SamplerShuffle, n=-1)::Int
    na = length(sp.a)
    if n < 0
        n = na
    end
    if n < 3
        sp.alg = typemin(Int)
    elseif na < n * 24 # TODO: do benchmarks
        if !isdefined(sp, :inds)
            sp.inds = fy_inds(sp.a)
        end
        fy_reset!(sp.inds, n)
        sp.alg = 3
    else
        if !isdefined(sp, :seen)
            sp.seen = Set{Int}()
        end
        sa_reset!(sp)
        sp.alg = 4
    end
end

function reset!(sp::SamplerShuffle, n...)
    isunivariate(sp) && _reset!(sp, n...)
    sp
end

function rand(rng::AbstractRNG, sp::SamplerShuffle{T}) where T
    isunivariate(sp) || return rand!(rng, T(undef, sp.n), sp, Val(1))
    alg = sp.alg
    if alg < 0
        local i
        while true
            i = rand(rng, sp.idx)
            -i != alg && break
        end
        sp.alg = -i
        return @inbounds sp.a[i]
    elseif alg == 3
        fy_rand(rng, sp.inds, sp.a)
    elseif alg == 4
        sa_rand(rng, sp)
    else
        rand(rng, reset!(sp))
    end
end

function rand!(rng::AbstractRNG, v::AbstractVector, sp::SamplerShuffle{T}, ::Val{1}) where T
    isunivariate(sp) &&
        throw(ArgumentError("can not create vector from univariate Shuffle"))
    n = sp.n
    length(v) != n && resize!(v, n) # conditional, as some vectors don't support resize!
    if n < 3
        n < 1 && return v
        j = rand(rng, sp.idx)
        @inbounds v[1] = sp.a[j]
        n < 2 && return v
        local i
        while true
            i = rand(rng, sp.idx)
            i != j && break
        end
        @inbounds v[2] = sp.a[i]
        return v
    end
    _reset!(sp, sp.n)
    if sp.alg == 3
        for i in eachindex(v)
            @inbounds v[i] = fy_rand(rng, sp.inds, sp.a)
        end
    else # alg == 4
        for i in eachindex(v)
            @inbounds v[i] = sa_rand(rng, sp)
        end
    end
    return v
end


text:
using CartesianGeneticProgramming
using LightGraphs
using MetaGraphs
using TikzGraphs
using TikzPictures
using LaTeXStrings

function to_graph(ind::CGPInd; active_outputs=trues(ind.n_out))
    actives = [n.active for n in ind.nodes]
    actives[1:ind.n_in] .= true
    vids = findall(actives)
    #pos = get_positions(c)
    mg = MetaDiGraph(SimpleDiGraph())
    add_vertices!(mg, length(vids)+sum(active_outputs))#c.nout)
    for i in 1:ind.n_in
        set_prop!(mg, i, :name, LaTeXString(string("\$in_{", i, "}\$")))
        set_prop!(mg, i, :type, 0)
    end
    for vi in (ind.n_in+1):length(vids)
        v = vids[vi]
        n = ind.nodes[v]
        f_name = split(split(repr(n.f), ".")[end], "_")[end]
        #if f_name == "const"
        #    set_prop!(mg, vi, :name, LaTeXString(@sprintf("%0.2f", n.p)))
        #else
        set_prop!(mg, vi, :name, LaTeXString(f_name))
        #end
        set_prop!(mg, vi, :function, n.f)
        set_prop!(mg, vi, :type, 2)
        #set_prop!(mg, vi, :param, n.p)
        cx = findfirst(x-> x==n.x,vids)
        cy = findfirst(x-> x==n.y,vids)
        if cx == cy
            add_edge!(mg, Edge(cx, vi))
            set_prop!(mg, cx, vi, :ci, 3)
        else
            add_edge!(mg, Edge(cx, vi))
            set_prop!(mg, cx, vi, :ci, 1)
            if typeof(cy)!=Nothing
                add_edge!(mg, Edge(cy, vi))
                set_prop!(mg, cy, vi, :ci, 2)
            end
        end
    end
    nid_count = 1
    for o in 1:ind.n_out
        if active_outputs[o]
            nid = length(vids)+nid_count
            set_prop!(mg, nid, :name, LaTeXString(string("\$out_{", o, "}\$")))
            set_prop!(mg, nid, :type, 1)
            oid = findfirst(x->x== ind.outputs[o],vids)
            add_edge!(mg, Edge(oid, nid))
            set_prop!(mg, nid, oid, :ci, 0)
            nid_count += 1
        end
    end
    mg
end

function chromo_draw(ind::CGPInd, file::String="graph.pdf"; active_outputs=trues(ind.n_out))
    mg = to_graph(ind, active_outputs=active_outputs)
    names = map(x->get_prop(mg, x, :name), 1:nv(mg))
    t = TikzGraphs.plot(mg.graph, names)
    TikzPictures.save(TikzPictures.PDF(file), t)
end


text:
include("read_write.jl")
include("pricing_specs.jl")
include("../../src/Dexter.jl")
#include("../../src/dd/construction.jl")

using Dexter

using LightGraphs
#using Gadfly
using GraphPlot
using JuMP
using Clp
using CPLEX

# read problem instance and construct the model
filename = joinpath(@__DIR__, "instances", "C_Input_Mon_s1.txt")
specs = read_pricing(filename)
m1 = create_pricing_model(specs)
m2 = create_pricing_LP(specs)

# This part computes the functional property of separable functions
aux_p = Array{Array{Array{Dexter.UVFuncProp, 1}, 1}, 1}(specs.pricing_constr_num)
for i=1:specs.pricing_constr_num
    aux_p[i] = Array{Array{Dexter.UVFuncProp, 1}, 1}(specs.var_num)
    for j=1:specs.var_num
        a = specs.ub[j]/(specs.constr_d[i,j]^(1/specs.constr_d[i,j]))
        rng1 = Dexter.UVFuncProp(0, 0, 1, 0, a)
        rng2 = Dexter.UVFuncProp(0, 0, 1, a, Inf)
        aux_p[i][j] = [rng1, rng2]
    end
end

bb_tree = Dexter.create_BB_tree(m1, m2)
root_time = @elapsed Dexter.create_BB_root!(bb_tree, prop = aux_p)

# Initializing parameters
time_limit0 = 200.0
iter_max_ecp0 = 1000        # max limit on ecp iterations
max_cut_num0 = 3            # max cut number added at EACH iteration of ecp
lp_to_ip_iter0 = 500         # number of ecp iterations to solve LP before switching to IP (if applicable)
lp_to_ip_tol0 = 0.00001      # constraint violation tolerance below which ecp switches to solving IP (if applicable)
tol_num_ecp0 = 10                # number of last iterations over which the relative objective improvement is computed
obj_improvement_tol0 = 1e-4         # relative objective value improvememtn over the last few (tol_num_ecp) iterations

# Running ECP algorithm
x_val, obj_val, ecp_status, initial_lp_obj, iter_ecp, total_cuts, ecp_time, relative_obj_improvement =
Dexter.ECP!(bb_tree.tree[1].model_info.m, bb_tree.tree[1].dd_info.dd, bb_tree.tree[1].dd_info.con_function, bb_tree.tree[1].dd_info.initial,
time_limit = time_limit0, iter_max_ecp = iter_max_ecp0, max_cut_num = max_cut_num0,
lp_to_ip_iter = lp_to_ip_iter0, lp_to_ip_tol = lp_to_ip_tol0, obj_improvement_tol = obj_improvement_tol0, tol_num_ecp = tol_num_ecp0)

# computing DD widths for all DDs
width = Array{Int, 1}(specs.pricing_constr_num)
for i=1:specs.pricing_constr_num
    width[i] = Dexter.get_width(bb_tree.tree[1].dd_info.dd[i])
end

# write results in txt file
fname = joinpath(@__DIR__, "instances", "Dexter_Result_Int_s1.txt")
write_pricing(fname, time_limit0, iter_max_ecp0, max_cut_num0, lp_to_ip_iter0, lp_to_ip_tol0, obj_val, obj_improvement_tol0, tol_num_ecp0,
 ecp_status, initial_lp_obj, iter_ecp, total_cuts, ecp_time, width, root_time, relative_obj_improvement)

#print("\n\n\n")
#println("Objective value: ", obj_val)
#println("ECP status: ", ecp_status)

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--part2"
    help = "Do part 2"
    action = :store_true
    "--test"
    help = "Run test data"
    action = :store_true
end

parsed_args = parse_args(ARGS, s)
test_suffix = parsed_args["test"] ? "_test" : ""

Extent = @NamedTuple{min::Int64, max::Int64}
Cuboid = @NamedTuple{x::Extent, y::Extent, z::Extent}
struct RebootStep
    on::Bool
    cuboid::Cuboid
end

function Base.show(io::IO, c::Cuboid)
    println(io, "[x=($(c.x.min), $(c.x.max)), y=($(c.y.min), $(c.y.max)), z=($(c.z.min), $(c.z.max)))]")
end

function parseline(l:: String) :: RebootStep
    parts = split(l," ")
    onOff = (parts[1] == "on")
    axes = split(parts[2], ",")

    function getminmax(axis::AbstractString) :: Extent
        p = split(axis, "=")
        mm=map(e -> parse(Int, e), split(p[2],".."))
        (min = mm[1], max = mm[2])
    end
    xyz = map(getminmax, axes)
    RebootStep(onOff, (x= xyz[1], y= xyz[2], z=xyz[3]))
end

function overlap(a::Extent, b::Extent)
    a.min <= b.max && a.max>= b.min
end

function intersect(a::Cuboid, b::Cuboid)
    overlap(a.x, b.x) && overlap(a.y, b.y) && overlap(a.z, b.z)
end

function len(e::Extent)
    BigInt(e.max) - BigInt(e.min) + BigInt(1)
end

function size(c::Cuboid)
    BigInt(len(c.x)) * BigInt(len(c.y)) * BigInt(len(c.z))
end

function covers(e1::Extent, e2::Extent)
    e1.min<=e2.min  && e1.max >= e2.max
end

function contains(c1::Cuboid, c2::Cuboid)
    covers(c1.x, c2.x) && covers(c1.y, c2.y) && covers(c1.z, c2.z)
end

steps = [parseline(l) for l in readlines("input/day_22$test_suffix.txt")]

function remove(cut::Cuboid, from::Cuboid) :: Vector{Cuboid}
    result = []
    if cut.x.min > from.x.min
        s = (x=(min=from.x.min, max=cut.x.min-1), y=from.y, z=from.z)
        from = (x=(min=cut.x.min, max=from.x.max), y=from.y, z=from.z)
        push!(result, s)
    end
    if cut.x.max < from.x.max
        s = (x=(min=cut.x.max+1, max=from.x.max), y=from.y, z=from.z)
        from = (x=(min=from.x.min, max=cut.x.max), y=from.y, z=from.z)
        push!(result, s)
    end
    if cut.y.min > from.y.min
        s = (x=from.x, y=(min=from.y.min, max=cut.y.min-1), z=from.z)
        from = (x=from.x, y=(min=cut.y.min, max=from.y.max), z=from.z)
        push!(result, s)
    end
    if cut.y.max < from.y.max
        s = (x=from.x, y=(min=cut.y.max+1, max=from.y.max), z=from.z)
        from = (x=from.x, y=(min=from.y.min, max=cut.y.max), z=from.z)
        push!(result, s)
    end
    if cut.z.min > from.z.min
        s = (x=from.x, y=from.y, z=(min=from.z.min, max=cut.z.min-1))
        from = (x=from.x, y=from.y, z=(min=cut.z.min, max=from.z.max))
        push!(result, s)
    end
    if cut.z.max < from.z.max
        s = (x=from.x, y=from.y, z=(min=cut.z.max+1, max=from.z.max))
        from = (x=from.x, y=from.y, z=(min=from.z.min, max=cut.z.max))
        push!(result, s)
    end
    result
end

function onstep(ons::Vector{Cuboid}, turning::Cuboid)::Vector{Cuboid}
    # println("onstep")
    result = [turning]
    overlaps = []
    for con in ons
        if !intersect(con, turning)
            push!(result, con) # keep it
        elseif contains(turning, con) # ignore it - I got you bae
            continue;
        else
            push!(overlaps, con)
        end
    end
    ## So now we have a collection of disjoint cuboids (overlaps) each of which intersects with the step.cuboid
    # println("ON: Overlaps - $(length(overlaps))")
    for c in overlaps
        # println(c)
        append!(result, remove(turning, c))
    end
    # println("done")
    result
end

function offstep(ons::Vector{Cuboid}, turning::Cuboid)::Vector{Cuboid}
    # println("offstep")
    if isempty(ons)
        return []
    end
    result = []
    overlaps = []
    for con in ons
        if !intersect(con, turning)
            push!(result, con) # keep it
        elseif contains(turning, con)
            continue;
        else
            push!(overlaps, con)
        end
    end
    # KNOW THE OVERLAPS ARE DISJOINT..
    # println("OFF: Overlaps - $(length(overlaps))")
    for c in overlaps
        append!(result,remove(turning, c))
    end
    result
end

function applystep(ons::Vector{Cuboid}, step::RebootStep)::Vector{Cuboid}
    stepper = step.on ? onstep : offstep
    res = stepper(ons, step.cuboid)
    # show_on(res)
    res
end

function show_on(ons::Vector{Cuboid})
    cubes = []
    for c in ons
        for x in c.x.min:c.x.max
            for y in c.y.min:c.y.max
                for z in c.z.min:c.z.max
                    push!(cubes, (x, y, z))
                end
            end
        end
    end
    sort!(cubes)
    for c in cubes
        println("$(c[1]), $(c[2]), $(c[3])")
    end
end

function part1()
    initarea = (x=(min=-50, max=50), y=(min=-50, max=50), z=(min=-50, max=50))
    relevant::Vector{RebootStep} = filter(s -> intersect(s.cuboid, initarea), steps)
    results = foldl(applystep, relevant; init=Vector{Cuboid}())
    sum(map(size, results))
end

function part2()
    results = foldl(applystep, steps; init=Vector{Cuboid}())
    sum(map(size, results))
end

p = parsed_args["part2"] ? part2 : part1
answer = p()
println("The answer is: $answer")


text:
function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{R, D},
    ::Type{S},
) where {
    R <: PSY.RenewableGen,
    D <: AbstractRenewableDispatchFormulation,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(R, sys)

    if !validate_available_devices(R, devices)
        return
    end

    #Variables
    add_variables!(psi_container, ActivePowerVariable, devices)
    add_variables!(psi_container, ReactivePowerVariable, devices)

    #Constraints
    add_constraints!(
        psi_container,
        RangeConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    add_constraints!(
        psi_container,
        RangeConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(psi_container, devices, model, get_feedforward(model))

    #Cost Function
    cost_function!(psi_container, devices, model, S)

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    model::DeviceModel{R, D},
    ::Type{S},
) where {
    R <: PSY.RenewableGen,
    D <: AbstractRenewableDispatchFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(R, sys)

    if !validate_available_devices(R, devices)
        return
    end

    #Variables
    add_variables!(psi_container, ActivePowerVariable, devices)

    #Constraints
    add_constraints!(
        psi_container,
        RangeConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
        get_feedforward(model),
    )
    feedforward!(psi_container, devices, model, get_feedforward(model))

    #Cost Function
    cost_function!(psi_container, devices, model, S)

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::DeviceModel{R, FixedOutput},
    ::Type{S},
) where {R <: PSY.RenewableGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(R, sys)

    if !validate_available_devices(R, devices)
        return
    end

    nodal_expression!(psi_container, devices, S)

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::DeviceModel{PSY.RenewableFix, D},
    ::Type{S},
) where {D <: AbstractRenewableDispatchFormulation, S <: PM.AbstractPowerModel}
    @warn("The Formulation $(D) only applies to FormulationControllable Renewable Resources, \n Consider Changing the Device Formulation to FixedOutput")

    construct_device!(psi_container, sys, DeviceModel(PSY.RenewableFix, FixedOutput), S)

    return
end

function construct_device!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::DeviceModel{PSY.RenewableFix, FixedOutput},
    ::Type{S},
) where {S <: PM.AbstractPowerModel}
    devices = get_available_components(PSY.RenewableFix, sys)

    if !validate_available_devices(PSY.RenewableFix, devices)
        return
    end

    nodal_expression!(psi_container, devices, S)

    return
end


text:
# See https://github.com/blegat/SwitchedSystems.jl

# Example 2.8 of
# P. Parrilo and A. Jadbabaie
# Approximation of the joint spectral radius using sum of squares
# Linear Algebra and its Applications, Elsevier, 2008, 428, 2385-2402
# Inspired from the construction of:
# Ando, T. and Shih, M.-h.
# Simultaneous Contractibility.
# SIAM Journal on Matrix Analysis & Applications, 1998, 19, 487
# The JSR is √2

@testset "[PJ08] Example 2.8 with $(factory.constructor)" for factory in sdp_factories
    isscs(factory) && continue
    @polyvar x[1:2]
    A1 = [1 0; 1 0]
    A2 = [0 1; 0 -1]
    expected_ub = [√2, 1]
    function testlyap(d, γ, feasible::Bool)
        m = SOSModel(factory)
        @variable m p Poly(monomials(x, 2d))
        # p strictly positive
        q = sum(x.^(2*d))
        @constraint m p >= q
        c1 = @constraint m p(x => A1 * vec(x)) <= γ^(2*d) * p
        c2 = @constraint m p(x => A2 * vec(x)) <= γ^(2*d) * p

        JuMP.optimize!(m)

        if feasible
            @test JuMP.primal_status(m) == MOI.FeasiblePoint
        else
            @test JuMP.dual_status(m) == MOI.InfeasibilityCertificate
            μ1 = JuMP.dual(c1)
            μ2 = JuMP.dual(c2)

            # The dual constraint should work on any polynomial.
            # Let's test it with q
            lhs = dot(μ1, q(x => A1 * vec(x))) + dot(μ2, q(x => A2 * vec(x)))
            rhs = dot(μ1, q) + dot(μ2, q)
            @test 1e-6 * max(abs(lhs), abs(rhs)) + lhs >= rhs
        end
    end
    testlyap(1, √2 - 1e-1, false)
    testlyap(1, √2 + 1e-1, true)
    testlyap(2, 1 - 1e-1, false)
    testlyap(2, 1 + 1e-1, true)
end


text:
@testset "assemble_mass_matrix.jl" begin

end


text:
import CSV
import TimeSeries

@testset "Test Timeseries formats" begin
    formats = [
        (
            IS.TimeseriesFormatYMDPeriodAsColumn,
            joinpath(FORECASTS_DIR, "YMDPeriodAsColumn.csv"),
            nothing,
        ),
        (
            IS.TimeseriesFormatYMDPeriodAsHeader,
            joinpath(FORECASTS_DIR, "YMDPeriodAsHeader.csv"),
            "fake",
        ),
        (
            IS.TimeseriesFormatComponentsAsColumnsNoTime,
            joinpath(FORECASTS_DIR, "ComponentsAsColumnsNoTime.csv"),
            nothing,
        ),
        (
            IS.TimeseriesFormatDateTimeAsColumn,
            joinpath(FORECASTS_DIR, "DateTimeAsColumn.csv"),
            nothing,
        ),
        (
            IS.TimeseriesFormatDateTimePeriodAsColumn,
            joinpath(FORECASTS_DIR, "DateTimePeriodAsColumn.csv"),
            nothing,
        ),
    ]

    for (format, filename, component_name) in formats
        file = CSV.File(filename)
        @test format == IS.get_timeseries_format(file)

        data = IS.read_time_series(filename, component_name)
        @test data isa TimeSeries.TimeArray
    end
end


text:
using NonlinearPDE
using Test

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


text:
module Load2
    include("load2_module.jl")
    using Interpolations
    using .Load2_module

    const LAMBDA = 3.886
    const THRESHOLD = 65.0
    const K = 0.190785707092222
    const XYSIZE = 150

    function construct()
        load2_param = Load2_module.Load2_module_param(LAMBDA, THRESHOLD, K)
        load2_val = Load2_module.Load2_module_variables(
            Array{Float64}(undef, XYSIZE),
            nothing)

        load2_val.x =
        [
            0.000, 0.010, 0.020, 0.030, 0.040, 0.050,
            0.060, 0.080, 0.100, 0.150, 0.200, 0.250,
            0.292, 0.300, 0.333, 0.350, 0.375, 0.400,
            0.417, 0.458, 0.500, 0.542, 0.584, 0.600,
            0.625, 0.667, 0.700, 0.709, 0.750, 0.792,
            0.800, 0.833, 0.875, 0.900, 0.917, 0.958,
            1.000, 1.042, 1.083, 1.100, 1.125, 1.167,
            1.200, 1.208, 1.250, 1.300, 1.400, 1.458,
            1.500, 1.600, 1.667, 1.700, 1.800, 1.875,
            1.900, 2.000, 2.083, 2.200, 2.292, 2.400,
            2.500, 2.600, 2.708, 2.800, 2.918, 3.000,
            3.125, 3.200, 3.333, 3.400, 3.500, 3.542,
            3.600, 3.750, 3.800, 3.960, 4.000, 4.167,
            4.375, 4.583, 4.792, 5.000, 5.209, 5.418,
            5.625, 5.834, 6.042, 6.250, 6.458, 6.667,
            6.875, 7.083, 7.292, 7.500, 7.708, 7.917,
            8.125, 8.333, 8.542, 8.750, 8.958, 9.167,
            9.375, 9.583, 9.792, 10.00, 10.22, 10.44,
            10.67, 10.92, 11.16, 11.43, 11.72, 12.01,
            12.31, 12.63, 12.97, 13.33, 13.72, 14.12,
            14.55, 15.01, 15.48, 16.00, 16.56, 17.14,
            17.78, 18.46, 19.20, 20.00, 20.87, 21.82,
            22.85, 24.00, 25.26, 26.67, 28.24, 30.00,
            32.00, 34.00, 34.29, 36.00, 36.92, 38.00,
            40.00, 45.00, 50.00, 55.00, 60.00, 65.00
        ]

        y =
        [
            1.000, 0.985, 0.972, 0.959, 0.947, 0.935,
            0.924, 0.902, 0.882, 0.835, 0.793, 0.755,
            0.727, 0.721, 0.700, 0.691, 0.675, 0.660,
            0.651, 0.627, 0.607, 0.582, 0.569, 0.562,
            0.552, 0.535, 0.521, 0.518, 0.502, 0.488,
            0.485, 0.475, 0.461, 0.453, 0.449, 0.436,
            0.425, 0.414, 0.406, 0.398, 0.393, 0.382,
            0.375, 0.374, 0.364, 0.353, 0.333, 0.322,
            0.315, 0.297, 0.287, 0.283, 0.268, 0.259,
            0.255, 0.244, 0.234, 0.221, 0.212, 0.202,
            0.193, 0.185, 0.176, 0.170, 0.162, 0.157,
            0.150, 0.145, 0.138, 0.134, 0.130, 0.127,
            0.125, 0.118, 0.116, 0.110, 0.108, 0.102,
            0.0956, 0.0895, 0.0837, 0.0788, 0.0739, 0.0695,
            0.0656, 0.0619, 0.0587, 0.0554, 0.0526, 0.0500,
            0.0473, 0.0450, 0.0430, 0.0408, 0.0389, 0.0371,
            0.0355, 0.0340, 0.0321, 0.0310, 0.0298, 0.0287,
            0.0275, 0.0265, 0.0255, 0.0244, 0.0235, 0.0225,
            0.0216, 0.0206, 0.0198, 0.0189, 0.0180, 0.0171,
            0.0163, 0.0155, 0.0147, 0.0139, 0.0131, 0.0123,
            0.0116, 0.0109, 0.0102, 0.0094, 0.0088, 0.0081,
            0.0075, 0.0069, 0.0064, 0.0058, 0.0053, 0.0048,
            0.0043, 0.0038, 0.0034, 0.0030, 0.0026, 0.0022,
            0.0019, 0.0017, 0.0016, 0.0015, 0.0011, 0.0011,
            0.0011, 0.00079, 0.00063, 0.00049, 0.00039, 0.00032
        ]

        load2_val.yitp = interpolate((load2_val.x,), y, Gridded(Linear()))

        return load2_param, load2_val
    end

    dy0(x, load2_param) = let
        a = -3.0 * load2_param.K ^ (3.0 / load2_param.LAMBDA) * x ^ (3.0 / load2_param.LAMBDA - 1.0)
        return a * (1.0 + (load2_param.K * x) ^ (3.0 / load2_param.LAMBDA)) ^ (-load2_param.LAMBDA - 1.0)
    end

    y0(x, load2_param) = let
        return (1.0 + (load2_param.K * x) ^ (3.0 / load2_param.LAMBDA)) ^ -load2_param.LAMBDA
    end

    make_vmax(x2, load2_param, load2_val) = let
        return x2 > load2_param.THRESHOLD ? dy0(x2, load2_param) : Interpolations.gradient(load2_val.yitp, x2)[1]
    end

    load2(x2, v2, load2_param, load2_val) = let
        return [x2 > load2_param.THRESHOLD ? y0(x2, load2_param) : load2_val.yitp(x2), v2]
    end
end


text:
"""
	Parameter(value::Real, name::Symbol)

Defines a parameter value (scalar), with an input value and a name.

# Examples
```jldoctest
julia> Parameter(3.14, :π)
Parameter(3.14, :π)
```
"""
struct Parameter <: DeterministicUQInput
    value::Real
    name::Symbol
end

function sample(p::Parameter, n::Integer=1)
    return DataFrame(p.name => ones(n) * p.value)
end

to_standard_normal_space!(p::Parameter, x::DataFrame) = nothing
to_physical_space!(p::Parameter, x::DataFrame) = nothing

mean(p::Parameter) = DataFrame(p.name => p.value)

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq