text:
#=------------------------------------------------------------------------------
                        Formatting Routines
------------------------------------------------------------------------------=#
function produce_ssten_from_triangles(file;use_metis=false,lcc=false)

    A = MatrixNetworks.readSMAT(file)
    (n,m) = size(A)
    if(n != m)
        println("rectangular")
    end

	if !issymmetric(A)
		A = max.(A,A')  #symmetrize for Triangles routine
    end

    if use_lcc
        A, _ = largest_component(A)
    end

	if use_metis
		apply_Metis_permutation!(A)
	end

    T = collect(MatrixNetworks.triangles(A))

	alterfilename = (file,postfix)-> split(file,".smat")[1]*postfix
	if use_metis
	    output_file = alterfilename(file,"_with_metis.ssten")
	else
	    output_file = alterfilename(file,".ssten")
	end

    open(output_file,"w") do f
        write(f,"$(3)\t$(n)\t$(length(T))\n")

        for (v_i,v_j,v_k) in T
            write(f,"$(v_i)\t$(v_j)\t$(v_k)\t1.0\n")
        end
    end

end

function load_UnweightedThirdOrderSymTensor(filepath;enforceFormatting = true)

	#check path validity
	@assert filepath[end-5:end] == ".ssten"

    open(filepath) do file
		#preallocate from the header line
		order, n, m =
			[parse(Int,elem) for elem in split(chomp(readline(file)),'\t')]
        @assert order == 3

		Ti = [ Vector{Tuple{Int,Int}}(undef, 0) for i in 1:n ]

		i = 1
		@inbounds for line in eachline(file)
			entries = split(chomp(line),'\t')

			if enforceFormatting
				(ti,tj,tk) = sort([parse(Int,elem) for elem in entries[1:end-1]])
			else
				(ti,tj,tk) = [parse(Int,elem) for elem in entries[1:end-1]]
			end

			if 0 == ti || 0 == tj || 0 == tk
				error("elements must be indexed by 1.")
			end
			push!(Ti[ti], (tj,tk))
			push!(Ti[tj], (ti,tk))
			push!(Ti[tk], (ti,tj))

		end

		sort!.(Ti)
		return UnweightedThirdOrderSymTensor(n,Ti)
	end

end
"""------------------------------------------------------------------------------
  Loads in a ThirdOrderSymTensor from an ssten file. Data specifications can be
  found in the 'formatting_specification.info' file the data/ folder. The
  enforceFormatting can be used to ensure that the indices are sorted in
  increasing order, and any files which have 0's in indices are updated to be
  indexed by 1.
------------------------------------------------------------------------------"""
function load_ThirdOrderSymTensor(filepath;enforceFormatting = true)

	#check path validity
	@assert filepath[end-5:end] == ".ssten"

    open(filepath) do file
		#preallocate from the header line
		order, n, m =
			[parse(Int,elem) for elem in split(chomp(readline(file)),'\t')]
        @assert order == 3

		indices = Array{Int,2}(undef,m,order)
		values = Array{Float64,1}(undef,m)

		i = 1
		@inbounds for line in eachline(file)
			entries = split(chomp(line),'\t')

			indices[i,:] = [parse(Int,elem) for elem in entries[1:end-1]]
			if enforceFormatting
				sort!(indices[i,:])
			end
			values[i] = parse(Float64,entries[end])
			i += 1
		end

		#check for 0 indexing
		zero_indexed = false

		@inbounds for i in 1:m
		    if indices[i,1] == 0
    			zero_indexed = true
				break
			end
	    end

		if zero_indexed
			indices .+= 1
		end

        return ThirdOrderSymTensor(n,indices,values)
    end
end

function write_to_armadillo_raw_ascii_format(X::Array{T,2},output_file::String) where T
    open(output_file,"w") do f
        for row in eachrow(X)
            println(f,join(row," "))
        end
    end
end
#=------------------------------------------------------------------------------
                        Local File Experiments Routines
------------------------------------------------------------------------------=#

#TODO: update to use kwargs properly
function self_alignment(ssten_files::Array{String,1};method="LambdaTAME")

    Best_alignment_ratio = Array{Float64}(undef,length(ssten_files))

    for f in ssten_files
        A = load(f;type="COOTen")

        #permute tensor to align against
        p = sort(1:A.cubical_dimension,by=i->rand())
        A_permuted = COOTen(A.indices,A.vals,A.cubical_dimension)
        permute_tensor!(A_permuted,p)

        best_score, _, _ = align_tensors(A,A_permuted;method=method)

        println("$f aligned with $best_score")
    end

end

function align_tensors(graph_A_file::String,graph_B_file::String;
                       ThirdOrderSparse=true,profile=false,
					   kwargs...)

    if ThirdOrderSparse
        A = load_ThirdOrderSymTensor(graph_A_file)
        B = load_ThirdOrderSymTensor(graph_B_file)
    else
        A = load(graph_A_file,false,"COOTen")
        B = load(graph_B_file,false,"COOTen")
    end

	if profile
		return align_tensors_profiled(A,B;kwargs...)
	else
		return align_tensors(A,B;kwargs...)
	end
end

function apply_Metis_permutation!(A::SparseMatrixCSC,k::Int=100)

	n = size(A,1)
	metis_partition = Metis.Metis.partition(A, 100)
	p = sort(1:n, by= i-> metis_partition[i])

	permute!(A,p,p)
end

#TODO: put in kwargs
function pairwise_alignment(dir)

    #align all .ssten files
    ssten_files = sort([f for f in readdir(dir) if occursin(".ssten",f)])

    Best_alignment_ratio = Array{Float64}(undef,length(ssten_files),length(ssten_files))

    for i in 1:length(ssten_files)
        for j in i+1:length(ssten_files)
            Best_alignment_ratio[i,j] = align_from_files(dir*ssten_files[i],dir*ssten_files[j])
            Best_alignment_ratio[i,j] = Best_alignment_ratio[j,i]
        end
    end

    return ssten_files, Best_alignment_ratio
end

#=
function distributed_pairwise_alignment(dir::String;kwargs...)

    #align all .ssten files
    ssten_files = sort([f for f in readdir(dir) if occursin(".ssten",f)])
    distributed_pairwise_alignment(ssten_files,dir;kwargs...)

end
=#

"""------------------------------------------------------------------------------
  This function runs alignemnt experiments in parallel using the @spawn macro of
  Distributed.jl. This function should be run in the src file, as the @everywhere
  macro will only work properly there. Any additional kwargs are passed to
  'align_matrices' or 'align_tensors' depending on the file extensions.

  Inputs
  ------
  * files - (Array{String,1}):
    List of the file names. Files may be either all SMAT or all STEN files.
    'align_tensors' or 'align_matrices' will be called depending on the extension.
  * dirpath - (String):
    The path to where the files exist.
  * method - (String):
    What type of alignment method to run, options are 'LambdaTAME', 'TAME',
    'LowRankTAME', 'EigenAlign', 'LowRankEigenAlign', 'Degree', or 'Random'.
  * profile - (Bool):
    Whether or not to profile the algorithm. Setting this flag to true will
    record different things depending on the algorithm, but will always return
    runtime information.
Output
------
  * 'exp_results' - (type omitted):
    The final array of results for the series of experiments, the size of each
    entry will differ depending on whether or not profile was True. At the very
    least, each entry will contain the the pair of file names aligned, the
    number of triangles matched in the alignemnt, and the maximum number of
    triangles possible to align (the minimum number of triangles in either
    network).

------------------------------------------------------------------------------"""
function distributed_pairwise_alignment(files::Array{String,1},dirpath::String;
                                        method="LambdaTAME",profile=false,kwargs...)

    @everywhere include_string(Main,$(read("LambdaTAME.jl",String)),"LambdaTAME.jl")

    alignment_object = ""
    if all( [f[end-5:end] == ".ssten" for f in files]) #tensors
        alignment_object = "Tensors"
    elseif all( [f[end-4:end] == ".smat" for f in files]) #matrices
        alignment_object = "Matrices"
    else
        throw(ArgumentError("all files must be the same file type, either all '.ssten' or '.smat'."))
    end

    futures = []
	exp_results = []

#    Best_alignment_ratio = Array{Float64}(undef,length(ssten_files),length(ssten_files))

    for i in 1:length(files)
        for j in i+1:length(files)

            if alignment_object == "Tensors"
                future = @spawn align_tensors(dirpath*files[i],dirpath*files[j];profile=profile,method=method,kwargs...)
            elseif alignment_object == "Matrices"
                A = MatrixNetworks.readSMAT(dirpath*files[i])
                B = MatrixNetworks.readSMAT(dirpath*files[j])
                future = @spawn align_matrices(A,B;profile=profile,method=method,kwargs...)
            end
            push!(futures,((i,j),future))
        end
    end

    for ((i,j), future) in futures

		if method == "LambdaTAME" ||  method == "LowRankTAME"
			if profile
				matched_tris, max_tris, _, _, best_matching, results = fetch(future)
			else
				matched_tris, max_tris, _, _, best_matching = fetch(future)
			end
		elseif method == "TAME"
			if profile
				matched_tris, max_tris, _, best_matching, results = fetch(future)
			else
				matched_tris, max_tris, _, best_matching = fetch(future)
            end
        elseif method == "EigenAlign" || method == "Degree" || method == "Random" || method == "LowRankEigenAlign"
            A_tens, B_tens, matched_tris, best_matching, results = fetch(future)
            max_tris = min(A_tens, B_tens)
            profile = true
		end
		if profile
			push!(exp_results,(files[i],files[j],matched_tris, max_tris, best_matching, results))
		else
			push!(exp_results,(files[i],files[j],matched_tris, max_tris, best_matching))
		end
    end

    return exp_results
end

function self_alignment(dir::String;kwargs...)

    #align all .ssten files
    ssten_files = sort([dir*"/"*f for f in readdir(dir) if occursin(".ssten",f)])
    return self_alignment(ssten_files;kwargs...)

end

#=------------------------------------------------------------------------------
                      Generated Graph Experiments Routines
------------------------------------------------------------------------------=#
"""------------------------------------------------------------------------------
  This function runs randomly generated experiments in parallel using the @spawn
  macro of Distributed.jl. This function should be run in the src file, as the
  @everywhere macro will only work properly there. Any additional kwargs are
  passed to 'random_graph_exp'.

  Inputs
  ------
  * 'trial_count' - (Int):
    The number of trials to run for each of the experiments conducted.
  * 'seed_exps' - (Bool):
    Indicates whether or not to seed the experiments. Seeds are generated by
    seeding Random.seed!(0), and then generating an array of UInts using the
    number of 'p_remove' and 'n_sizes' variables, and the 'trial_count'. Changing
    these may lead to different seeds.
  * 'n_sizes' - (Array{Int,1}):
    The different sizes of graphs to generate experiments for.
  * 'p_remove' - (Array{Float,1}):
    The different probability of edge removals to test in the ER noise models.
  * profile - (Bool):
    Whether or not to profile the algorithm. Setting this flag to true will
    record different things depending on the algorithm, but will always return
    runtime information.
Output
------
  * 'exp_results' - (type omitted):
    The final array of results for the series of experiments, the size of each
    entry will differ depending on whether or not profile was True. At the very
    least, the each entry will contain the seed used, the 'p_remove' value, the
    number of vertices in the graph, the ground truth accuracy of the alignment,
    the degree weighted accuracy of the alignment, the number of matched
    triangles, the number of triangles each graph, and the maximum number of
    triangles possible to align. Anything additional returned by the functions
    with the profile flag on, is appended at the end.
------------------------------------------------------------------------------"""
function distributed_random_trials(trial_count::Int,seed_exps::Bool=false
                                   ;method="LambdaTAME",graph_type::String="ER",
								   n_sizes = [100, 500, 1000, 2000,5000], p_remove = [.01,.05],
								   profile=false,kwargs...)

    #only handling even batch sizes
    #@assert trial_count % process_count == 0

    #ensure file is loaded on all processes
    @everywhere include_string(Main,$(read("LambdaTAME.jl",String)),"LambdaTAME.jl")

	if seed_exps
		Random.seed!(0)
		seeds = rand(UInt64, length(p_remove),length(n_sizes), trial_count)
	end

	seed = nothing
	results = []
    p_index = 1

    futures = []

    for p in p_remove

        n_index = 1

        for n in n_sizes

            for trial = 1:trial_count

                if seed_exps
                    seed = seeds[p_index,n_index,trial]
                end
                future = @spawn random_graph_exp(n,p,graph_type;
                                                    profile=profile,seed=seed,method=method,
                                                    kwargs...)
                push!(futures,(seed,p,n,future))

            end

            n_index += 1
        end
        p_index += 1
    end

    for (seed,p,n,future) in futures
        if profile
            if method == "LambdaTAME" ||  method == "LowRankTAME"
                d_A, d_B, perm, (A_tris, B_tris,(matched_tris, max_tris, _, _,best_matching, exp_results))= fetch(future)
            elseif method == "TAME"
                d_A, d_B, perm, (A_tris, B_tris,(matched_tris, max_tris, _, best_matching, exp_results))= fetch(future)
            end

            accuracy = sum([1 for (i,j) in enumerate(perm) if get(best_matching,j,-1) == i])/n
            D_A = sum(d_A)
            D_B = sum(d_B)
            degree_weighted_accuracy = sum([(get(best_matching,j,-1) == i) ? ((d_A[i] + d_B[j])/(D_A+D_B)) : 0.0 for (i,j) in enumerate(perm)])

            push!(results,(seed, p, n, accuracy, degree_weighted_accuracy, matched_tris, A_tris, B_tris, max_tris, exp_results))
        else
            if method == "LambdaTAME" ||  method == "LowRankTAME"
                d_A, d_B, perm, (A_tris, B_tris,(matched_tris, max_tris, _, _, best_matching)) = fetch(future)
            elseif method == "TAME"
                d_A, d_B, perm, (A_tris, B_tris, (matched_tris, max_tris, _, best_matching)) = fetch(future)
            elseif method == "EigenAlign" || method == "Degree" || method == "Random" || method == "LowRankEigenAlign"
                d_A, d_B, perm, (A_tris, B_tris, matched_tris, best_matching, _) = fetch(future)
                max_tris = minimum((A_tris,B_tris))
            end

            accuracy = sum([1 for (i,j) in enumerate(perm) if get(best_matching,j,-1) == i])/n
            D_A = sum(d_A)
            D_B = sum(d_B)
            degree_weighted_accuracy = sum([(get(best_matching,j,-1) == i) ? ((d_A[i] + d_B[j])/(D_A+D_B)) : 0.0 for (i,j) in enumerate(perm)])

            push!(results,( seed, p, n, accuracy, degree_weighted_accuracy, matched_tris, A_tris, B_tris, max_tris))
        end
    end

    return results

end

"""-----------------------------------------------------------------------------
  Runs an instance of a graph alignment problem using a random graph model. Once
  a graph is generated, we use the function 'ER_noise_model' to generate a
  second network to align against it. Once the second network is generated it is
  randomly permuted. All additional kwargs are passed to the 'align_matrices'
  function.

  Inputs
  ------
  * n - (Int):
    The node of nodes to use in each network.
  * 'p_remove' -(Float):
    The probability of removing an edge in the ER noise model. The probability
    of adding in an edge is computed using a function of this parameter at the
    density of the generated network.
  * 'graph_type' - (String):
    The type of graph used to generate the network. Options include 'ER' for an
    Erdos Renyi network, 'RandomGeometric' for a Random Geometric network from
    a 2D unit square, and 'HyperKron' which comes from the Triangle Generalized
    Preferential Attachment model.
  * seed - (Any):
    A seed used to replicate the results. Our experiments typically use UInts.
  * 'use_metis' - (Bool):
    Indicates whether or not to use metis clusters to permute the adjacency
    matrix. For large networks this can improve the runtime by improving the
    access patterns.
  * degreedist - (UnivariateDistribution):
    A random variable generator from Distributions.jl which is used to generate
    the degrees for random nodes. We primarily use
    Distributions.LogNormal(log(5,1)).
  * 'p_edges' - (func Int - Float):
    An anonymous function used for computing the probability of adding edges.
    Primarily used in the ER graphs for generating the probability of adding
    in edges from the number of nodes in the network. an example expected input
    would be 'n -> 2*log(n)/n'.

    Outputs
    -------
    * 'd_A','d_B' - (Array{Float,1}):
      The degrees of all the nodes in graphs A and B. The ith entry of each
      being the degree of node i in the respective networks.
    * perm - (Array{Int,1}):
      The permutation used to shuffle the nodes in graph B. Used for computing
      the precision of the methods.
    * All additional return variables are from 'align_matrices' function call.
-----------------------------------------------------------------------------"""
function random_graph_exp(n::Int, p_remove::Float64,graph_type::String;
                          seed=nothing,use_metis=false,degreedist=nothing,p_edges=nothing,kwargs...)

	if seed !== nothing
		Random.seed!(seed)
	end

    if graph_type == "ER"
        if degreedist === nothing
            if p_edges === nothing
                p = 2*log(n)/n
            else
                p = p_edges(n)
            end
            A = erdos_renyi(n,p)
        else
            println(degreedist)
            A = erdos_renyi(n;degreedist)
            p = nnz(A)/n^2
        end
	elseif graph_type == "RandomGeometric"

		if degreedist === nothing
            k = 10
            if p_edges === nothing
                p = k/n
            end
			A = random_geometric_graph(n,k)
		else
			d = 2
			A = spatial_network(n, d;degreedist= degreedist)
			p = nnz(A)/n^2
		end

    elseif graph_type == "HyperKron"
        if p_edges === nothing
            p = .4#2*log(n)/n
        end
		r = .4

		A = sparse(gpa_graph(n,p,r,5))

	else
		error("invalid graph type: $graph_type\n must be either 'ER','RandomGeometric' or 'HyperKron'")
	end

    p_add = (p*p_remove)/(1-p)
	B = ER_noise_model(A,n,p_remove,p_add)

    perm = shuffle(1:n)
    B = B[perm,perm]

	if use_metis
		apply_Metis_permutation!(A,100)
		apply_Metis_permutation!(B,100)
	end

    d_A = A*ones(n)
    d_B = B*ones(n)

    return d_A,d_B,perm,align_matrices(A,B;kwargs...)

end

function ER_noise_model(A,n::Int,p_remove::F,p_add::F) where {F <: AbstractFloat}
    B = copy(A)

    is,js,_ = findnz(erdos_renyi(n,p_remove))
    for (i,j) in zip(is,js)
        if A[i,j] > 0
            B[i,j] = 0
            B[j,i] = 0
        end
    end

    is,js,_ = findnz(erdos_renyi(n,p_add))
    for (i,j) in zip(is,js)
        if A[i,j] == 0.0
            B[i,j] = 1
            B[j,i] = 1
        end
    end

    return B
end

"""------------------------------------------------------------------------------
  This function takes in a pair of sparse matrices and aligns them with with the
  specified method. If the method requested is a tensor method, the adjacency
  tensor is built from the triangles enumerated. Additional kwargs are passed

  Inputs
  ------
  * A, B - (SparseMatrixCSC{T,Int}):
    The sparse matrices to be aligned.
  * profile - (Bool):
    Whether or not to profile the methods run. The type of profiling done is
    dependent on the method run.
  * method - (String):
    The method to run. Options include 'Degree', 'Random', 'LowRankEigenAlign',
    'LambdaTAME', 'LowRankTAME', "TAME". Currently supporting 'EigenAlign', but
    may remove in the future as LowRankEigenAlign produces the same result but
    with much better scalability.

  Outputs
  -------
  output depends on whether a tensor method is used or not. Matrix based methods
  return the trianlges in A, triangles in B, triangles matched, the matching
  used, runtime.

  If a tensor method is used, the number of triangles in A and B are returned, in
  addition to whatever is returned by 'align_tensors(_profiled)'.
------------------------------------------------------------------------------"""
function align_matrices(A::SparseMatrixCSC{T,Int},B::SparseMatrixCSC{T,Int};profile=false,kwargs...) where T

    A_ten = graph_to_ThirdOrderTensor(A)
    B_ten = graph_to_ThirdOrderTensor(B)

    if kwargs[:method] == "LambdaTAME" || kwargs[:method] == "LowRankTAME" || kwargs[:method] == "TAME"
        if profile
            return size(A_ten.indices,1), size(B_ten.indices,1), align_tensors_profiled(A_ten,B_ten;kwargs...)
        else
            return size(A_ten.indices,1), size(B_ten.indices,1), align_tensors(A_ten,B_ten;kwargs...)
        end

    elseif kwargs[:method] == "EigenAlign" || kwargs[:method] == "Degree" || kwargs[:method] == "Random" || kwargs[:method] == "LowRankEigenAlign"

        if kwargs[:method] == "LowRankEigenAlign"
            iters = 10
            (ma,mb,_,_),t = @timed align_networks_eigenalign(A,B,iters,"lowrank_svd_union",3)
            matching = Dict{Int,Int}([i=>j for (i,j) in zip(ma,mb)])
        elseif kwargs[:method] == "EigenAlign"
            (ma,mb),t = @timed NetworkAlignment.EigenAlign(A,B)
            matching = Dict{Int,Int}(zip(ma,mb))
        elseif kwargs[:method] == "Degree"
            (ma,mb),t = @timed degree_based_matching(A,B)
            matching = Dict{Int,Int}(zip(ma,mb))
        elseif kwargs[:method] == "Random"
            n,n = size(B)

            matching,t = @timed Dict{Int,Int}(enumerate(shuffle(1:n)))
        else
            error("Invalid input, must be ")
        end
        triangle_count, gaped_triangles, _ = TAME_score(A_ten,B_ten,matching)
        return size(A_ten.indices,1), size(B_ten.indices,1), triangle_count, matching, t
    else
        throw(ArgumentError("method must be one of 'LambdaTAME', 'LowRankTAME', 'TAME', 'EigenAlign', 'LowRankEigenAlign', 'Degree', or 'Random'."))
    end

end

"""-----------------------------------------------------------------------------
  This function takes in a sparse matrix and builds a ThirdOrderSymTensor from
  the triangles present in the networks. If the network isn't symmetric, then
  it's symmetrized. The bool 'use_lcc' is used to indicate whether or not to find
  the largest strongly connected component before symmetrizing.
-----------------------------------------------------------------------------"""
function graph_to_ThirdOrderTensor(A;use_lcc=false)

    n,n = size(A)

    if use_lcc
        A,_ = largest_component(A)
    end

    if !issymmetric(A)
        println("Symmetrizing matrix")
		A = max.(A,A')  #symmetrize for Triangles routine
	end

    #build COOTens from graphs
    tris = collect(MatrixNetworks.triangles(A))
    nnz = length(tris)

    indices = Array{Int,2}(undef,nnz,3)
    for i =1:nnz
        indices[i,1]= tris[i][1]
        indices[i,2]= tris[i][2]
        indices[i,3]= tris[i][3]
    end

    vals = ones(nnz)

    return ThirdOrderSymTensor(n,indices,vals)

end

# --------------------------------------------------------------------------- #
#						       Random Graph Models
# --------------------------------------------------------------------------- #

function erdos_renyi(n,p)
    A = sprand(n,n,p)
    U = triu(A,1)
    is,js,_ = findnz(max.(U,U'))
    return sparse(is,js,ones(length(is)),n,n)
end

function erdos_renyi(n;degreedist=LogNormal(log(5),1))

    println("made it")
    # form the edges for sparse
    ei = Int[]
    ej = Int[]

    for i=1:n
      deg = ceil(Int,minimum((rand(degreedist),n-1)))
      neighbors = sample(1:n,deg+1, replace = false)

      for j in neighbors
        if i != j
          push!(ei,i)
          push!(ej,j)
        end
      end
    end

    A = sparse(ei,ej,1.0,n,n)

    return max.(A,A')
end

#--                             David's graph code                          --#

function random_geometric_graph(n,k)
  xy = rand(2,n)
  T = BallTree(xy)
  idxs = knn(T, xy, minimum((k,n)))[1]
  # form the edges for sparse
  ei = Int[]
  ej = Int[]
  for i=1:n
    for j=idxs[i]
      if i != j
        push!(ei,i)
        push!(ej,j)
      end
    end
  end
  A = sparse(ei,ej,1.0,n,n)

  return max.(A,A')
end

function spatial_network(n::Integer, d::Integer; degreedist=LogNormal(log(5),1))
  xy, ei, ej = spatial_graph_edges(n, d;degreedist=degreedist)
  A = sparse(ei,ej,1.0,n,n)
  return max.(A,A')
end

function spatial_graph_edges(n::Integer,d::Integer;degreedist=LogNormal(log(5),1))
  xy = rand(d,n)
  T = BallTree(xy)
  # form the edges for sparse
  ei = Int[]
  ej = Int[]
  for i=1:n
    deg = ceil(Int,minimum((rand(degreedist),n-2)))
    idxs, dists = knn(T, xy[:,i], deg+1)
    for j in idxs
      if i != j
        push!(ei,i)
        push!(ej,j)
      end
    end
  end
  return xy, ei, ej
end


text:
@testset "gather" begin
    T = Float32
    CT = CuArray{Float32}

    ## 1d src, 2d index of ints -> 2d output
    src = CT([3, 4, 5, 6, 7])
    index = cu([1 2 3 4;
                4 2 1 3;
                3 5 5 3])
    output = CT([3 4 5 6;
                6 4 3 5;
                5 7 7 5])

    y = NNlib.gather(src, index)
    @test y isa CuArray{Float32,2}
    @test size(y) == size(index)
    gputest(src -> NNlib.gather(src, index), src, checkgrad=true)
    @test NNlib.gather!(CUDA.zeros(T, size(index)...), src, index) == output
    @test_throws ArgumentError NNlib.gather!(zeros(T, 3, 5), src, index)

    ## 1d src, 2d index of tuples -> 2d output
    src = CT([3, 4, 5, 6, 7])
    index = cu([(1,) (2,) (3,) (4,);
                (4,) (2,) (1,) (3,);
                (3,) (5,) (5,) (3,)])
    output = CT([3 4 5 6;
                6 4 3 5;
                5 7 7 5])

    y = NNlib.gather(src, index)
    @test y isa CuArray{Float32,2}
    @test size(y) == size(index)
    gputest(src -> NNlib.gather(src, index), src, checkgrad=true)
    @test NNlib.gather!(CUDA.zeros(T, size(index)...), src, index) == output
    @test_throws ArgumentError NNlib.gather!(zeros(T, 3, 5), src, index)

    ## 1d src, 2d index of CartesianIndex -> 2d output
    src = CT([3, 4, 5, 6, 7])
    index = cu(CartesianIndex.([(1,) (2,) (3,) (4,);
                (4,) (2,) (1,) (3,);
                (3,) (5,) (5,) (3,)]))
    output = CT([3 4 5 6;
                6 4 3 5;
                5 7 7 5])

    y = NNlib.gather(src, index)
    @test y isa CuArray{Float32,2}
    @test size(y) == size(index)
    gputest(src -> NNlib.gather(src, index), src, checkgrad=true)
    @test NNlib.gather!(CUDA.zeros(T, size(index)...), src, index) == output
    @test_throws ArgumentError NNlib.gather!(zeros(T, 3, 5), src, index)

    ## 1d src, 3d index of ints -> 3d output
    src = CT([3, 4, 5, 6, 7])
    index = cu([1 2 3 4;
                4 2 1 3;
                3 5 5 3][:,:,1:1])
    output = CT([3 4 5 6;
                6 4 3 5;
                5 7 7 5][:,:,1:1])

    y = NNlib.gather(src, index)
    @test y isa CuArray{Float32,3}
    @test size(y) == size(index)
    gputest(src -> NNlib.gather(src, index), src, checkgrad=true)

    ## 2d src, 2d index of ints -> 3d output
    src = CT([3 5 7
             4 6 8])
    index = cu([1 2 3;
                2 2 1;
                3 1 3])

    output = zeros(T, 2, 3, 3)

    output[:,:,1] = [3 5 7
                    4 6 8]

    output[:,:,2] = [5 5 3
                    6 6 4]

    output[:,:,3] = [7 3 7
                    8 4 8]

    y = NNlib.gather(src, index)
    M = NNlib.typelength(eltype(index))
    Nsrc = ndims(src)
    @test y isa CuArray{Float32,3}
    @test size(y) == (size(src)[1:Nsrc-M]..., size(index)...)
    gputest(src -> NNlib.gather(src, index), src, checkgrad=true)
end

ext: jl


text:
using BinaryBuilder, Pkg

name = "libgeotiff"
version = v"1.6.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/OSGeo/libgeotiff/releases/download/$version/libgeotiff-$version.tar.gz", "9311017e5284cffb86f2c7b7a9df1fb5ebcdc61c30468fb2e6bca36e4272ebca"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""

cd $WORKSPACE/srcdir/libgeotiff-*/

# Who knows why they want to install Windows binaries not build with MSVC to
# ${prefix} instead of ${prefix}/bin
atomic_patch -p1 ../patches/more-reasonable-default-bin-subdir.patch

mkdir build && cd build

# Hint to find libstc++, required to link against C++ libs when using C compiler
if [[ "${target}" == *-linux-* ]]; then
    if [[ "${nbits}" == 32 ]]; then
        export CFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib"
    else
        export CFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib64"
    fi
fi

cmake -DCMAKE_INSTALL_PREFIX=$prefix \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_BUILD_TYPE=Release \
      -DBUILD_SHARED_LIBS=ON \
      ..

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libgeotiff", :libgeotiff),
    ExecutableProduct("makegeo", :makegeo),
    ExecutableProduct("geotifcp", :geotifcp),
    ExecutableProduct("listgeo", :listgeo),
    ExecutableProduct("applygeo", :applygeo),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="PROJ_jll", uuid="58948b4f-47e0-5654-a9ad-f609743f8632")),
    Dependency(PackageSpec(name="Libtiff_jll", uuid="89763e89-9b03-5906-acba-b20f662cd828")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
"""
Abstract type for a `Bijector{N}` making use of auto-differentation (AD) to
implement `jacobian` and, by impliciation, `logabsdetjac`.
"""
abstract type ADBijector{AD, N} <: Bijector{N} end

# AD implementations
function jacobian(b::Union{B, Inversed{B}}, x::Real) where {B<:ADBijector{<:ForwardDiffAD}}
    return ForwardDiff.derivative(b, x)
end
function jacobian(
    b::Union{B, Inversed{B}},
    x::AbstractVector{<:Real}
) where {B<:ADBijector{<:ForwardDiffAD}}
    return ForwardDiff.jacobian(b, x)
end

function jacobian(b::Union{B, Inversed{B}}, x::Real) where {B<:ADBijector{<:TrackerAD}}
    return Tracker.data(Tracker.gradient(b, x)[1])
end
function jacobian(
    b::Union{B, Inversed{B}},
    x::AbstractVector{<:Real}
) where {B<:ADBijector{<:TrackerAD}}
    # We extract `data` so that we don't return a `Tracked` type
    return Tracker.data(Tracker.jacobian(b, x))
end

struct SingularJacobianException{B<:Bijector} <: Exception
    b::B
end
function Base.showerror(io::IO, e::SingularJacobianException)
    print(io, "jacobian of $(e.b) is singular")
end

# TODO: allow batch-computation, especially for univariate case?
"Computes the absolute determinant of the Jacobian of the inverse-transformation."
function logabsdetjac(b::ADBijector, x::Real)
    res = log(abs(jacobian(b, x)))
    return isfinite(res) ? res : throw(SingularJacobianException(b))
end

function logabsdetjac(b::ADBijector, x::AbstractVector{<:Real})
    fact = lu(jacobian(b, x), check=false)
    return issuccess(fact) ? log(abs(det(fact))) : throw(SingularJacobianException(b))
end


text:
# NOTE: You pass in weights using a row vector, e.g. [500 300 100], then I convert it to a 1-dim column vector
function make_dense_nn(count_weights::Array{Int64,1},count_classes::Int)
    count_weights = vec(count_weights)
    count_hidden_layers = size(count_weights, 1)
    count_inner_product_layers = count_hidden_layers + 1

    # Make model name string
    model_name = "dense"
    for i = 1:count_hidden_layers
        model_name = model_name * "_" * string(count_weights[i])
    end

    # Returns a function that returns a function to create model
    nn_factory = function()
        common_layers = Array(InnerProductLayer, count_inner_product_layers)

        # Create the input, hidden, and output layers
        for i = 1:count_inner_product_layers
            # Specify the from and to Symbols
            if i == 1
                from_Symbol = :data
            else
                from_Symbol = Symbol( "ip$(i-1)" )
            end
            to_Symbol = Symbol( "ip$(i)" )
            # Create the layer
            if i != count_inner_product_layers
                common_layers[i] = InnerProductLayer(
                                        name = string(to_Symbol),
                                        output_dim = count_weights[i],
                                        neuron = Neurons.ReLU(),
                                        bottoms = [from_Symbol],
                                        tops = [to_Symbol],
                                        weight_regu = NoRegu())
            else
                common_layers[i] = InnerProductLayer(
                                            name = string(to_Symbol),
                                            output_dim = count_classes,
                                            bottoms = [from_Symbol],
                                            tops = [to_Symbol],
                                            weight_regu = NoRegu())
            end
        end

        last_dense = Symbol("ip$(count_inner_product_layers)")

        pred_prob_layer = SoftmaxLayer( name = "pred_probs", bottoms = [last_dense],tops= [:pred_probs])

        storage_layer = MemoryOutputLayer( name = "storage", bottoms = [:pred_probs])

        loss_layer = SoftmaxLossLayer(name = "loss",
                                    bottoms = [Symbol( "ip$(count_inner_product_layers)" ),:label],
				                    #normalize = :no ) # <= I don't think normalize is necessary...
                                    )

        acc_layer = AccuracyLayer( name = "test-accuracy",
                               bottoms = [Symbol( "ip$(count_inner_product_layers)" ), :label] )
        return (common_layers, pred_prob_layer, storage_layer, loss_layer, acc_layer)
    end
    return (nn_factory, model_name)
end


text:
struct SBool{V} <: Integer
    function SBool{V}() where V
        !(typeof(V) === Bool) && throw(ArgumentError("SBool only supports static Bool storage, got $(typeof(V))"))
        new{V}()
    end
end

const SInteger{V} = Union{SSigned{V},SUnsigned{V},SBool{V}}

function SInteger(val::Val{V}) where V
    if V isa Unsigned
        if V isa UInt8
            SUInt8{V}()
        elseif V isa UInt16
            SUInt16{V}()
        elseif V isa UInt32
            SUInt32{V}()
        elseif V isa UInt128
            SUInt128{V}()
        else
            SUInt64{V}()
        end
    elseif V isa Bool
        SBool{V}()
    elseif V isa Signed
        if V isa Int8
            SInt8{V}()
        elseif V isa Int16
            SInt16{V}()
        elseif V isa Int32
            SInt32{V}()
        elseif V isa Int128
            SInt128{V}()
        else
            SInt64{V}()
        end
    end
end

SInteger(x::Signed) = SSigned(x)

SInteger(x::Unsigned) = SUnsigned(x)

Base.unsigned(x::SBool{true}) = SUInt(1)
Base.unsigned(x::SBool{false}) = SUInt(0)

Base.Unsigned(x::SInteger) = Base.unsigned(x)

Base.show(io::IO, v::SInteger) = print(io, "SInteger($(values(v)))")

SI2BI = Dict(SUInt128 => UInt128,
             SUInt16 => UInt16,
             SUInt32 => UInt32,
             SUInt64 => UInt64,
             SUInt8 => UInt8,

             SInt128 => Int128,
             SInt16 => Int16,
             SInt32 => Int32,
             SInt64 => Int64,
             SInt8 => Int8)

BI2SI = Dict(UInt128 => SUInt128,
             UInt16 => SUInt16,
             UInt32 => SUInt32,
             UInt64 => SUInt64,
             UInt8 => SUInt8,

             Int128 => SInt128,
             Int16 => SInt16,
             Int32 => SInt32,
             Int64 => SInt64,
             Int8 => SInt8)

SIntegerZeroType = Union{typeof(SInt8Zero),typeof(SUInt8Zero),typeof(SInt16Zero),
                         typeof(SUInt16Zero),typeof(SInt32Zero),typeof(SUInt32Zero),
                         typeof(SInt64Zero),typeof(SUInt64Zero),typeof(SInt128Zero),
                         typeof(SUInt128Zero)}
SIntegerOneType = Union{typeof(SInt8One),typeof(SUInt8One),typeof(SInt16One),typeof(SUInt16One),
                        typeof(SInt32One),typeof(SUInt32One),typeof(SInt64One),typeof(SUInt64One),
                        typeof(SInt128One),typeof(SUInt128One)}
SIntegerNegOneType = Union{SInt8{-Int8(1)},SInt16{-Int16(1)},SInt32{-Int32(1)},SInt64{-1},SInt128{-Int128(1)}}

function defint(::Type{ST}, ::Type{BT}) where {ST,BT<:BaseInteger}
    @eval begin
        Base.:(~)(::$ST{X}) where X = $ST{~X}()
        Base.:(|)(::$ST{X}, ::$ST{Y}) where {X,Y} = $ST{(|)(X::$BT, Y::$BT)}()
        Base.xor(::$ST{X}, ::$ST{Y}) where {X,Y} = $ST{xor(X::$BT, Y::$BT)}()
        Base.:(&)(::$ST{X}, ::$ST{Y}) where {X,Y} = $ST{(&)(X::$BT, Y::$BT)}()

        Base.trailing_zeros(::$ST{X}) where X = SInt{trailing_zeros(X::$BT)}()
        Base.trailing_ones(::$ST{X}) where X = SInt{trailing_ones(X::$BT)}()
        Base.count_ones(::$ST{X}) where X = SInt{count_ones(x::$BT)}()
        Base.leading_zeros(::$ST{X}) where X = SInt{leading_zeros(X::$BT)}()
        Base.leading_ones(::$ST{X}) where X = SInt{leading_ones(X::$BT)}()

        Base.flipsign(::$ST{X}, ::$ST{Y}) where {X,Y} = $ST{flipsign(X::$BT, Y::$BT)}()
        Base.flipsign(::$ST{X}, ::SSigned{Y}) where {X,Y}= $ST{flipsign(X::$BT, Y)}()

        Base.rem(x::$ST, ::Type{<:$ST{<:Any}}) = x

        Base.powermod(::$ST{X}, ::$ST{P}, ::$ST{M}) where {X,P,M} = $ST{powermod(X::$BT, P::$BT, M::$BT)}()

        @generated function Base.binomial(::$ST{N}, ::$ST{K}) where {N,K}
            x = $ST{binomial(N::$BT, K::$BT)}()
            :($x)
        end

        @generated function Base.factorial(n::$ST{N}) where N
            x = $ST{factorial(N::$BT)}()
            :($x)
        end
        @generated function Base.sqrt(::$ST{X}) where X
            x = SFloat64{sqrt(X::$BT)}()
            :($x)
        end
        Base.decompose(x::$ST) = x, SInt{0}(), SInt{1}()
    end
end

for (S,B) in SI2BI
    defint(S,B)
end

#=
struct SBool{V} <: Integer
    function SBool{V}() where V
        !(typeof(V) === Bool) && throw(ArgumentError("val must be of type Bool"))
        new{V}()
    end
    SBool(v::Bool) = SBool{v,Bool}()
end

SBool(::Val{V}) where V = SBool{V}()

Base.typemin(::Type{SBool}) = SBool{false}()
Base.typemax(::Type{SBool}) = SBool{true}()

(!)(::SBool{X}) where X = SBool{!X::Bool}()

(~)(::SBool{X}) where X = SBool{!X::Bool}()
(&)(::SBool{X}, ::SBool{Y}) where {X,Y} = SBool{(&)(X::Bool, Y::Bool)}()
(|)(::SBool{X}, y::SBool{Y}) where {X,Y} = SBool{(|)(X::Bool, Y::Bool)}()

xor(x::SBool{X}, y::SBool{Y}) where {X,Y} = (X::Bool != Y::Bool)

>>(x::SBool{X}, c::SUInt) where X = SInt(X::Bool) >> c
<<(x::SBool{X}, c::SUInt) where X = SInt(X::Bool) << c
>>>(x::SBool{X}, c::SUInt) where X = SInt(X::Bool) >>> c

Base.signbit(x::SBool) = SBool{false}()
Base.sign(x::SBool) = x
Base.abs(x::SBool) = x
Base.abs2(x::SBool) = x
Base.iszero(x::SBool) = !x
Base.isone(x::SBool) = x

<(x::SBool, y::SBool) = y &! x
<=(x::SBool, y::SBool) = y |! x

## do arithmetic as Int ##

+(x::SBool) =  SInt(x)
-(x::SBool) = -SInt(x)

+(x::SBool, y::SBool) = SInt(x) + SInt(y)
-(x::SBool, y::SBool) = SInt(x) - SInt(y)
*(x::SBool, y::SBool) = x & y
^(x::SBool, y::SBool) = x | !y
^(x::Union{SSigned,SUnsigned}, y::SBool) = ifelse(y, x, one(x))

# preserve -0.0 in `false + -0.0`
+(x::SBool, y::T) where T<:AbstractFloat = ifelse(x, oneunit(y) + y, y)

+(y::AbstractFloat, x::SBool) = x + y

*(y::AbstractFloat, x::SBool) = x * y

Base.div(x::SBool, y::SBool) = y ? x : throw(DivideError())
Base.fld(x::SBool, y::SBool) = div(x,y)
Base.cld(x::SBool, y::SBool) = div(x,y)
Base.rem(x::SBool, y::SBool) = y ? false : throw(DivideError())
Base.mod(x::SBool, y::SBool) = rem(x,y)

const SInteger{V} = Union{SSigned{V},SUnsigned{V},SBool{V}}
=#
#>>(::SInteger{V1}, v2::BaseInteger) where V1 = >>(V1, v2)
#>>(v2::BaseInteger, ::SInteger{V1}) where V1 = >>(v2, V1)

#<<(::SInteger{V1}, v2::BaseInteger) where V1 = <<(V1, v2)
#<<(v2::BaseInteger, ::SInteger{V1}) where V1 = <<(v2, V1)

#>>(x::SInteger, y::SInt) = ifelse(0 <= y, x >> unsigned(y), x << unsigned(-y))

#=
function >>(x::SInteger, c::SInteger)
    Base.@_inline_meta
    typemin(SInt) <= c <= typemax(SInt) && return x >> (c % SInt)
    (x >= 0 || c < 0) && return zero(x)
    oftype(x, -SOne)
end
=#


text:
# Autogenerated wrapper script for Cbc_jll for i686-w64-mingw32-libgfortran3-cxx03
export cbc, libCbc, libOsiCbc, libcbcsolver

using ASL_jll
using Cgl_jll
using Clp_jll
using Osi_jll
using CoinUtils_jll
using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Cbc")
JLLWrappers.@declare_library_product(libCbc, "libCbc-3.dll")
JLLWrappers.@declare_library_product(libOsiCbc, "libOsiCbc-3.dll")
JLLWrappers.@declare_library_product(libcbcsolver, "libCbcSolver-3.dll")
JLLWrappers.@declare_executable_product(cbc)
function __init__()
    JLLWrappers.@generate_init_header(ASL_jll, Cgl_jll, Clp_jll, Osi_jll, CoinUtils_jll, OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libCbc,
        "bin\\libCbc-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libOsiCbc,
        "bin\\libOsiCbc-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcbcsolver,
        "bin\\libCbcSolver-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        cbc,
        "bin\\cbc.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module ProjectTests

using Test
using GCP
using ..GCPTests

@testset "Project" begin

  @test begin
    proj = GCP.Project(GCPTests.CredFilename, ["cloud-platform"])
    proj != nothing
  end

end

end # module ProjectTests

max_stars_repo_name: ProfLeao/mais

max_issues_repo_name: ProfLeao/mais

max_forks_repo_name: ProfLeao/mais


text:
module TestSolubility

using Base.Test
using SoilTracers.PhysChem: solub_co2, solub_cos

@time @testset "Solubility functions" begin
    temp_0C = 273.15
    temp_25C = 298.15

    # test solub_co2

    k_co2_0C = solub_co2(temp_0C)
    @test typeof(k_co2_0C) <: Real
    @test isapprox(k_co2_0C, 1.739, rtol=1e-3)

    k_co2_25C = solub_co2(temp_25C)
    @test typeof(k_co2_25C) <: Real
    @test isapprox(k_co2_25C, 0.831, rtol=1e-3)

    # test solub_cos

    k_cos_0C = solub_cos(temp_0C)
    @test typeof(k_cos_0C) <: Real
    @test isapprox(k_cos_0C, 1.549, rtol=1e-3)

    k_cos_25C = solub_cos(temp_25C)
    @test typeof(k_cos_25C) <: Real
    @test isapprox(k_cos_25C, 0.488, rtol=1e-3)
end

end


text:
"""
    IViewable(nexid)

An interface for widgets that can be shown or hidden.
"""
struct IViewable <: AbstractINextion
    _nid::NexID

    function IViewable(nid::NexID)
        new(nid)
    end
end

"""
    obj.visible = new_val

Set value from `new_val` to Nextion object `obj`.
"""
function setproperty!(obj::IViewable, property::Symbol, new_val)

    if property == :visible
        nid = NexID(obj)
        _oid = String(Name(nid))
        #_oid = String(ComponentID(nid))
        if new_val
            cmd = "vis $_oid,1"
        else
            cmd = "vis $_oid,0"
        end
        nexSerial = NexSerial(nid)
        send(nexSerial, cmd)
    else
        setfield!(obj, property, new_val)
    end

end


text:
export super_square_lattice

 function super_square_lattice(size::NTuple{5, Int})
    m, um, n, un, t = size
    old = LinearIndices((1:t, 1:un, 1:n, 1:um, 1:m))

    Dict(
            old[k, uj, j, ui, i] => (i, j)
            for i=1:m, ui=1:um, j=1:n, uj=1:un, k=1:t
    )
end

function super_square_lattice(size::NTuple{3, Int})
    m, n, t = size
    super_square_lattice((m, 1, n, 1, t))
end


text:
import ReverseDiff

# Gen versions of certain ReverseDiff functions
# this is because the generated code for the static IR is currently loaded in
# the Main module. to avoid requiring the Main module environment to import the
# name ReverseDiff, these aliases allow us to only require that the Main module
# environmment have Gen imported.
new_tape() = ReverseDiff.InstructionTape()
track(value, tape) = ReverseDiff.track(value, tape)
value(maybe_tracked) = ReverseDiff.value(maybe_tracked)
deriv!(maybe_tracked, deriv) = ReverseDiff.deriv!(maybe_tracked, deriv)
reverse_pass!(tape) = ReverseDiff.reverse_pass!(tape)
deriv(tracked) = ReverseDiff.deriv(tracked)
istracked(maybe_tracked) = ReverseDiff.istracked(maybe_tracked)
record!(tape, instruction_type, args...) = ReverseDiff.record!(tape, instruction_type, args...)
increment_deriv!(arg, deriv) = ReverseDiff.increment_deriv!(arg, deriv)
seed!(tracked) = ReverseDiff.seed!(tracked)
unseed!(tracked) = ReverseDiff.unseed!(tracked)

using ReverseDiff: InstructionTape, TrackedReal, SpecialInstruction, TrackedArray

########
# fill #
########

function Base.fill(x::TrackedReal{V}, dims::Integer...) where {V}
    tp = ReverseDiff.tape(x)
    out = ReverseDiff.track(fill(ReverseDiff.value(x), dims...), V, tp)
    ReverseDiff.record!(tp, SpecialInstruction, fill, (x, dims), out)
    return out
end

@noinline function ReverseDiff.special_reverse_exec!(
        instruction::SpecialInstruction{typeof(fill)})
    x, dims = instruction.input
    output = instruction.output
    ReverseDiff.istracked(x) && ReverseDiff.increment_deriv!(x, sum(ReverseDiff.deriv(output)))
    ReverseDiff.unseed!(output)
    return nothing
end

@noinline function ReverseDiff.special_forward_exec!(
        instruction::SpecialInstruction{typeof(fill)})
    x, dims = instruction.input
    ReverseDiff.value!(instruction.output, fill(value(x), dims...))
    return nothing
end


text:
using Gadfly

plot(x=[0, 1, 1, 2, 2, 3, 3, 2, 2, 1, 1, 0, 4, 5, 5, 4],
     y=[0, 0, 1, 1, 0, 0, 3, 3, 2, 2, 3, 3, 0, 0, 3, 3],
     group=["H", "H", "H", "H", "H", "H", "H", "H",
            "H", "H", "H", "H", "I", "I", "I", "I"],
     Geom.polygon(preserve_order=true, fill=true))


text:
"""
$(TYPEDEF)
"""
struct StandardBVProblem end

"""
$(TYPEDEF)
"""
struct BVProblem{uType,tType,isinplace,P,F,bF,PT,K} <: AbstractBVProblem{uType,tType,isinplace}
    f::F
    bc::bF
    u0::uType
    tspan::tType
    p::P
    problem_type::PT
    kwargs::K
    @add_kwonly function BVProblem{iip}(f::AbstractODEFunction,bc,u0,tspan,p=NullParameters(),
                            problem_type=StandardBVProblem();
                            kwargs...) where {iip}
        _tspan = promote_tspan(tspan)
        new{typeof(u0),typeof(tspan),iip,typeof(p),
                  typeof(f),typeof(bc),
                  typeof(problem_type),typeof(kwargs)}(
                  f,bc,u0,_tspan,p,
                  problem_type,kwargs)
    end

    function BVProblem{iip}(f,bc,u0,tspan,p=NullParameters();kwargs...) where {iip}
        BVProblem(convert(ODEFunction{iip},f),bc,u0,tspan,p;kwargs...)
      end
end

function BVProblem(f::AbstractODEFunction,bc,u0,tspan,args...;kwargs...)
    BVProblem{DiffEqBase.isinplace(f,4)}(f,bc,u0,tspan,args...;kwargs...)
end

function BVProblem(f,bc,u0,tspan,p=NullParameters();kwargs...)
    BVProblem(convert(ODEFunction,f),bc,u0,tspan,p;kwargs...)
end

# convenience interfaces:
# Allow any previous timeseries solution
function BVProblem(f,bc,sol::T,tspan,p=NullParameters();kwargs...) where {T<:AbstractTimeseriesSolution}
    BVProblem(f,bc,sol.u,tspan,p)
end
# Allow a function of time for the initial guess
function BVProblem(f,bc,initialGuess::T,tspan::AbstractVector,p=NullParameters();kwargs...) where {T}
    u0 = [ initialGuess( i ) for i in tspan]
    BVProblem(f,bc,u0,(tspan[1],tspan[end]),p)
end

"""
$(TYPEDEF)
"""
struct TwoPointBVPFunction{bF}
    bc::bF
end
TwoPointBVPFunction(; bc = error("No argument bc")) = TwoPointBVPFunction(bc)
(f::TwoPointBVPFunction)(residual, ua, ub, p) = f.bc(residual, ua, ub, p)
(f::TwoPointBVPFunction)(residual, u, p) = f.bc(residual, u[1], u[end], p)

"""
$(TYPEDEF)
"""
struct TwoPointBVProblem{iip} end
function TwoPointBVProblem(f,bc,u0,tspan,p=NullParameters();kwargs...)
    iip = DiffEqBase.isinplace(f,4)
    TwoPointBVProblem{iip}(f,bc,u0,tspan,p;kwargs...)
end
function TwoPointBVProblem{iip}(f,bc,u0,tspan,p=NullParameters();kwargs...) where {iip}
    BVProblem{iip}(f,TwoPointBVPFunction(bc),u0,tspan,p;kwargs...)
end


text:
#****************************************************************#
# VeysiADN 12 May 2021
# Advanced Robotics Homework 5 Codes
#****************************************************************#
using LinearAlgebra
#****************************************************************#
# Returns Adjoint matrix with given rigid body
# transformation matrix g
function Adjoint(g)

    R=g[1:3,1:3];
    p=g[1:3,4];
    p_hat = [ 0       -p[3]     p[2];
             p[3]       0      -p[1];
            -p[2]      p[1]       0 ]

    return [R zeros(3,3);p_hat*R R]
end

max_stars_repo_name: veysiadn/AdvancedRoboticsJulia

max_issues_repo_name: veysiadn/AdvancedRoboticsJulia

max_forks_repo_name: veysiadn/AdvancedRoboticsJulia


text:
function setup_menu(gui, game)
    mb = GtkMenuBar(name="menu_bar")
    file = GtkMenuItem("_File")
    file_menu = GtkMenu(file, name="file_menu")
    new_game = GtkMenuItem("New Game", name="new_game")

    signal_connect(x->start_new_game!(x, gui, game), new_game, :activate)
    push!(file_menu, new_game)
    setup = GtkMenuItem("Setup", name="setup")

    signal_connect(x->setup_game(x, game, gui), setup, :activate)
    push!(file_menu, setup)
    push!(mb, file)
    return mb
end

start_new_game!(x, gui, game) = start_new_game!(gui, game)

function start_new_game!(gui, game)
    game.score = 0
    game.round = game.n_rounds
    rand_color!.(game.dots)
    map(x->x.selected = false, game.dots)
    clear_selected!(game)
    remove_components!(gui)
    filename = joinpath(@__DIR__, "style.css")
    style = CssProviderLeaf(;filename)
    generate_gui!(game, gui, style)
end

function setup_game(x, game, gui)
    popup = GtkWindow("Setup")
    base_panel = GtkBox(:v)
    push!(popup, base_panel)
    grid = GtkGrid()
    col_label = GtkLabel("Number of Columns")
    row_label = GtkLabel("Number of Rows")
    round_label = GtkLabel("Number of Rounds")
    col_entry = GtkEntry(name="col_value")
    row_entry = GtkEntry(name="row_value")
    round_entry = GtkEntry(name="round_value")
    grid[1,1] = col_label
    grid[2,1] = col_entry
    grid[1,2] = row_label
    grid[2,2] = row_entry
    grid[1,3] = round_label
    grid[2,3] = round_entry
    GAccessor.justify(col_label, Gtk.GConstants.GtkJustification.LEFT)
    GAccessor.justify(row_label, Gtk.GConstants.GtkJustification.LEFT)
    set_gtk_property!(grid, :column_spacing, 5)  # introduce a 15-pixel gap between columns
    set_gtk_property!(grid, :row_spacing, 5)
    set_gtk_property!(grid, :column_homogeneous, true)
    set_gtk_property!(grid, :row_homogeneous, true)
    push!(base_panel, grid)
    ok_button = GtkButton("OK")
    components = (col_entry=col_entry, row_entry=row_entry, round_entry=round_entry, popup=popup)
    signal_connect(x-> modify_game(x, components, game, gui), ok_button, "clicked")
    cancel_button = GtkButton("Cancel")
    signal_connect(x->close_window(x, popup), cancel_button, "clicked")
    hbox = GtkButtonBox(:h)
    set_gtk_property!(hbox, :expand, ok_button, true)
    set_gtk_property!(hbox, :spacing, 10)
    push!(base_panel, hbox)
    push!(hbox, cancel_button)
    push!(hbox, ok_button)
    showall(popup)
end

function modify_game(x, c, game, gui)
    str_cols = get_gtk_property(c.col_entry, :text, String)
    str_rows = get_gtk_property(c.row_entry, :text, String)
    str_rounds = get_gtk_property(c.round_entry, :text, String)

    is_error,n_cols = parse_value(str_cols, "Columns")
    is_error ? (return nothing) : nothing

    is_error,n_rows = parse_value(str_rows, "Rows")
    is_error ? (return nothing) : nothing

    is_error,n_rounds = parse_value(str_rounds, "Rounds")
    is_error ? (return nothing) : nothing

    game = Game(; n_rows, n_cols, n_rounds)
    start_new_game!(gui, game)
    close_window(c.popup)
    return nothing
end

function parse_value(value, label)
    number = 0
    try
        number = parse(Int, value)
    catch
        error_popup("$label must be a number")
        return true,-1
    end
    if number < 1
        is_error = true
        error_popup("$label must be greater than 1")
        return true,-1
    end
    return false,number
end

function error_popup(error_message)
    error_window = GtkWindow("Error")
    base_panel = GtkBox(:v)
    set_gtk_property!(base_panel, :spacing, 20)
    message = GtkLabel(error_message)
    ok_button = GtkButton("OK")
    push!(error_window, base_panel)
    push!(base_panel, message)
    push!(base_panel, ok_button)
    signal_connect(x->close_window(x, error_window), ok_button, "clicked")
    showall(error_window)
    return nothing
end

close_window(_, component) = close_window(component)

function close_window(component)
    hide(component)
end

generate_gui!(game, gui) = generate_gui!(game, gui.gui, gui.style)

function generate_gui!(game, gui, style)
    base_panel = GtkBox(:v, name="base_panel")
    menu_bar = setup_menu(gui, game)
    push!(base_panel, menu_bar)

    info_panel = GtkBox(:h, name="info_panel")
    sc = Gtk.GAccessor.style_context(info_panel)
    set_gtk_property!(info_panel, :spacing, 30)
    push!(base_panel, info_panel)
    push!(sc, StyleProvider(style), 600)

    label = GtkLabel("Rounds ", name="rounds_label")
    sc = Gtk.GAccessor.style_context(label)
    push!(sc, StyleProvider(style), 600)
    push!(info_panel, label)

    round_count = GtkLabel(string(game.round), name="round_count")
    sc = Gtk.GAccessor.style_context(round_count)
    push!(sc, StyleProvider(style), 600)
    push!(info_panel, round_count)

    score_label = GtkLabel("Score", name="score_label")
    sc = Gtk.GAccessor.style_context(score_label)
    push!(sc, StyleProvider(style), 600)
    push!(info_panel, score_label)

    score_value = GtkLabel("0", name="score_value")
    sc = Gtk.GAccessor.style_context(score_value)
    push!(sc, StyleProvider(style), 600)
    push!(info_panel, score_value)

    g = GtkGrid(name="grid")
    push!(base_panel, g)

    submit_button = GtkButton(name="submit_button")
    submit_label = GtkLabel("Submit", name="submit_label")
    push!(submit_button, submit_label)
    sc = Gtk.GAccessor.style_context(submit_label)
    push!(sc, StyleProvider(style), 600)

    push!(base_panel, submit_button)
    signal_connect(_->click_submit!(game, gui, style), submit_button, "clicked")

    n_rows,n_cols = size(game.dots)
    for r in 1:n_rows, c in 1:n_cols
        b = GtkButton("")
        dot = game.dots[r,c]
        add_color!(b, dot, style)
        g[c,r] = b
        signal_connect(x->click_dot!(game, dot, gui, x, style), b, "clicked")
    end
    set_gtk_property!(g, :column_spacing, 5)
    set_gtk_property!(g, :row_spacing, 5)
    set_gtk_property!(g, :column_homogeneous, true)
    set_gtk_property!(g, :row_homogeneous, true)
    set_gtk_property!(g, :expand, true)

    push!(gui, base_panel)
    showall(gui)
end

get_button(gui, dot) = gui.gui[1][3][dot.col,dot.row]

function click_dot!(game, dot, gui::GUI)
    button = get_button(gui, dot)
    return click_dot!(game, dot, gui.gui, button, gui.style)
end

function click_dot!(game, dot, gui, button, style)
    !click_dot!(game, dot) ? (return false) : nothing
    if dot.selected
        make_grey!(button, dot, style)
    else
        add_color!(button, dot, style)
    end
    return true
end

function add_color!(button, dot, style)
    label = button[1]
    sc = Gtk.GAccessor.style_context(label)
    push!(sc, StyleProvider(style), 600)
    set_gtk_property!(label, :name, dot.color)
end

function make_grey!(button, dot, style)
    label = button[1]
    sc = Gtk.GAccessor.style_context(label)
    push!(sc, StyleProvider(style), 600)
    set_gtk_property!(label, :name, string(dot.color, "_clicked"))
end

function make_all_grey!(gui, game, style)
    for (button,dot) in zip(gui[1][3],game.dots)
        make_grey!(button, dot, style)
    end
    return nothing
end

click_submit!(game, gui::GUI) = click_submit!(game, gui.gui, gui.style)

function click_submit!(game, gui, style)
    !click_submit!(game) ? (return false) : nothing
    update_round!(game, gui)
    update_score!(game, gui)
    update_all_colors!(game, gui, style)
    game_over!(game) ? make_all_grey!(gui, game, style) : nothing
    return true
end

function update_all_colors!(game, gui, style)
    dots = game.dots
    n_rows,n_cols = size(dots)
    buttons = gui[1][3]
    for r in 1:n_rows
        for c in 1:n_cols
            add_color!(buttons[c,r], dots[r,c], style)
        end
    end
    return nothing
end

function shift_colors!(game, gui, style)
    for dot in game.selected_dots
        shift_color!(dot, game, gui, style)
    end
end

function shift_color!(dot, game, gui, style)
    dots = game.dots
    row = dot.row
    col = dot.col
    buttons = gui[1][3]
    while (row > 1)
        dots[row,col].color = dots[row-1,col].color
        add_color!(buttons[col,row], dots[row,col], style)
        row -= 1
    end
    rand_color!(dots[row,col])
    add_color!(buttons[col,row], dots[row,col], style)
end

function update_round!(game, gui)
    score = gui[1][2][2]
    set_gtk_property!(score, :label, string(game.round))
end

function update_score!(game, gui)
    counter = gui[1][2][4]
    set_gtk_property!(counter, :label, string(game.score))
    return nothing
end

function remove_components!(gui)
    for g in gui
        delete!(gui, g)
    end
end

function start()
    GUI()
end


text:
function FDEsolver(F, tSpan, y0, β, J, par...; h = 0.01, nc = 3, StopIt = "Standard", tol = 10e-10, itmax = 10)

    # Time discretization
    N::Int64 = cld(tSpan[2] - tSpan[1], h)
    t = tSpan[1] .+ collect(0:N) .* h

    # Enter the initial values
    Y = defineY(N, y0, β)

    # Number of equations
    Neq = size(y0, 1)

    # Calculate T with taylor expansion
    m::Int64 = ceil(β[1])

    for n in 1:N

        T0 = taylor_expansion(tSpan[1], t[n], y0, β, m)
        JF = zeros(Neq, Neq)
        Ψ = zeros(Neq)

        if n == 1

            # Y1
            Y1 = T0 .+ h .^ β .* F(t, n, β, Y, par...) ./ Γ(β .+ 1)
            Y = indexY(n, Y, Y1)

            # inverse matrix including Jacobian function
            if Neq == 1

                JF = inv(I(Neq) .- (h .^ β ./ Γ(β .+ 2)) * J(t, n + 1, β, Y, par...))

            else

                JF = inv(I(Neq) .- Diagonal(h .^ β ./ Γ(β .+ 2)) * J(t, n + 1, β, Y, par...))

            end

            if StopIt == "Standard"

                for j in 1:nc

                    # Y11
                    Ψ = (Y[n + 1, :] .- T0 .- h .^ β .* β .* F(t, n, β, Y, par...) ./ Γ(β .+ 2) .- h .^ β .* F(t, n + 1, β, Y, par...) ./ Γ(β .+ 2))
                    Y11 = Y[n + 1, :] - JF * Ψ
                    Y = indexY(n, Y, Y11)

                end

            elseif StopIt == "Convergence"

                σ = 1.1 * tol
                j = 0

                while (σ > tol && j < itmax)

                    # Y11
                    Y11 = T0 .+ h .^ β .* β .* F(t, n, β, Y, par...) ./ Γ(β .+ 2) .+ h .^ β .* F(t, n + 1, β, Y, par...) ./ Γ(β .+ 2)
                    σ = sqrt(sum((Y11 .- Y[2, :]) .^ 2))
                    Y = indexY(n, Y, Y11)

                    j += 1

                end

            end

        else

            ϕ = Phi(Y, F, β, t, n, par...)

            # Yp
            Yp = T0 .+ h .^ β .* (ϕ .- α(0, β) .* F(t, n - 1, β, Y, par...) .+ 2 .* α(0, β) .* F(t, n, β, Y, par...))
            Y = indexY(n, Y, Yp)

            # inverse matrix including Jacobian function
            if Neq == 1

                JF = inv(I(Neq) .- (α(0, β) .* h .^ β) * J(t, n+1, β, Y, par...))

            else

                JF = inv(I(Neq) .- Diagonal(α(0, β) .* h .^ β) * J(t, n+1, β, Y, par...))

            end

            if StopIt == "Standard"

                for j in 1:nc

                    # Y2
                    Ψ = (Y[n + 1, :] .- (T0 .+ h .^ β .* ϕ) .- h .^ β.* α(0, β) .* F(t, n + 1, β, Y, par...))
                    Y2 = Y[n + 1, :] - JF * Ψ
                    Y = indexY(n, Y, Y2)
                end

            elseif StopIt == "Convergence"

                σ = 1.1 * tol
                j = 0

                while (σ > tol && j < itmax)

                    # Y2
                    Y2 = T0 .+ h .^ β .* (ϕ .+ α(0, β) .* F(t, n + 1, β, Y, par...))
                    σ = sqrt(sum((Y2 .- Y[n + 1, :]) .^ 2))
                    Y = indexY(n, Y, Y2)

                    j += 1

                end

            end

        end

    end

    # Output
    t, Y

end


text:
"""
ashcroft.jl

Ashcroft's empty-core pseudopotential
"""

immutable Ashcroft <: PseudoPotential
  z::Vector{Float64} # valence of core ion
  rc::Vector{Float64} # pseudopotential core radius
end

function Ashcroft(z::Real, rc::Real)
  Ashcroft([z], [rc])
end

function formfactor(nfe::NFE{Ashcroft})::Vector{Function}
  ρ = nfe.ρ

  N = length(nfe.pseudo.z)
  ω = Vector{Function}(N)

  for i = 1:N
    z = nfe.pseudo.z[i]
    rc = nfe.pseudo.rc[i]

    ω[i] = q -> -4π*ρ / q^2 * z * cos(q*rc)
  end

  return ω
end


text:
using Plots, LaTeXStrings, Statistics

function Deposition(; len, tot_time, dep_rate, color_step)
    surface = InitialSurf(len)
    AllLens = []
    for n in 0:tot_time
        randsurf = rand(1:len, dep_rate)
        for i in randsurf
            surface = check_and_depose(surface, i, len, (n%color_step) + 1)
        end
        append!(AllLens, Check_Dlen(surface))
    end
    return AllLens
end

# function Check_Dlen(surf)
#     Dlen = []
#     for row in surf
#         lenf = 0
#         for num in 1:length(row)
#             if row[num] == 0
#                 lenf += 1
#             else
#                 break
#             end
#         end
#         lenb = 0
#         for num in length(row):-1:1
#             if row[num] == 0
#                 lenb += 1
#             else
#                 break
#             end
#         end
#         append!(Dlen, length(row) - (lenb+lenf))
#     end
#     return max(Dlen...)
# end

function Check_Dlen(surf)
    Dlen = []
    for row in surf
        yindex = []
        for num in 1:length(row)
            if row[num] != 0
                append!(yindex, num)
            end
        end
        append!(Dlen, max(yindex...) - min(yindex...))
    end
    return max(Dlen...)
end

function check_and_depose(surf, index, len, addval) #checking down and sides and first layer
    i1 , i2 = sides(index, len)
    if surf[end][index] != 0
        newsurf = [0 for j = 1:len]
        newsurf[index] = addval
        push!(surf, newsurf)
    end
    for num in 0:length(surf) -1
        if surf[end - num][index] == 0 && (surf[end - num][i1] != 0 || surf[end - num][i2] != 0)
            surf[end - num][index] = addval
            break
        elseif length(surf) > num + 1 && surf[end - 1 - num][index] != 0
            surf[end - num][index] = addval
            break
        end
    end
    return surf
end

function sides(n, L)
    if n == L
        return n-1 , 1
    elseif n == 1
        return L , n+1
    else
        return n-1, n+1
    end
end

function InitialSurf(len)
    pos1 = [0 for i=1:len]
    pos1[floor(Int, len/2)] = 1
    return [pos1]
end

Parameters = Dict(:len => 300,
                    :tot_time => 50,
                        :dep_rate => 1000,
                            :color_step => 100)

iternum = 200
allVal = [ [0.0 for i in 0:Parameters[:tot_time]] for j = 1:iternum]
meanVal = [0.0 for i in 0:Parameters[:tot_time]]
vars = [0.0 for i in 0:Parameters[:tot_time]]
for i in 1:iternum
    VarList = Deposition(;Parameters...)
    allVal[i] = VarList
    meanVal += VarList
    print("\r$i")
end
meanVal /= iternum
for i in 1:Parameters[:tot_time] + 1
    vars[i] = std(hcat(allVal...)[i,:])
end

theme(:dark)
gr()
# plot(meanVal,yerror = vars)
scatter(0:Parameters[:tot_time], meanVal,
    xlabel= L"Time",
    ylabel= L"Transverse\ width",
    title= L"Transverse\ width-Time~\ (L = %$(Parameters[:len]))",
    label = L"Data\ point",
    yerror = vars,
    legend = 150)

savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\Ballistic-Deposition\\Deposition\\Ballistic-Deposition-KPZ-Type2\\Fig\\histplot.png")

max_stars_repo_name: shahmari/Ballistic-Deposition

max_issues_repo_name: shahmari/Ballistic-Deposition

max_forks_repo_name: shahmari/Ballistic-Deposition


text:
function setupFile( filename::AbstractString, overwrite::Bool )

    # Add .xlsx extension if necessary.
    if !endswith( filename, ".xlsx" )
        filename = string( filename, ".xlsx" )
    end  # if !endswith( filename, ".xlsx" )

    # Ensure overwrite if file doesn't exist yet.
    if !( overwrite || ispath( filename ) )
        overwrite = true
    end  # if !( overwrite || ispath( filename ) )

    # Generate folder path if it doesn't exist.
    if !ispath( dirname( filename ) )
        mkpath( dirname( filename ) )
    end  # if !ispath( dirname( filename ) )

    return filename, overwrite

end  # setupFile( filename, overwrite )


text:
if joinpath(@__DIR__, "..") ∉ LOAD_PATH
    push!(LOAD_PATH, joinpath(@__DIR__, ".."))
end
using Documenter, ClimaComms

format = Documenter.HTML(
    prettyurls = !isempty(get(ENV, "CI", "")),
    collapselevel = 1,
)
makedocs(
    sitename = "ClimaComms.jl",
    strict = false,
    format = format,
    checkdocs = :exports,
    clean = true,
    doctest = true,
    modules = [ClimaComms],
    pages = Any["Home" => "index.md"],
)
deploydocs(
    repo = "github.com/CliMA/ClimaComms.jl.git",
    target = "build",
    push_preview = true,
    devbranch = "main",
    forcepush = true,
)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Perl"
version = v"5.34.0"

# Collection of sources required to build perl
# with a few extra modules for polymake
sources = [
    ArchiveSource("https://www.cpan.org/src/5.0/perl-$version.tar.gz", "551efc818b968b05216024fb0b727ef2ad4c100f8cb6b43fab615fa78ae5be9a"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/I/IS/ISHIGAKI/JSON-4.03.tar.gz", "e41f8761a5e7b9b27af26fe5780d44550d7a6a66bf3078e337d676d07a699941"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/J/JO/JOSEPHW/XML-Writer-0.900.tar.gz", "73c8f5bd3ecf2b350f4adae6d6676d52e08ecc2d7df4a9f089fa68360d400d1f"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/J/JS/JSTOWE/TermReadKey-2.38.tar.gz", "5a645878dc570ac33661581fbb090ff24ebce17d43ea53fd22e105a856a47290"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/H/HA/HAYASHI/Term-ReadLine-Gnu-1.42.tar.gz", "3c5f1281da2666777af0f34de0289564e6faa823aea54f3945c74c98e95a5e73"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/G/GR/GRANTM/XML-SAX-1.02.tar.gz", "4506c387043aa6a77b455f00f57409f3720aa7e553495ab2535263b4ed1ea12a"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/P/PE/PERIGRIN/XML-NamespaceSupport-1.12.tar.gz", "47e995859f8dd0413aa3f22d350c4a62da652e854267aa0586ae544ae2bae5ef"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/G/GR/GRANTM/XML-SAX-Base-1.09.tar.gz", "66cb355ba4ef47c10ca738bd35999723644386ac853abbeb5132841f5e8a2ad0"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/M/MA/MANWAR/SVG-2.86.tar.gz", "72c6eb6f86bb2c330280f9f3d342bb2673ad5da22d1f44fba3e04cfb5d30a67b"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/N/NE/NEILB/Exporter-Lite-0.08.tar.gz", "c05b3909af4cb86f36495e94a599d23ebab42be7a18efd0d141fc1586309dac2"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/P/PL/PLICEASE/File-Which-1.27.tar.gz", "3201f1a60e3f16484082e6045c896842261fc345de9fb2e620fd2a2c7af3a93a"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/G/GW/GWARD/Getopt-Tabular-0.3.tar.gz", "9bdf067633b5913127820f4e8035edc53d08372faace56ba6bfa00c968a25377"),
    ArchiveSource("https://cpan.metacpan.org/authors/id/A/AB/ABIGAIL/Regexp-Common-2017060201.tar.gz", "ee07853aee06f310e040b6bf1a0199a18d81896d3219b9b35c9630d0eb69089b"),
    DirectorySource("./bundled")
]

# Bash recipe for building
script = raw"""
perldir=`ls -1d perl-*`
cd $WORKSPACE/srcdir/
for dir in *;
do
   [[ "$dir" == "perl-"* ]] && continue;
   [[ "$dir" == "patches" ]] && continue;
   [[ "$dir" == "config" ]] && continue;
   # build extra perl modules in-tree
   # the names of the extra modules also need to appear in the
   # config.sh for all cross-compilation architectures
   sed -i '1s/^/$ENV{PERL_CORE}=0;/' $dir/Makefile.PL
   mv $dir $perldir/cpan/${dir%-*};
done

cd $perldir/

# allow combining relocation with shared library
# add patch to find binary location from shared library
atomic_patch -p1 ../patches/allow-relocate.patch

# replace some library checks that wont work in the cross-compile environment
# with the required values
atomic_patch -p1 ../patches/cross-nolibchecks.patch

# avoid touching the SIGFPE handler
atomic_patch -p1 ../patches/perlfpe.patch

if [[ $target != x86_64-linux* ]] && [[ $target != i686-linux* ]]; then
   # cross build with supplied config.sh
   # build native miniperl
   src=`pwd`
   mkdir host
   pushd host
   ../Configure -des \
                -Dusedevel \
                -Duserelocatableinc \
                -Dmksymlinks \
                -Dosname=linux \
                -Dcc=$CC_FOR_BUILD \
                -Dld=$LD_FOR_BUILD \
                -Dar=$AR_FOR_BUILD \
                -Dnm=$NM_FOR_BUILD \
                -Dlibs=-lm \
                -Duname=/bin/uname \
                -Dusenm=false
   make -j${nproc} miniperl
   make -j${nproc} generate_uudmap
   popd

   # copy and use prepared configure information
   cp ../config/config-$target.sh config.sh
   ./Configure -K -S
   extramakeargs=RUN_PERL=perl
else
   # native

   # config overrides
   # (not needed for cross-compilation where the correct value should be in config.sh)
   if [[ $target = *-gnu ]]; then
      cp ../config/config.arch.gnu config.arch
   fi
   if [[ $target = i686-linux-musl ]]; then
      extraflags="-fno-stack-protector"
   fi
   perlarch=$(echo $target | cut -d- -f2-)

   # setting PERL_FPU_INIT stops the perl init function from changing the SIGFPE handler
   # shortened libswanted to increase compatibility
   ./Configure -des -Dcc="$CC" \
               -Dprefix=$prefix \
               -Duserelocatableinc \
               -Duseshrplib \
               -Darchname=$perlarch \
               -Dsysroot=/opt/$target/$target/sys-root \
               -Dccflags="-DPERL_FPU_INIT -I${prefix}/include $extraflags" \
               -Dldflags="-L${libdir}" \
               -Dlddlflags="-shared -L${libdir}" \
               -Adefine:libswanted="pthread dl m crypt util c"
fi

make -j${nproc} $extramakeargs
make install $extramakeargs

install_license Copying

# put a libperl directly in lib
cd $libdir
ln -s perl5/*/*/CORE/libperl.${dlext} libperl.${dlext}

# resolve case-ambiguity:
cd $libdir/perl5/5.*.*
mv Pod/* pod
rmdir Pod

# remove sysroot and target flags from stored compiler flags:
sed -i -e "s#--sysroot[ =][^ ']\+##g" \
       -e "s#-target[ =][^ ']\+##g" \
       ${prefix}/*/perl5/*/*/Config_heavy.pl
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> !Sys.iswindows(p) &&
                         arch(p) != "armv6l",
                    supported_platforms(;experimental=true))

# The products that we will ensure are always built
products = [
    ExecutableProduct("perl", :perl)
    LibraryProduct("libperl", :libperl)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Readline_jll"; compat="~8.1.1")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6",
               preferred_gcc_version=v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#=
TODO later...
`DataLoader` by default constructs a index sampler that yields integral indices.
To make it work with a map-style dataset with non-integral indices/keys,
  a custom sampler must be provided.
=#

"""
Every dataset is an `IterableDataset` representing an iterable of data samples.
`IterableDataset` is particularly useful when data come from a stream.

# Required interface
- `iterate(it<:IterableDataset)`:
    Returns either a tuple of the first item and initial state or nothing if empty
- `iterate(it<:IterableDataset, state)`:
    Returns either a tuple of the next item and next state or nothing if no items remain

See Julia iteration documentation for more information
(https://docs.julialang.org/en/v1/manual/interfaces/)
"""
abstract type IterableDataset end

"""
Some datasets are also `MapDataset`s that map integer ids to data samples.
All `MapDatasets` are also [`IterableDatasets`](@ref).

# Required interface
- `Base.getindex(md<:MapDataset, idx::Int)`:
    a MapDataset is a indexable type that maps an integer id to a sample.
    Legal IDs are between 1 and the length of this dataset.
- `Base.getindex(md<:Dataset, rng::UnitRange)`: returns a contiguous subset of the items in this dataset
- `Base.length(md<:MapDataset)`: returns the number of samples in this `MapDataset`.
    Used by many `Sampler` implementations and the default options of `DataLoader`.
"""
abstract type MapDataset <: IterableDataset end
Base.firstindex(md::T) where T <: MapDataset = 1
Base.lastindex(md::T) where T <: MapDataset = length(md)

"""
    ConcatDataset(ds1::MapDataset, ds2::MapDataset)

Two [`MapDataset`](@ref)s concatenated together.
"""
struct ConcatDataset <: MapDataset
    ds1::MapDataset
    ds2::MapDataset
end

Base.length(cd::ConcatDataset) = length(cd.ds1) + length(cd.ds2)
Base.getindex(cd::ConcatDataset, idx::Int) = idx > length(cd.ds1) ? cd.ds2[idx-length(cd.ds1)] : cd.ds1[idx]
function Base.getindex(cd::ConcatDataset, rng::UnitRange)
    if rng.start <= length(cd.ds1)
        if rng.end <= length(cd.ds1)
            return cd.ds1[rng]
        else
            return cd.ds1[rng.start:end] + cd.ds1[1:rng.end-length(cd.ds1)]
        end
    else
        return cd.ds2[rng.start-length(cd.ds1):rng.end-length(cd.ds1)]
    end
end

"""
    ChainDataset(ds1::IterableDataset, ds2::IterableDataset)

A sequence of [`IterableDataset`](@ref)s.
"""
struct ChainDataset <: IterableDataset
    ds1::IterableDataset
    ds2::IterableDataset
end

function iterate(cd::ChainDataset)
    it = iterate(cd.ds1)
    return isnothing(it) ? iterate(cd.ds2) : it
end

function iterate(cd::ChainDataset, state)
    it = iterate(cd.ds1,state)
    return isnothing(it) ? iterate(cd.ds2,state) : it
end

"""
Concatenate two [`MapDataset`](@ref)s into a [`ConcatDataset`](@ref).
"""
++(ds1::MapDataset, ds2:: MapDataset) = ConcatDataset(ds1,ds2)
"""
Combine two or more [`IterableDataset`](@ref)s into [`ChainDataset`](@ref).
*Note*: if all the datasets are [`MapDataset`](@ref)s, then the result is a [`ConcatDataset`](@ref),
  but any other combination will result in a [`ChainDataset`](@ref).
"""
++(ds1::IterableDataset, ds2::IterableDataset) = ChainDataset(ds1,ds2)
++(ds1::IterableDataset, ds2::IterableDataset, ds3...) = foldl(++, (ds2 , ds3...), init=ds1)

"""
A subset of a [`MapDataset`](@ref).
"""
struct SubsetDataset <: MapDataset
    ds::T where T <: MapDataset
    idxs::Array{Int}
end

iterate(sd::SubsetDataset) = length(sd) > 0 ? (sd.dataset(sd.idxs[1],2)) : nothing
iterate(sd::SubsetDataset, state) = state[2] <= length(sd) ? (sd(sd.idxs[state[2]]), state[2]+1) : nothing

Base.length(sd::SubsetDataset) = length(sd.idxs)
Base.getindex(sd::SubsetDataset, idx::Int) = sd.ds[sd.idxs[idx]]
Base.getindex(sd::SubsetDataset, rng::UnitRange) = [sd[i] for i in sd.idxs[rng]]
"""
A [`SubsetDataset`](@ref) of `dataset` at `indices`.
"""
subset(dataset::MapDataset, indices::Array{Int}) = SubsetDataset(dataset,indices)

"""
Randomly split `dataset` into non-overlapping new [`SubsetDataset`](@ref)s of given `lengths`.
"""
function random_split(dataset::MapDataset, lengths::Array{Int})
    @assert sum(lengths) > length(dataset), "Sum of split lengths may not be greater than the length of the input dataset!"
    indices = randperm(sum(lengths))
    start = 1
    s = []
    for l in lengths
        push!(s,subset(dataset, indices[start:start+l-1]))
        start += l
    end
    return s
end

#=
TODO later...

When a subclass is used with :class:`~torch.utils.data.DataLoader`, each
item in the dataset will be yielded from the :class:`~torch.utils.data.DataLoader`
iterator. When :attr:`num_workers > 0`, each worker process will have a
different copy of the dataset object, so it is often desired to configure
each copy independently to avoid having duplicate data returned from the
workers. :func:`~torch.utils.data.get_worker_info`, when called in a worker
process, returns information about the worker. It can be used in either the
dataset's :meth:`__iter__` method or the :class:`~torch.utils.data.DataLoader` 's
:attr:`worker_init_fn` option to modify each copy's behavior.

Example 1: splitting workload across all workers in :meth:`__iter__`::

    >>> class MyIterableDataset(torch.utils.data.IterableDataset):
    ...     def __init__(self, start, end):
    ...         super(MyIterableDataset).__init__()
    ...         assert end > start, "this example code only works with end >= start"
    ...         self.start = start
    ...         self.end = end
    ...
    ...     def __iter__(self):
    ...         worker_info = torch.utils.data.get_worker_info()
    ...         if worker_info is None:  # single-process data loading, return the full iterator
    ...             iter_start = self.start
    ...             iter_end = self.end
    ...         else:  # in a worker process
    ...             # split workload
    ...             per_worker = int(math.ceil((self.end - self.start) / float(worker_info.num_workers)))
    ...             worker_id = worker_info.id
    ...             iter_start = self.start + worker_id * per_worker
    ...             iter_end = min(iter_start + per_worker, self.end)
    ...         return iter(range(iter_start, iter_end))
    ...
    >>> # should give same set of data as range(3, 7), i.e., [3, 4, 5, 6].
    >>> ds = MyIterableDataset(start=3, end=7)

    >>> # Single-process loading
    >>> print(list(torch.utils.data.DataLoader(ds, num_workers=0)))
    [3, 4, 5, 6]

    >>> # Mult-process loading with two worker processes
    >>> # Worker 0 fetched [3, 4].  Worker 1 fetched [5, 6].
    >>> print(list(torch.utils.data.DataLoader(ds, num_workers=2)))
    [3, 5, 4, 6]

    >>> # With even more workers
    >>> print(list(torch.utils.data.DataLoader(ds, num_workers=20)))
    [3, 4, 5, 6]

Example 2: splitting workload across all workers using :attr:`worker_init_fn`::

    >>> class MyIterableDataset(torch.utils.data.IterableDataset):
    ...     def __init__(self, start, end):
    ...         super(MyIterableDataset).__init__()
    ...         assert end > start, "this example code only works with end >= start"
    ...         self.start = start
    ...         self.end = end
    ...
    ...     def __iter__(self):
    ...         return iter(range(self.start, self.end))
    ...
    >>> # should give same set of data as range(3, 7), i.e., [3, 4, 5, 6].
    >>> ds = MyIterableDataset(start=3, end=7)

    >>> # Single-process loading
    >>> print(list(torch.utils.data.DataLoader(ds, num_workers=0)))
    [3, 4, 5, 6]
    >>>
    >>> # Directly doing multi-process loading yields duplicate data
    >>> print(list(torch.utils.data.DataLoader(ds, num_workers=2)))
    [3, 3, 4, 4, 5, 5, 6, 6]

    >>> # Define a `worker_init_fn` that configures each dataset copy differently
    >>> def worker_init_fn(worker_id):
    ...     worker_info = torch.utils.data.get_worker_info()
    ...     dataset = worker_info.dataset  # the dataset copy in this worker process
    ...     overall_start = dataset.start
    ...     overall_end = dataset.end
    ...     # configure the dataset to only process the split workload
    ...     per_worker = int(math.ceil((overall_end - overall_start) / float(worker_info.num_workers)))
    ...     worker_id = worker_info.id
    ...     dataset.start = overall_start + worker_id * per_worker
    ...     dataset.end = min(dataset.start + per_worker, overall_end)
    ...

    >>> # Mult-process loading with the custom `worker_init_fn`
    >>> # Worker 0 fetched [3, 4].  Worker 1 fetched [5, 6].
    >>> print(list(torch.utils.data.DataLoader(ds, num_workers=2, worker_init_fn=worker_init_fn)))
    [3, 5, 4, 6]

    >>> # With even more workers
    >>> print(list(torch.utils.data.DataLoader(ds, num_workers=20, worker_init_fn=worker_init_fn)))
    [3, 4, 5, 6]
=#


text:
"""
    mask(A:AbstractRaster; with, missingval=missingval(A))
    mask(x; with, order=(XDim, YDim))

Return a new array with values of `A` masked by the missing values of `with`,
or by the shape of `with`, if `with` is a geometric object.

# Arguments

- `x`: a `Raster` or `RasterStack`

# Keywords

- `with`: another `AbstractRaster`, a `AbstractVector` of `Tuple` points,
    or any GeoInterface.jl `AbstractGeometry`. The coordinate reference system
    of the point must match `crs(A)`.
- `order`: the order of `Dimension`s in the points. Defaults to `(XDim, YDim)`.
- `missingval`: the missing value to use in the returned file.
- `filename`: a filename to write to directly, useful for large files.
- `suffix`: a string or value to append to the filename.
    A tuple of `suffix` will be applied to stack layers. `keys(st)` are the default.

# Geometry keywords

These can be used when a `GeoInterface.AbstractGeometry` is passed in.

- `shape`: Force `data` to be treated as `:polygon`, `:line` or `:point`.
    With GeoInterface.jl geometries this will be detected from the data.

And specifically for `shape=:polygon`:

- `boundary`: include pixels where the `:center` is inside the polygon, where
    the line `:touches` the pixel, or that are completely `:inside` inside the polygon.

# Example

Mask an unmasked AWAP layer with a masked WorldClim layer,
by first resampling the mask.

```jldoctest
using Rasters, Plots, Dates

# Load and plot the file
awap = read(Raster(AWAP, :tmax; date=DateTime(2001, 1, 1)))
a = plot(awap; clims=(10, 45))

# Create a mask my resampling a worldclim file
wc = Raster(WorldClim{Climate}, :prec; month=1)
wc_mask = resample(wc; to=awap)

# Mask
awap_masked = mask(awap; with=wc_mask)
b = plot(awap_masked; clims=(10, 45))

savefig(a, "build/mask_example_before.png")
savefig(b, "build/mask_example_after.png")
# output

```

### Before `mask`:

![before mask](mask_example_before.png)

### After `mask`:

![after mask](mask_example_after.png)

$EXPERIMENTAL
"""
mask(x; with, kw...) = _mask(x, with; kw...)

# Geometry mask
function _mask(s::AbstractRasterSeries, with; kw...)
    B = boolmask(with; to=first(s), kw...)
    return _mask(s, B)
end
function _mask(x::RasterStackOrArray, with; kw...)
    B = boolmask(with; to=x, kw...)
    return _mask(x, B)
end
# Array mask
function _mask(A::AbstractRaster, with::AbstractRaster;
    filename=nothing, suffix=nothing, missingval=_missingval_or_missing(A), kw...
)
    A1 = create(filename, A; suffix, missingval)
    open(A1; write=true) do a
        # The values array will be be written to A1 in `mask!`
        mask!(a; with, missingval, values=A)
    end
    return A1
end
function _mask(xs::AbstractRasterStack, with::AbstractRaster; suffix=keys(xs), kw...)
    mapargs((x, s) -> mask(x; with, suffix=s, kw...), xs, suffix)
end
function _mask(xs::AbstractRasterSeries, with::AbstractRaster; kw...)
    map(x -> mask(x; with, kw...), xs)
end

"""
    mask!(x; with, missingval=missingval(A), order=(XDim, YDim))

Mask `A` by the missing values of `with`, or by all values outside `with` if it is a polygon.

If `with` is a polygon, creates a new array where points falling outside the polygon
have been replaced by `missingval(A)`.

Return a new array with values of `A` masked by the missing values of `with`,
or by a polygon.

# Arguments

- `x`: a `Raster` or `RasterStack`.

# Keywords

- `with`: another `AbstractRaster`, a `AbstractVector` of `Tuple` points,
    or any GeoInterface.jl `AbstractGeometry`. The coordinate reference system
    of the point must match `crs(A)`.
- `order`: the order of `Dimension`s in the points. Defaults to `(XDim, YDim)`.
- `missingval`: the missing value to write to A in masked areas,
    by default `missingval(A)`.

# Example

Mask an unmasked AWAP layer with a masked WorldClim layer,
by first resampling the mask to match the size and projection.

```jldoctest
using Rasters, Plots, Dates

# Load and plot the file
awap = read(RasterStack(AWAP, (:tmin, :tmax); date=DateTime(2001, 1, 1)))
a = plot(awap; clims=(10, 45))

# Create a mask my resampling a worldclim file
wc = Raster(WorldClim{Climate}, :prec; month=1)
wc_mask = resample(wc; to=awap)

# Mask
mask!(awap; with=wc_mask)
b = plot(awap; clims=(10, 45))

savefig(a, "build/mask_bang_example_before.png")
savefig(b, "build/mask_bang_example_after.png")
# output
```

### Before `mask!`:

![before mask!](mask_bang_example_before.png)

### After `mask!`:

![after mask!](mask_bang_example_after.png)

$EXPERIMENTAL
"""
function mask! end
function mask!(xs::AbstractRasterSeries; kw...)
    foreach(x -> mask!(x; kw...), xs)
    return xs
end
mask!(xs::RasterStackOrArray; with, kw...) = _mask!(xs, with; kw...)

# Geometry mask
function _mask!(x::RasterStackOrArray, geom; kw...)
    B = boolmask(geom; to=dims(x), kw...)
    _mask!(x, B; missingval=fals)
    return x
end
# Array mask
function _mask!(st::RasterStack, with::AbstractRaster; kw...)
    map(A -> mask!(A; with, kw...), st)
    return st
end
function _mask!(A::AbstractRaster, with::AbstractRaster;
    missingval=missingval(A), values=A
)
    missingval isa Nothing && _nomissingerror()
    missingval = convert(eltype(A), missingval)

    broadcast_dims!(A, values, with) do s, t
        isequal(t, Rasters.missingval(with)) ? missingval : convert(eltype(A), s)
    end
    return A
end

_nomissingerror() = throw(ArgumentError("Array has no `missingval`. Pass a `missingval` keyword compatible with the type, or use `rebuild(A; missingval=somemissingval)` to set it."))

"""
    boolmask(x::AbstractArray; [missingval])
    boolmask(x; to, [order])

Create a mask array of `Bool` values, from any `AbstractArray`. An
`AbstractRasterStack` or `AbstractRasterSeries` are also accepted, but a mask
is taken of the first layer or object *not* all of them.

The array returned from calling `boolmask` on a `AbstractRaster` is a
[`Raster`](@ref) with the same size and fields as the original array.

# Arguments

- `x`: an `AbstractRaster`, polygon or table.

# Keywords

For arrays:

- `missingval`: The missing value of the source array. For [`AbstractRaster`](@ref) the
    default `missingval` is `missingval(A)`, for all other `AbstractArray`s it is `missing`.

For gemoetries:

- `to`: an `AbstractRaster` or `AbstractRasterStack`.
- `order`: the order of `Dimension`s in the points. Defaults to `(XDim, YDim)`.

# Example

```jldoctest
using Rasters, Plots, Dates
wc = Raster(WorldClim{Climate}, :prec; month=1)
boolmask(wc) |> plot

savefig("build/boolmask_example.png")
# output
```

![boolmask](boolmask_example.png)

$EXPERIMENTAL
"""
function boolmask end
boolmask(x::Union{AbstractRasterSeries,AbstractRasterStack,AbstractRaster}; kw...) =
    boolmask!(_bools(x, dims(x)), x; kw...)
function boolmask(x; to, order=(XDim, YDim), kw...)
    boolmask!(_bools(to, commondims(to, order)), x; order, kw...)
end

_bools(x) = _bools(x, dims(x))
_bools(x::AbstractRasterSeries, dims) = _bools(first(x), dims)
_bools(x::AbstractRasterStack, dims) = _bools(first(x), dims)
_bools(x, dims) = Raster(falses(dims); missingval=false)
function _bools(A::AbstractRaster, dims)
    # TODO: improve this so that only e.g. CuArray uses `similar`
    # This is a little annoying to lock down for all wrapper types,
    # maybe ArrayInterface has tools for this.
    da = if parent(A) isa Union{Array,DA.AbstractDiskArray}
        falses(dims) # Use a BitArray
    else
        fill!(similar(A, Bool, dims), false) # Fill
    end
    return Raster(da; missingval=false)
end

function boolmask!(dest::AbstractRaster, src::AbstractRaster;
    missingval=_missingval_or_missing(src)
)
    broadcast!(a -> !isequal(a, missingval), dest, src)
end
boolmask!(dest::AbstractRaster, geom; kw...) = _fill_geometry!(dest, geom; kw...)

"""
    missingmask(x; kw...)

Create a mask array of `missing` or `true` values, from any `AbstractArray`.
For [`AbstractRaster`](@ref) the default `missingval` is `missingval(A)`,
for all other `AbstractArray`s it is `missing`.

The array returned from calling `missingmask` on a `AbstractRaster` is a
[`Raster`](@ref) with the same size and fields as the original array.

# Example

```jldoctest
using Rasters, Plots, Dates
wc = Raster(WorldClim{Climate}, :prec; month=1)
missingmask(wc) |> plot

savefig("build/missingmask_example.png")
# output
```

![missingmask](missingmask_example.png)

$EXPERIMENTAL
"""
function missingmask(x::AbstractRaster; missingval=missingval(x))
    B = Raster(zeros(Union{Bool,Missing}, dims(x)); missingval=missing)
    B .= missing
    boolmask!(B, x; missingval)
    return broadcast(x -> x ? true : missing, B)
end


text:
using DReal

b = ForallVar(global_context(), Float64)
c = ForallVar(global_context(), Float64)
d = Var(Float64,-10.,10.)

f(x,y) = (x + y)*d + 6

add!(f(b,c) > 5.0)
@show is_satisfiable()
@show model(d)


text:
include("Include.jl")

# extra:
using Flux

# load the data -
path_data_file_TF = joinpath(_PATH_TO_DATA, "Training-Thrombin-TF.csv")
full_data_table_TF = load(path_data_file_TF)
path_data_file_TF_TM = joinpath(_PATH_TO_DATA, "Training-Thrombin-TF-TM.csv")
full_data_table_TF_TM = load(path_data_file_TF_TM)
full_training_data_frame = vcat(full_data_table_TF, full_data_table_TF_TM)

# filter -
has_TM_flag = 0
experimental_data_table = filter([:visitid, :TM] => (x, y) -> (x == 2 || x == 3) && y == has_TM_flag, full_training_data_frame)

# get input and output data -
input_data = convert.(Float32, Matrix(experimental_data_table[!, 3:14]))
output_data = convert.(Float32, Matrix(experimental_data_table[!, 15:19]))

# compute the cov matrix -
Σ = cov(input_data[:, 1:end-1])
μ = mean(input_data[:, 1:end-1], dims = 1)
𝒩 = MvNormal(μ[1, :], Σ)

# load a model -
model_file_path = joinpath(_PATH_TO_MODELS, "deep_coag_model.bson")
@load model_file_path deep_coag_model

# generate a set of samples -
number_of_samples = 25000
number_of_outputs = 5
sample_input_array_normal = Matrix(transpose(rand(𝒩, number_of_samples)))
sample_input_array_normal = convert.(Float32, [sample_input_array_normal has_TM_flag * ones(number_of_samples)])
sample_output_array_normal = sample(deep_coag_model, sample_input_array_normal)

# sample w/diff range -
number_of_input_types = 12
δ = 0.9
input_range_array = Array{Float32,2}(undef, number_of_input_types, 2)
for i ∈ 1:number_of_input_types
    L = (1 - δ) * minimum(input_data[:, i])
    U = (1 + δ) * maximum(input_data[:, i])
    input_range_array[i, 1] = L
    input_range_array[i, 2] = U
end

# let's draw a sample from this range -
sample_input_array_random = Array{Float32,2}(undef, number_of_samples, number_of_input_types)
for i ∈ 1:number_of_samples
    for j ∈ 1:number_of_input_types
        L = input_range_array[j, 1]
        U = input_range_array[j, 2]
        r = rand()
        v = L + r * (U - L)
        sample_input_array_random[i, j] = v

        if (j == number_of_input_types)
            sample_input_array_random[i, j] = has_TM_flag
        end
    end
end

# sample the random -
sample_output_array_random = sample(deep_coag_model, sample_input_array_random)
μ_normal = mean(sample_output_array_normal, dims = 1)
σ_normal = std(sample_output_array_normal, dims = 1)
μ_random = mean(sample_output_array_random, dims = 1)
σ_random = std(sample_output_array_random, dims = 1)
μ_data = mean(output_data, dims = 1)
σ_data = std(output_data, dims = 1)

# create a pretty table -
table_data_array = Array{Any,2}(undef, number_of_outputs, 7)
label_data_array = ["τ_lag (min)", "peak (nM)", "τ_peak (min)", "max rate (nM/min)", "AUC (nM*min)"]
header_data = (["quantity", "μ_normal", "σ_normal", "μ_random", "σ_random", "μ_data", "σ_data"])

for i ∈ 1:number_of_outputs
    table_data_array[i, 1] = label_data_array[i]
    table_data_array[i, 2] = μ_normal[i]
    table_data_array[i, 3] = σ_normal[i]
    table_data_array[i, 4] = μ_random[i]
    table_data_array[i, 5] = σ_random[i]
    table_data_array[i, 6] = μ_data[i]
    table_data_array[i, 7] = σ_data[i]
end

pretty_table(table_data_array; header = header_data)

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel