kernel_write.jl

include(joinpath(@__DIR__, "math_kernels.jl"))
include(joinpath(@__DIR__, "string_math_kernels.jl"))

# The kernel-generating function, analogous to fgen.
kgen(num::Num; constants::Vector{Num}=Num[], overwrite::Bool=false, splitting::Symbol=:default, affine_quadratic::Bool=true) = kgen(num, setdiff(pull_vars(num), constants), [:all], constants, overwrite, splitting, affine_quadratic)
kgen(num::Num, gradlist::Vector{Num}; constants::Vector{Num}=Num[], overwrite::Bool=false, splitting::Symbol=:default, affine_quadratic::Bool=true) = kgen(num, gradlist, [:all], constants, overwrite, splitting, affine_quadratic)
kgen(num::Num, raw_outputs::Vector{Symbol}; constants::Vector{Num}=Num[], overwrite::Bool=false, splitting::Symbol=:default, affine_quadratic::Bool=true) = kgen(num, setdiff(pull_vars(num), constants), raw_outputs, constants, overwrite, splitting, affine_quadratic)
kgen(num::Num, gradlist::Vector{Num}, raw_outputs::Vector{Symbol}; constants::Vector{Num}=Num[], overwrite::Bool=false, splitting::Symbol=:default, affine_quadratic::Bool=true) = kgen(num, gradlist, raw_outputs, constants, overwrite, splitting, affine_quadratic)
function kgen(num::Num, gradlist::Vector{Num}, raw_outputs::Vector{Symbol}, constants::Vector{Num}, overwrite::Bool, splitting::Symbol, affine_quadratic::Bool)
    # Create a hash of the expression and check if the function already exists
    expr_hash = string(hash(string(num)*string(gradlist)), base=62)
    if (overwrite==false) && (isfile(joinpath(@__DIR__, "storage", "f_"*expr_hash*".jl")))
        try func_name = eval(Meta.parse("f_"*expr_hash))
            return func_name
        catch
            include(joinpath(@__DIR__, "storage", "f_"*expr_hash*".jl"))
            func_name = eval(Meta.parse("f_"*expr_hash))
            @eval $(func_name)(CUDA.zeros(Float64, 2, 4+2*length($gradlist)), [CUDA.zeros(Float64, 2,3) for i = 1:length($gradlist)]...)
            return func_name
        end
    end
    
    # If we reach this comment, we're going to be creating/modifying the kernel file.
    # Create/open the file in write mode.
    file = open(joinpath(@__DIR__, "storage", "f_"*expr_hash*".jl"), "w")
    write(file, "# Generated at $(Dates.now())\n\n")
    write(file, "# Kernel(s) generated for the expression: $(string(num))\n\n")
    close(file)

    # Parse the list of requested outputs (Not currently useful, since all outputs are generated by default)
    func_outputs = Symbol[]
    for output in raw_outputs
        output == :cv     && push!(func_outputs, :cv)
        output == :cc     && push!(func_outputs, :cc)
        output == :lo     && push!(func_outputs, :lo)
        output == :hi     && push!(func_outputs, :hi)
        output == :MC     && push!(func_outputs, :cv, :cc, :lo, :hi)
        output == :cvgrad && push!(func_outputs, :cvgrad)
        output == :ccgrad && push!(func_outputs, :ccgrad)
        output == :grad   && push!(func_outputs, :cvgrad, :ccgrad)
        output == :all    && push!(func_outputs, :cv, :cc, :lo, :hi, :cvgrad, :ccgrad)
        if ~(output in [:cv, :cc, :lo, :hi, :MC, :cvgrad, :ccgrad, :grad, :all])
            error("Output list contains an invalid output symbol: :$output. Acceptable symbols
                   include [:cv, :cc, :lo, :hi, :MC, :cvgrad, :ccgrad, :grad, :all]")
        end
    end
    if isempty(func_outputs)
        error("No outputs specified.")
    end

    # Check the number of independent variables in the expression. Note that
    # even if gradlist has 32+ elements, if the expression itself has fewer
    # than 32 independent variables, it can be handled in a single kernel.
    # Note also that constants are still required as inputs, so we can't
    # exclude them from the list of participating variables
    indep_vars = get_name.(pull_vars(num))

    # Deal with quadratic functions differently, if `affine_quadratic` is true.
    # This uses SCIP's method of handling nonconvex quadratic terms, which is explained
    # in section 2.4.3.2 of:
    # Vigerske, S. and Gleixner, A. "SCIP: global optimization of mixed-integer 
    # nonlinear programs in a branch-and-cut framework". Optimization Methods 
    # and Software, 33:3, 563-593 (2018). DOI: 10.1080/10556788.2017.1335312
    # This method is also used in EAGO's `affine_relax_quadratic!` function.
    if affine_quadratic==true && is_quadratic(num) # NOTE: When switching to MOI variables, this will be easy to detect
        func_name = kgen_affine_quadratic(expr_hash, num, gradlist, func_outputs, constants)
        return func_name
    end

    # Now we need to split the input expression, if it's too complicated
    # to handle quickly in a single kernel. Complexity is determined by
    # the complexity() function, which (slightly conservatively) estimates
    # the number of lines the kernel will require. The two options for 
    # splitting based on complexity are:
    # Default: Split at 1500 or greater, unless the total remaining is under 2k
    #          (Also, split if >31 variables needed)
    # Conservative: Slower compilation, faster execution. Split at 10k or greater,
    #               unless total remaining is under 13k. (Also, split if >31 
    #               variables needed)

    # Begin by factoring and substituting the expression
    factored = perform_substitutions(factor(num, split_div=true))

    # Get the number of lines and variables associated with each factor
    n_lines = complexity(factored)
    n_vars = var_counts(factored)

    # Prepare input and output information
    kernel_nums = Int[]
    inputs = Vector{String}[]
    outputs = String[]

    # Pick the split points based on the "splitting" input
    if splitting==:low # Low, formerly none
        split_point = Inf
        max_size = Inf
    elseif splitting==:default # Default, formerly low
        split_point = 10000
        max_size = 12000
    elseif splitting==:high # Formerly default
        split_point = 1500
        max_size = 2000
    elseif splitting==:max # Extremely small
        split_point = 500
        max_size = 750
    else
        error("Splitting must be one of: {:low, :default, :high, :max}")
    end

    # Pull out sparsity information in the factorization
    sparsity = detect_sparsity(factored, gradlist)

    # Decide if the kernel needs to be split
    if (n_vars[end] < 31) && ((n_lines[end] <= max_size) || (findfirst(x -> x > split_point, n_lines)==length(n_lines)))
        # Complexity is fairly low; only a single kernel needed
        create_kernel!(expr_hash, 1, num, get_name.(gradlist), func_outputs, constants, factored, sparsity)
        push!(kernel_nums, 1)
        push!(inputs, string.(indep_vars))
        push!(outputs, "OUT")
    else
        # Complexity is not low enough; need multiple kernels
        complete = false
        kernel_count = 1
        # structure_list = String[] # Experimental
        while !complete
            # Determine which line to break at
            line_ID = findfirst(x -> x > split_point, n_lines)
            vars_ID = findfirst(x -> (x == 30) || (x == 31), n_vars)
            if isnothing(vars_ID)
                new_ID = line_ID
            elseif isnothing(line_ID)
                new_ID = vars_ID
            else
                new_ID =  min(line_ID, vars_ID)
            end

            # =============================================================================
            # =============================================================================
            # EXPERIMENTAL: 

            # We want to make the element at `new_ID` into a kernel, unless the structure
            # is exactly the same as a kernel we made previously.
            # new_term = extract(factored, new_ID)
            # new_structure, order = structure(new_term)
            # @show new_structure
            # # @show structure_list
            # if new_structure in structure_list
            #     # We've already made this structure. Identify which kernel that was and use that
            #     # one instead of making a new kernel
            #     kernel_ID = findfirst(x -> x==new_structure, structure_list)
            #     push!(kernel_nums, kernel_ID)
            #     push!(inputs, string.(get_name.(order)))
            #     push!(outputs, string(factored[new_ID].lhs))          
            #     println("Was in the list")  
            # else
            #     # Send the element at `new_ID` to create_kernel!()
            #     create_kernel!(expr_hash, kernel_count, new_term, get_name.(gradlist), func_outputs, constants)
            #     push!(structure_list, new_structure)
            #     push!(kernel_nums, kernel_count)
            #     push!(inputs, string.(get_name.(pull_vars(extract(factored, new_ID)))))
            #     push!(outputs, string(factored[new_ID].lhs))
            #     kernel_count += 1
            # end

            # =============================================================================
            # =============================================================================
            #### Start of alternative to experimental section

            # Send the element at `new_ID` to create_kernel!()
            create_kernel!(expr_hash, kernel_count, extract(factored, new_ID), get_name.(gradlist), func_outputs, constants, factored, sparsity)
            push!(kernel_nums, kernel_count)
            push!(inputs, string.(get_name.(pull_vars(extract(factored, new_ID)))))
            push!(outputs, string(factored[new_ID].lhs))
            kernel_count += 1

            #### End of alternative to experimental section

            # Eliminate this part of the factored list, since we've already calculated
            # it from this kernel
            factored[new_ID] = factored[new_ID].lhs ~ factored[new_ID].lhs

            # Re-calculate line and variable counts
            n_lines = complexity(factored)
            n_vars = var_counts(factored)

            # If the total number of lines (not including the final line) is below the max size 
            # and the number of variables is below 32, we can make the final kernel and be done
            if (n_vars[end] < 32) && (all(n_lines[1:end-1] .<= max_size))
                create_kernel!(expr_hash, kernel_count, extract(factored), get_name.(gradlist), func_outputs, constants, factored, sparsity)
                push!(kernel_nums, kernel_count)
                push!(inputs, string.(get_name.(pull_vars(extract(factored)))))
                push!(outputs, "OUT")
                complete = true
            end
        end
    end

    # Include all the kernels that were just generated
    include(joinpath(@__DIR__, "storage", "f_"*expr_hash*".jl"))

    # We can assume that the newly created kernels have sufficiently
    # high register usage that they all have a max number of 256 threads.
    # All that we need is to figure out the maximum number of blocks
    # on the user's machine, and we have some decent parameters with
    # which to call @cuda. This command checks the number of Streaming
    # Multiprocessors (SMs) in the user's GPU.
    blocks = Int32(CUDA.attribute(CUDA.device(), CUDA.DEVICE_ATTRIBUTE_MULTIPROCESSOR_COUNT))

    # Now we need to write the main CPU function that calls all of
    # the generated kernels and append it to the file.
    file = open(joinpath(@__DIR__, "storage", "f_"*expr_hash*".jl"), "a")
    write(file, outro(expr_hash, kernel_nums, inputs, outputs, blocks, get_name.(gradlist)))
    close(file)

    # Compile the function and kernels
    include(joinpath(@__DIR__, "storage", "f_"*expr_hash*".jl"))
    func_name = eval(Meta.parse("f_"*expr_hash))
    @eval $(func_name)(CUDA.zeros(Float64, 2, 4+2*length($gradlist)), [CUDA.zeros(Float64, 2, 3) for i = 1:length($gradlist)]...)
    return func_name
end

# This is a quick function to detect if a Num object is quadratic. This function
# should be unnecessary when using MOI variables, since the object will already
# indicate if the term is quadratic.
function is_quadratic(num::Num)
    expr = num.val

    # The function is quadratic if it's composed of one simple term like:
    # [const], [var], [coeff]*[var], [var]*[var], [coeff]*[var]*[var]
    # or if it's an `Add` type, where each term is one of the above.
    if typeof(expr) <: Real
        # If the expr has some Real-valued type (Float64, Int64, etc.), treat
        # it like it's quadratic
        nothing
    elseif base_term(expr)
        # Easy enough to treat this as quadratic, since the normal process doesn't make sense
        nothing
    elseif exprtype(expr)==MUL
        if length(expr.dict)==1
            # One element in the dict. Check that the term is a base-level variable and that
            # its exponent is 1 or 2
            for (key, val) in expr.dict
                if !base_term(key)
                    # Inner term is not a base-level variable, so it's not quadratic
                    return false
                end
                if !(val in [1, 2])
                    # Inner term's exponent is not 1 or 2, so it's not quadratic
                    return false
                end
            end
        elseif length(expr.dict)==2
            # Two elements in the dict. Check that both are base-level variables and that
            # their exponents are exactly 1
            for (key, val) in expr.dict
                if !base_term(key)
                    # At least one term isn't a base-level variable, so it's not quadratic
                    return false
                end
                if !isone(val)
                    # At least one term has a non-one exponent, so it's not quadratic
                    return false
                end
            end
        else
            # More than 2 elements in the dict? Not quadratic.
            return false 
        end
    elseif exprtype(expr)==ADD
        # The coefficient is irrelevant; we only need to explore the dictionary.
        # The keys are important, but the values are irrelevant.
        for key in keys(expr.dict)
            if base_term(key)
                # SYM is fine
                nothing
            elseif exprtype(key)==MUL
                # It must be exactly [var]*[var], or it's not quadratic
                if length(key.dict) != 2
                    return false
                end
                for (subkey, subval) in key.dict
                    if !base_term(subkey)
                        # Both terms must be base-level variables
                        return false
                    end
                    if !isone(subval)
                        # Both terms must have exponents of 1
                        return false
                    end
                end
            elseif exprtype(key)==POW
                # It must be exactly [var]^2, or it's not quadratic
                if !base_term(key.base)
                    return false
                end
                if key.exp != 2
                    return false
                end
            else
                # If it's not one of the above, it's not quadratic
                return false
            end
        end
    else
        return false # If it's not one of the above, it's not quadratic
    end
    return true
end

# A special version of kgen that only applies to quadratic functions. Instead of
# doing McCormick relaxations, this returns either affine bounds or secant line
# bounds, depending on where on the quadratic function the point of interest is.
function kgen_affine_quadratic(expr_hash::String, num::Num, gradlist::Vector{Num}, func_outputs::Vector{Symbol}, constants::Vector{Num})
    # Since it's quadratic, we can construct the kernel according to
    # `affine_relax_quadratic!` in EAGO.

    # Extract the variables that are participating in this expression
    expr = num.val
    vars = get_name.(pull_vars(expr))
    varlist = string.(get_name.(gradlist))

    # Open the file again
    file = open(joinpath(@__DIR__, "storage", "f_"*expr_hash*".jl"), "a")

    # Put in the preamble.
    if isempty(vars)
        write(file, preamble_string(expr_hash, ["OUT";], 1, 1, length(gradlist)))
    else
        write(file, preamble_string(expr_hash, ["OUT"; string.(vars)], 1, 1, length(gradlist)))
    end


    # Depending on the format of the expression, compose the kernel differently
    if typeof(expr) <: Real
        # It's something like a Float64. We initialize with the value of the expr
        write(file, SCMC_quadaff_initialize(expr))
        # And that's pretty much it... Wrap up afterwards as normal

    elseif base_term(expr)
        # If it's a SYM type, we initialize with zero
        write(file, SCMC_quadaff_initialize(0.0))

        # We can now add to the temporary subgradient storages
        # for this variable with a coefficient of 1.0
        write(file, SCMC_quadaff_affine(string(get_name(expr)), 1.0, varlist))

    elseif exprtype(expr)==MUL
        # If it's a MUL type, we initialize with zero
        write(file, SCMC_quadaff_initialize(0.0))

        # For this to be quadratic, there's either 1 term with an exponent
        # of 1 or 2, or there's 2 terms with exponents of 1. If it's 1 term
        # with an exponent of 1, it's essentially affine
        if length(expr.dict)==1
            for (key, val) in expr.dict
                if isone(val)
                    write(file, SCMC_quadaff_affine(string(get_name(key)), expr.coeff, varlist))
                else # Must be something squared
                    write(file, SCMC_quadaff_squared(string(get_name(key)), expr.coeff, varlist))
                end
            end
        else # There must be two elements in the dictionary
            binary_vars = string.(get_name.(keys(expr.dict)))
            binary_vars = binary_vars[sort_vars(binary_vars)]
            write(file, SCMC_quadaff_binary(binary_vars..., expr.coeff, varlist))
        end

    elseif exprtype(expr)==ADD
        # Final option. If it's ADD type, we initialize with the addition coefficient
        write(file, SCMC_quadaff_initialize(expr.coeff))

        # Now we must go through every term in ADD's dictionary and add it as appropriate
        for (key, val) in expr.dict
            if base_term(key)
                write(file, SCMC_quadaff_affine(string(get_name(key)), val, varlist))
            elseif exprtype(key)==MUL
                binary_vars = string.(get_name.(keys(key.dict)))
                binary_vars = binary_vars[sort_vars(binary_vars)]
                write(file, SCMC_quadaff_binary(binary_vars..., val, varlist))
            elseif exprtype(key)==POW
                write(file, SCMC_quadaff_squared(string(get_name(key.base)), val, varlist))
            else
                error("How did you get here? Why is $key in a quadratic term?")
            end
        end
    else
        error("How did you get here? Somehow this is marked as quadratic: $num")
    end

    # Wrap up the kernel by calculating bounds and relaxations
    # (NOTE: Technically, calculating the intercepts of the relaxations,
    #  which is what the kernel is doing normally, is more useful. ParBB
    #  will need to convert from relaxation values to intercepts anyway.
    #  EAGO already does this and bypasses the need to calculate relaxations.
    #  But, for compatibility with McCormick-style relaxations in ParBB,
    #  it's easier to simply calculate what ParBB is expecting.)
    if isempty(varlist)
        write(file, postamble_quadaff(String[], String[]))
    elseif isempty(vars)
        write(file, postamble_quadaff(String[], varlist))
    else
        write(file, postamble_quadaff(string.(vars), varlist))
    end
    close(file)

    # Include this kernel so SCMC knows what it is
    include(joinpath(@__DIR__, "storage", "f_"*expr_hash*".jl"))

    # Add onto the file the "main" CPU function that calls the kernel
    blocks = Int32(CUDA.attribute(CUDA.device(), CUDA.DEVICE_ATTRIBUTE_MULTIPROCESSOR_COUNT))
    file = open(joinpath(@__DIR__, "storage", "f_"*expr_hash*".jl"), "a")
    if isempty(gradlist)
        write(file, outro(expr_hash, [1], [String[]], ["OUT"], blocks, Symbol[]))
    elseif isempty(vars)
        write(file, outro(expr_hash, [1], [String[]], ["OUT"], blocks, get_name.(gradlist)))
    else
        write(file, outro(expr_hash, [1], [string.(vars)], ["OUT"], blocks, get_name.(gradlist)))
    end
    close(file)

    # Include the file again to get the final kernel
    include(joinpath(@__DIR__, "storage", "f_"*expr_hash*".jl"))
    func_name = eval(Meta.parse("f_"*expr_hash))
    @eval $(func_name)(CUDA.zeros(Float64, 2, 4+2*length($gradlist)), [CUDA.zeros(Float64, 2,4) for i = 1:length($gradlist)]...)
    return func_name
end


# Note: Currently doesn't consider specific outputs or constants.
# Outputs can be modified by passing that through to the string
# kernels, which will then replace, e.g., "OUT" with "#OUT".
# Not sure yet how I want to handle constants.
# This function takes information about the file name, kernel ID, and
# the expression that a SINGLE kernel is being created for, and creates
# that kernel in the specified file.
create_kernel!(expr_hash::String, kernel_ID::Int, num::Num, gradlist::Vector{Symbol}, func_outputs::Vector{Symbol}, constants::Vector{Num}, orig_factored::Vector{Equation}, orig_sparsity::Vector{Vector{Int}}) = create_kernel!(expr_hash, kernel_ID, num.val, gradlist, func_outputs, constants, orig_factored, orig_sparsity)
function create_kernel!(expr_hash::String, kernel_ID::Int, num::BasicSymbolic{Real}, gradlist::Vector{Symbol}, func_outputs::Vector{Symbol}, constants::Vector{Num}, orig_factored::Vector{Equation}, orig_sparsity::Vector{Vector{Int}})
    # This function will create a kernel for `num`, with the name:
    # "f_" * expr_hash * "_$n". This name will be pushed to `kernels`,
    # and a vector of the required inputs variables will be pushed to
    # `inputs`. Other inputs are needed to determine portions of the
    # kernel itself.

    # Start by factorizing the input expression
    factorized = factor(num, split_div=true)

    # Perform substitutions if possible
    factorized = perform_substitutions(factorized)

    # Collect all the LHS terms and participating variables
    LHS = string.(getfield.(factorized, :lhs))
    vars = get_name.(pull_vars(num))

    # Put the factorized expression into directed acyclic graph form
    edgelist, varids = eqn_edges(factorized) #varids includes all aux variables also
    g = SimpleDiGraph(edgelist)

    # Perform a topological sort to get the order in which we should
    # perform calculations (i.e., the final entry in "varorder" is the
    # full original expression)
    varorder = varids[topological_sort(g)]

    # Now, we need sparsity information for all the variables. We can pull sparsity information
    # normally, unless it's a temporary variable, in which case we have to refer to the original
    # factorization and sparsity.
    string_gradlist = string.(gradlist)
    sparsity = Vector{Vector{Int}}(undef, length(varorder))
    for i in eachindex(varorder)
        if varorder[i] in string_gradlist
            # Mark sparsity if the variable is already in gradlist
            sparsity[i] = [findfirst(==(varorder[i]), string_gradlist)]
        else
            # Find out what index we're on
            idx = findfirst(x -> isequal(string(x.lhs), varorder[i]), factorized)

            if isnothing(idx)
                sparsity[i] = orig_sparsity[findfirst(x -> isequal(string(x.lhs), varorder[i]), orig_factored)]
            else
                # Extract all the variables for this index
                vars = pull_vars(extract(factorized, idx))

                # For each variable in the expanded expression, add in sparsity information
                curr_sparsity = Int[]
                for var in vars
                    ID = findfirst(==(string(get_name(var))), string_gradlist)
                    if isnothing(ID)
                        # If we didn't find the variable, we need to scan the original factorization,
                        # and then pull sparsity info from the original sparsity list
                        ID = findfirst(x -> isequal(string(x.lhs), string(var)), orig_factored)
                        push!(curr_sparsity, orig_sparsity[ID]...)
                    else
                        # If we do find the variable, we can add this variable directly into the sparsity
                        push!(curr_sparsity, ID)
                    end
                end

                # Add a sorted, unique list to the sparsity tracker
                sparsity[i] = sort(unique(curr_sparsity))
            end
        end
    end

    # Check if we need temporary variables at all. We don't need
    # temporary variables if we only have addition, or if we have
    # the addition of single-rule terms, since we can just keep adding
    # new information to the existing output space. E.g.:
    #   x +   y +     z : No temporary variables needed
    #  3x + x*y +   z^2 : No temporary variables needed
    # x*y + y*z + x*y*z : Temporary variable needed because x*y*z is two rules
    # need_temps = depth_check(factorized)
    # NOTE: Alternatively, saving data to global GPU memory will definitely
    #       be slower than saving it to a temporary variable and copying it
    #       to global memory only at the end. Though, that would mean more
    #       local storage for each thread, which would limit the number of
    #       threads per SM that could be used (and affect occupancy). So,
    #       in some cases it may be better to store data in temporary variables,
    #       and in other cases it might be better to store directly to the
    #       final output location. This may require some testing, and then
    #       perhaps a flag that overrides "need_temps" and the subtraction
    #       of 1 from the temp count later on.
    #  (Disabling entirely for now)

    # Glossary:
    # varids: ALL variables including base variables and aux variables
    # varorder: A topologically sorted list of varids
    # vars: Variables that participate in the original expression 
    #       (i.e., NOT including ones produced through factorization)

    # g.fadjlist[i]: Contains indices of varids that depend on varids[i]
    # g.badjlist[i]: Contains indices of varids that are needed to compute varids[i]

    # Calculate the number of temporary variables needed. 
    temp_endlist = []
    maxtemp = 0
    # if need_temps #(Skip for now)
        for i in eachindex(varorder) # Loop through every participating variable
            if (varorder[i] in string.(get_name.(vars)))
                # Skip the variable if it's an input
                continue
            end
            # Find which index varorder[i] is in `varids`
            ID = findfirst(x -> occursin(varorder[i], x), varids)
            tempID = 0

            # If we are not already keeping track of temporary variables,
            # this becomes the first one
            if isempty(temp_endlist)
                # Keep track of what varids depend on this temporary variable
                push!(temp_endlist, copy(g.fadjlist[ID]))
                tempID = 1
            else
                # Check if this variable's expression uses addition. If so,
                # check if either of the RHS variables appeared earlier
                # in varorder.
                factorized_ID = findfirst(x -> isequal(string(x.lhs), varorder[i]), factorized)
                if exprtype(factorized[factorized_ID].rhs) == ADD
                    # Check through the temp_endlist to see if any temporary variables
                    # ONLY point to this ID (i.e., they aren't used elsewhere). If so,
                    # we can re-use that temporary variable and overwrite the results
                    # with addition.
                    for j in eachindex(temp_endlist)
                        if (length(temp_endlist[j])==1) && (temp_endlist[j][1]==ID)
                            temp_endlist[j] = copy(g.fadjlist[ID])
                            tempID = j
                            break
                        end
                    end
                end

                # Check if there are any temporary variables we can override
                for j in eachindex(temp_endlist)
                    if isempty(temp_endlist[j]) # Then we can override this one
                        temp_endlist[j] = copy(g.fadjlist[ID])
                        tempID = j
                        break
                    end
                end
                if tempID==0 # Then we haven't found one we can override
                    push!(temp_endlist, copy(g.fadjlist[ID]))
                    tempID = length(temp_endlist)
                end
            end

            # Now that we're done with this variable, look over other
            # temporary variables to see if they're no longer needed
            for j in eachindex(temp_endlist)
                if ID in temp_endlist[j]
                    filter!(x -> x!=ID, temp_endlist[j])
                end
            end
            if tempID > maxtemp
                maxtemp = tempID
            end
        end
    # end # Skipping outer loop for now

    # [Deprecating, using temporary variables to decrease global memory accesses]
    # # We have one more temporary variable than we need, since at
    # # least one result could have been stored in the final output
    # # of this kernel.
    # maxtemp -= 1

    # At this point, we should be ready to write the kernel. Open
    # the file in "append" mode so that other info that was written
    # won't be impacted.
    file = open(joinpath(@__DIR__, "storage", "f_"*expr_hash*".jl"), "a")

    # Put in the preamble.
    write(file, preamble_string(expr_hash, ["OUT"; string.(get_name.(vars))], kernel_ID, maxtemp, length(gradlist)))

    # Loop through the topological list to add calculations in order
    temp_endlist = []
    outvar = ""
    name_tracker = copy(varids)
    for i in eachindex(varorder) # Order in which variables are calculated
        # Skip calculation if the variable is one of the inputs
        if (varorder[i] in string.(get_name.(vars)))
            continue
        end

        # Determine the corresponding ID of the variable in varids
        ID = findfirst(x -> occursin(varorder[i], x), varids)

        # Get the inputs for this operation by checking the name
        # tracker (there might be a better way to do this... ah well)
        factorized_ID = findfirst(x -> isequal(string(x.lhs), varorder[i]), factorized)
        participants = get_name.(pull_vars(factorized[factorized_ID].rhs))
        inputs = []
        for p in string.(participants)
            # Find the corresponding element in varids
            varids_ID = findfirst(x -> isequal(x, p), varids)
            # Push the name_tracker name to the input list
            push!(inputs, name_tracker[varids_ID])
        end

        # [Deprecating; I'll use temporary variables the whole way and then set
        # the output variable at the end for final copying]
        # # If this is the final variable, it'll be called "OUT". No need
        # # for temp variables
        # if i==length(varorder)
        #     name_tracker[ID] = "OUT"
        # else

        # Determine which tempID to use/override. temp_endlist keeps
        # track of where variables will be used in the future (stored
        # as g.fadjlist), with elements removed as they are used. If
        # there is an empty row in temp_endlist, we can re-use that
        # tempID. If there isn't an empty row, we add a new row.
        tempID = 0

        # If we are not already keeping track of temporary variables,
        # this becomes the first one
        if isempty(temp_endlist)
            # Keep track of what varids depend on this temporary variable
            push!(temp_endlist, copy(g.fadjlist[ID]))
            tempID = 1
        else
            # Check if this variable's expression uses addition. If so,
            # check if we can reuse a temporary variable.
            factorized_ID = findfirst(x -> isequal(string(x.lhs), varorder[i]), factorized)
            if exprtype(factorized[factorized_ID].rhs) == ADD
                # Check through the temp_endlist to see if any temporary variables
                # ONLY point to this ID (i.e., they aren't used elsewhere). If so,
                # we can re-use that temporary variable and overwrite the results
                # with addition.
                for j in eachindex(temp_endlist)
                    if (length(temp_endlist[j])==1) && (temp_endlist[j][1]==ID)
                        temp_endlist[j] = copy(g.fadjlist[ID])
                        tempID = j
                        break
                    end
                end
            end

            # Check if there are any temporary variables we can override
            for j in eachindex(temp_endlist)
                if isempty(temp_endlist[j]) # Then we can override this one
                    temp_endlist[j] = copy(g.fadjlist[ID])
                    tempID = j
                    break
                end
            end
            if tempID==0 # Then we haven't found one we can override
                push!(temp_endlist, copy(g.fadjlist[ID]))
                tempID = length(temp_endlist)
            end
        end

        # When we refer to this variable in the future, we need to know what tempID
        # the variable is using
        name_tracker[ID] = "temp$(tempID)"

        # Now we can append this temporary variable to the list of inputs
        # for the correct operation
        inputs = [name_tracker[ID]; inputs]

        # Now we can pass the equation's RHS and the inputs to call the correct
        # writer function
        write_operation(file, factorized[factorized_ID].rhs, inputs, string.(gradlist), sparsity[i])

        # Now that we're done with this variable, eliminate this variable
        # from the lists of temporary variables' requirements
        for j in eachindex(temp_endlist)
            if ID in temp_endlist[j]
                filter!(x -> x!=ID, temp_endlist[j])
            end
        end
        
        # Keep track of the name of the output variable for setting the output
        if i==length(varorder)
            outvar = name_tracker[ID]
        end
    end

    # Now that the function is complete, we can close out the function
    write(file, postamble(outvar))
    close(file)

    # And here's the function, if we want to print it out...
    # @show joinpath(@__DIR__, "storage", "f_"*expr_hash*".jl")
    return nothing
end

# This function detects subexpressions in the factored list of equations
# that have known/improved relaxations, and replaces them with their
# improved forms. All LHS forms are replaced by the RHS equivalent.
# Some reformulations must be checked before others, since structures
# may change and hide possible reformulations. E.g., if we have the
# expression 2^(log(x*y)), the ideal reformulation is: (x*y)^log(2).
# But, based on the ordering of factors, "log(x*y)" will appear earlier
# than "2^(log(x*y)), and "log(x*y)" reformulates to "log(x) + log(y)".
# This reformulation would disable the ideal reformulation above,
# since we'd then have 2^(log(x) + log(y)). Larger-scale reformulations
# like this example are marked by "HIGH PRIORITY" and are searched for
# first, before any other reformulations.
# Included forms:
# 1) exp(x1)*exp(x2) = exp(x1+x2) [EAGO paper/Misener2014]
# 2) log(x1^a) = a*log(x1) [EAGO paper/Misener2014]
# 3) log(a^x1) = x1*log(a) [EAGO paper]
# 4) log(x1*x2) = log(x1) + log(x2) [EAGO paper]
# 5) (x1^a)^b = x1^(a*b) [EAGO paper]
# 6) a^(log(x1)) = x1^(log(a)) [EAGO paper] (HIGH PRIORITY)
# 7) log(inv(x1)) = -log(x1) [EAGO paper]
# 8) CONST1*CONST2*x1 = (CONST1*CONST2)*x1
# 9) 1 / (1 + exp(-x)) = Sigmoid(x)
# 10) sin(x) = cos(x - pi/2)
#
# Forms that aren't relevant yet:
# 1) (a^x1)^b = (a^b)^x1 [EAGO paper] (Can't do powers besides integers)
function perform_substitutions(old_factored::Vector{Equation})
    factored = deepcopy(old_factored)

    # Register any terms we want to substitute
    @eval @register_symbolic SCMC_sigmoid(x)

    scan_flag = true
    while scan_flag
        scan_flag = false

        # Higher priority reformulations
        for index0 in eachindex(factored)
            # 6) a^(log(x1)) = x1^(log(a)) [EAGO paper]
            if exprtype(factored[index0].rhs)==POW
                # We only apply this rule if the base of the power is Real-valued
                if typeof(factored[index0].rhs.base) <: Real
                    index1 = findfirst(x -> isequal(x.lhs, factored[index0].rhs.exp), factored)
                    if !isnothing(index1) && exprtype(factored[index1].rhs)==TERM
                        if factored[index1].rhs.f==log
                            # We also only want to apply this rule if the argument of 
                            # log() is a variable
                            if !(typeof(arguments(factored[index1].rhs)[]) <: Real)
                                # We have:
                                # aux1 = log(x1)
                                # aux2 = CONST^aux1
                                # 
                                # We will replace with:
                                # aux1 = log(x1)
                                # aux2 = x1^[log(CONST)]
                                scan_flag = true

                                # Identify the log variable and base constant value
                                log_var = arguments(factored[index1].rhs)[]
                                base_val = factored[index0].rhs.base

                                # Modify index0
                                ID_1 = findfirst(x -> isequal(x.rhs, log_var^log(base_val)), factored)
                                if isnothing(ID_1)
                                    @eval $factored[$index0] = $factored[$index0].lhs ~ $log_var^log($base_val)
                                else
                                    for i in eachindex(factored)
                                        @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $factored[$ID_1].lhs))
                                    end
                                    deleteat!(factored, index0)
                                end
                                break
                            end
                        end
                    end
                end
            end
        end

        # Lower priority reformulations
        if scan_flag==false
            for index0 in eachindex(factored)
                # 1) exp(x1)*exp(x2) = exp(x1+x2) [EAGO paper/Misener2014]
                if exprtype(factored[index0].rhs)==MUL
                    # We need to check all the multiplication arguments. We create
                    # a list of all the args using a recursive function.
                    args = pull_mult(factored, index0)
    
                    # Count the number of args that are "exp"
                    exp_count = 0
                    exp_args = []
                    if length(args) > 1
                        for arg in args
                            if !(typeof(arg) <: Real) && exprtype(arg)==TERM
                                if arg.f==exp
                                    exp_count += 1
                                    push!(exp_args, arguments(arg)[])
                                end
                            end
                        end
                    end
    
                    # If we have more than 1 "exp", we'll apply this rule. 
                    # Otherwise, we don't do anything.
                    if exp_count > 1
                        # We're going to create one big term and factor it, since
                        # the arrangement of multiplications in the original factorization
                        # might not line up with the changes
                        scan_flag = true
                        # Create the initial exp()
                        new_expr = exp(sum(exp_args))
                        # Multiply by all the non-exp terms
                        for arg in args
                            if exprtype(arg)!=TERM || arg.f!=exp
                                new_expr *= arg
                            end
                        end
                        # Create a factorization of this new expr
                        new_factorization = factor(new_expr, split_div=true)
                        # Scan through the new factorization to see if we can merge elements
                        # with the original factored list
                        done = false
                        last_match = nothing
                        while !done
                            done = true
                            for i in eachindex(new_factorization)
                                ID_1 = findfirst(x -> isequal(x.rhs, new_factorization[i].rhs), factored)
                                if !isnothing(ID_1)
                                    # Match was found in the main factorization. Remove this element
                                    # from the new factorization and replace all references to it with
                                    # the auxiliary variable of ID_1
                                    for j in eachindex(new_factorization)
                                        new_factorization[j] = new_factorization[j].lhs ~ substitute(new_factorization[j].rhs, Dict(new_factorization[i].lhs => factored[ID_1].lhs))
                                    end
                                    last_match = factored[ID_1].lhs
                                    deleteat!(new_factorization, i)
                                    done = false
                                    break
                                end
                            end
                        end
                        # Add remaining elements of the factorization to the location of index0
                        if length(new_factorization) >= 1
                            for j = length(new_factorization)-1:-1:1
                                insert!(factored, index0, new_factorization[j])
                            end
                            @eval $factored[$index0+length($new_factorization)-1] = $factored[$index0+length($new_factorization)-1].lhs ~ $new_factorization[end].rhs
                        else
                            # Or, if no elements exist, replace all instances of this new auxiliary variable with the last known match
                            for i in eachindex(factored)
                                @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $last_match))
                            end
                            deleteat!(factored, index0)
                        end
                        break
                    end
    
                end
    
                # 2) log(x1^a) = a*log(x1) [EAGO paper/Misener2014]
                # Equivalent in this code to:
                # 3) log(a^x1) = x1*log(a) [EAGO paper]
                if exprtype(factored[index0].rhs)==TERM
                    if factored[index0].rhs.f==log
                        index1 = findfirst(x -> isequal(x.lhs, arguments(factored[index0].rhs)[]), factored)
                        if !isnothing(index1) && exprtype(factored[index1].rhs)==POW
                            # If the argument of log() is POW type, we can separate based 
                            # on the `base` and `exp`
                            # start: 
                            # aux1 = log(x^a)
                            # 
                            # Convert to:
                            # aux2 = log(x)
                            # aux1 = a*aux2
                            scan_flag = true
                            factor_base = factored[index1].rhs.base
                            factor_exp = factored[index1].rhs.exp
                            
                            newsym = gensym(:aux)
                            newsym = Symbol(string(newsym)[3:5] * string(newsym)[7:end])
                            newvar = genvar(newsym)
                            if typeof(factor_base) <: Real
                                # Check whether the new term exists, and then either add it or only refer to that term
                                ID_1 = findfirst(x -> isequal(x.rhs, factor_exp*log(factor_base)), factored)
                                if isnothing(ID_1)
                                    @eval $factored[$index0] = $factored[$index0].lhs ~ $factor_exp*log($factor_base)
                                else
                                    for i in eachindex(factored)
                                        @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $factored[$ID_1].lhs))
                                    end
                                    deleteat!(factored, index0)
                                end
                            else
                                # Check whether terms already exist, and then either add them or only refer to them
                                ID_1 = findfirst(x -> isequal(x.rhs, log(factor_base)), factored)
                                if isnothing(ID_1)
                                    insert!(factored, index0, Equation(Symbolics.value(newvar), log(factor_base)))
                                    @eval $factored[$index0+1] = $factored[$index0+1].lhs ~ $factor_exp*$factored[$index0].lhs
                                else
                                    # The log term exists... does factor_exp*log_term also exist?
                                    ID_2 = findfirst(x -> isequal(x.rhs, factor_exp*factored[ID_1].lhs), factored)
                                    if isnothing(ID_2)
                                        @eval $factored[$index0] = $factored[$index0].lhs ~ $factor_exp*$factored[$ID_1].lhs
                                    else
                                        # Both the inner and outer terms already exist. All terms that refer to factored[index0]
                                        # should instead refer to the term that already exists
                                        for i in eachindex(factored)
                                            @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $factored[$ID_2].lhs))
                                        end
                                        deleteat!(factored, index0)
                                    end
                                end
                            end
                            break
                        end
                    end
                end
    
                # 4) log(x1*x2) = log(x1) + log(x2) [EAGO paper]
                if exprtype(factored[index0].rhs)==TERM
                    if factored[index0].rhs.f==log
                        index1 = findfirst(x -> isequal(x.lhs, arguments(factored[index0].rhs)[]), factored)
                        if !isnothing(index1) && exprtype(factored[index1].rhs)==MUL
                            # If the argument of log() is MUL type, we can create two new
                            # log() terms and add them together
                            # start: 
                            # aux1 = log(x*y)
                            # 
                            # Convert to:
                            # aux2 = log(x)
                            # aux3 = log(y)
                            # aux1 = aux2+aux3
                            #
                            # Or:
                            # aux1 = log(CONST*x)
                            # to:
                            # aux2 = log(x)
                            # aux1 = log(CONST) + aux2
                            scan_flag = true
                            args = [keys(factored[index1].rhs.dict)...]
                            coeff = factored[index1].rhs.coeff

                            if length(args)==1
                                # If there's only one arg, it's the CONST*x case
                                # Check whether log(x) already exists
                                ID_1 = findfirst(x -> isequal(x.rhs, log(args[1])), factored)
                                if isnothing(ID_1)
                                    # If we don't already have log(x), create a new auxiliary variable
                                    # and update the old expression
                                    newsym = gensym(:aux)
                                    newsym = Symbol(string(newsym)[3:5] * string(newsym)[7:end])
                                    newvar = genvar(newsym)
                                    insert!(factored, index0, Equation(Symbolics.value(newvar), log(args[1])))
                                    @eval $factored[$index0+1] = $factored[$index0+1].lhs ~ log($coeff) + $factored[$index0].lhs
                                else
                                    # If log(x) already exists, check if log(COEFF) + log(x) exists
                                    ID_2 = findfirst(x -> isequal(x.rhs, log(coeff) + factored[ID_1].lhs), factored)
                                    if isnothing(ID_2)
                                        @eval $factored[$index0] = $factored[$index0].lhs ~ log($coeff) + $factored[$ID_1].lhs
                                    else
                                        for i in eachindex(factored)
                                            @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $factored[$ID_2].lhs))
                                        end
                                        deleteat!(factored, index0)
                                    end
                                end
                            else
                                # If there isn't one argument, there are two, and it's
                                # the x*y case. Check whether log(x) or log(y) already exist
                                ID_1 = findfirst(x -> isequal(x.rhs, log(args[1])), factored)
                                ID_2 = findfirst(x -> isequal(x.rhs, log(args[2])), factored)

                                if isnothing(ID_1) && isnothing(ID_2)
                                    # Neither log(x) nor log(y) exist, so we need to make them
                                    newsym = gensym(:aux)
                                    newsym = Symbol(string(newsym)[3:5] * string(newsym)[7:end])
                                    newvar1 = genvar(newsym)
                                    insert!(factored, index0, Equation(Symbolics.value(newvar1), log(args[1])))
                                    newsym = gensym(:aux)
                                    newsym = Symbol(string(newsym)[3:5] * string(newsym)[7:end])
                                    newvar2 = genvar(newsym)
                                    insert!(factored, index0, Equation(Symbolics.value(newvar2), log(args[2])))
                                    @eval $factored[$index0+2] = $factored[$index0+2].lhs ~ $factored[$index0].lhs + $factored[$index0+1].lhs
                                elseif isnothing(ID_1)
                                    # Only log(y) exists, so we need to make log(x)
                                    newsym = gensym(:aux)
                                    newsym = Symbol(string(newsym)[3:5] * string(newsym)[7:end])
                                    newvar1 = genvar(newsym)
                                    insert!(factored, index0, Equation(Symbolics.value(newvar1), log(args[1])))
                                    # Find ID_2 again since it may have shifted
                                    ID_2 = findfirst(x -> isequal(x.rhs, log(args[2])), factored)
                                    @eval $factored[$index0+1] = $factored[$index0+1].lhs ~ $factored[$index0].lhs + $factored[$ID_2].lhs
                                elseif isnothing(ID_2)
                                    # Only log(x) exists, so we need to make log(y)
                                    newsym = gensym(:aux)
                                    newsym = Symbol(string(newsym)[3:5] * string(newsym)[7:end])
                                    newvar2 = genvar(newsym)
                                    insert!(factored, index0, Equation(Symbolics.value(newvar2), log(args[2])))
                                    # Find ID_1 again since it may have shifted
                                    ID_1 = findfirst(x -> isequal(x.rhs, log(args[1])), factored)
                                    @eval $factored[$index0+1] = $factored[$index0+1].lhs ~ $factored[$ID_1].lhs + $factored[$index0].lhs
                                else
                                    # Both log(x) and log(y) exist, so we can use those auxiliary variables
                                    @eval $factored[$index0] = $factored[$index0].lhs ~ $factored[$ID_1].lhs + $factored[$ID_2].lhs
                                end
                            end
                            break
                        end
                    end
                end
    
                # 5) (x1^a)^b = x1^(a*b) [EAGO paper]
                if exprtype(factored[index0].rhs)==POW
                    # Two cases we care about: either we've found (thing1)^(thing2) 
                    # and thing1 is POW, or (thing1)^(thing2) and thing2 is POW.
                    index1 = findfirst(x -> isequal(x.lhs, factored[index0].rhs.base), factored)
                    if !isnothing(index1) && exprtype(factored[index1].rhs)==POW
                        # If we have:
                        # aux1 = x1^[any1] (index1)
                        # aux2 = aux1^[any2] (index0)
                        #
                        # Go to:
                        # aux1 = x1^[any1]
                        # aux3 = [any1]*[any2]
                        # aux2 = x1^aux3
                        # 
                        # Or, if [any1] and [any2] are both numbers,
                        # aux1 = x1^[any1]
                        # aux2 = x1^([any1]*[any2])
                        scan_flag = true
    
                        # Identify the lower-level base, "x1"
                        factor_base = factored[index1].rhs.base
    
                        # Identify the two exponents
                        exp1 = factored[index0].rhs.exp
                        exp2 = factored[index1].rhs.exp
    
                        # Create new terms based on the exponents
                        if typeof(exp1) <: Real && typeof(exp2) <: Real
                            ID_1 = findfirst(x -> isequal(x.rhs, factor_base^(exp1*exp2)), factored)
                            if isnothing(ID_1)
                                @eval $factored[$index0] = $factored[$index0].lhs ~ $factor_base^($exp1*$exp2)
                            else
                                for i in eachindex(factored)
                                    @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $factored[$ID_1].lhs))
                                end
                                deleteat!(factored, index0)
                            end
                        else
                            ID_1 = findfirst(x -> isequal(x.rhs, exp1*exp2), factored)
                            if isnothing(ID_1)
                                newsym = gensym(:aux)
                                newsym = Symbol(string(newsym)[3:5] * string(newsym)[7:end])
                                newvar = genvar(newsym)
                                insert!(factored, index0, Equation(Symbolics.value(newvar), exp1*exp2))
                                @eval $factored[$index0+1] = $factored[$index0+1].lhs ~ $factor_base^$factored[$index0].lhs
                            else
                                ID_2 = findfirst(x -> isequal(x.rhs, factor_base^factored[ID_1].lhs), factored)
                                if isnothing(ID_2)
                                    @eval $factored[$index0] = $factored[$index0].lhs ~ $factor_base^$factored[$ID_1].lhs
                                else
                                    for i in eachindex(factored)
                                        @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $factored[$ID_2].lhs))
                                    end
                                    deleteat!(factored, index0)
                                end
                            end
                        end
                        break
                    end
                    # Try again with the second case
                    index1 = findfirst(x -> isequal(x.lhs, factored[index0].rhs.exp), factored)
                    if !isnothing(index1) && exprtype(factored[index1].rhs)==POW
                        # If we have:
                        # aux1 = [any1]^[any2] (index1)
                        # aux2 = x1^[aux1] (index0)
                        #
                        # Go to:
                        # aux1 = [any1]^[any2]
                        # aux3 = [any1]*[any2]
                        # aux2 = x1^aux3
                        scan_flag = true
    
                        # Identify the lower-level base, "x1"
                        factor_base = factored[index0].rhs.base
    
                        # Identify the correct exponents
                        exp1 = factored[index1].rhs.base
                        exp2 = factored[index1].rhs.exp
    
                        # Create new terms based on the exponents
                        ID_1 = findfirst(x -> isequal(x.rhs, exp1*exp2), factored)
                        if isnothing(ID_1)
                            newsym = gensym(:aux)
                            newsym = Symbol(string(newsym)[3:5] * string(newsym)[7:end])
                            newvar = genvar(newsym)
                            insert!(factored, index0, Equation(Symbolics.value(newvar), exp1*exp2))
                            @eval $factored[$index0+1] = $factored[$index0+1].lhs ~ $factor_base^$factored[$index0].lhs
                        else
                            ID_2 = findfirst(x -> isequal(x.rhs, factor_base^factored[ID_1].lhs), factored)
                            if isnothing(ID_2)
                                @eval $factored[$index0] = $factored[$index0].lhs ~ $factor_base^$factored[$ID_1].lhs
                            else
                                for i in eachindex(factored)
                                    @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $factored[$ID_2].lhs))
                                end
                                deleteat!(factored, index0)
                            end
                        end
                        break
                    end
                end
    
                # 7) log(inv(x1)) = -log(x1) [EAGO paper]
                if exprtype(factored[index0].rhs)==TERM
                    if factored[index0].rhs.f==log
                        index1 = findfirst(x -> isequal(x.lhs, arguments(factored[index0].rhs)[]), factored)
                        if !(isnothing(index1)) && exprtype(factored[index1].rhs)==DIV
                            # Since division is split into multiplication and inversion,
                            # if we have division, that means we have 1/[thing]. 
                            # We have:
                            # aux1 = 1 / [thing]
                            # aux2 = log(aux1)
                            #
                            # We want to convert this to:
                            # aux1 = 1 / [thing]
                            # aux3 = log([thing])
                            # aux2 = -aux3
                            scan_flag = true
                            ID_1 = findfirst(x -> isequal(x.rhs, log(factored[index1].rhs.den)), factored)
                            if isnothing(ID_1)
                                newsym = gensym(:aux)
                                newsym = Symbol(string(newsym)[3:5] * string(newsym)[7:end])
                                newvar = genvar(newsym)
                                insert!(factored, index0, Equation(Symbolics.value(newvar), log(factored[index1].rhs.den)))
                                @eval $factored[$index0+1] = $factored[$index0+1].lhs ~ -$factored[$index0].lhs
                            else
                                ID_2 = findfirst(x -> isequal(x.rhs, -factored[ID_1].lhs), factored)
                                if isnothing(ID_2)
                                    @eval $factored[$index0] = $factored[$index0].lhs ~ -$factored[$ID_1].lhs
                                else
                                    for i in eachindex(factored)
                                        @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $factored[$ID_2].lhs))
                                    end
                                    deleteat!(factored, index0)
                                end
                            end
                            break
                        end
                    end
                end
                
                # 8) CONST1*CONST2*x1 = (CONST1*CONST2)*x1
                if exprtype(factored[index0].rhs)==MUL
                    # We need to check all the multiplication arguments. We create
                    # a list of all the args using a recursive function.
                    args = pull_mult(factored, index0)
    
                    # Count the number of args that are Real-valued
                    real_count = 0
                    if length(args) > 1
                        for arg in args
                            if typeof(arg) <: Real
                                real_count += 1
                            end
                        end
                    end
    
                    # If we have more than 1 Real-valued argument, we'll apply this rule. 
                    # Otherwise, we don't do anything.
                    if real_count > 1
                        # We're going to create one big term and factor it, since
                        # the arrangement of multiplications in the original factorization
                        # might not line up with the changes
                        scan_flag = true
                        # Create the initial object
                        new_expr = Num(1.0)
                        # Multiply by all the terms
                        for arg in args
                            new_expr *= arg
                        end
                        # Create a factorization of this new expr
                        new_factorization = factor(new_expr, split_div=true)

                        
                        # Scan through the new factorization to see if we can merge elements
                        # with the original factored list
                        done = false
                        last_match = nothing
                        while !done
                            done = true
                            for i in eachindex(new_factorization)
                                ID_1 = findfirst(x -> isequal(x.rhs, new_factorization[i].rhs), factored)
                                if !isnothing(ID_1)
                                    # Match was found in the main factorization. Remove this element
                                    # from the new factorization and replace all references to it with
                                    # the auxiliary variable of ID_1
                                    for j in eachindex(new_factorization)
                                        new_factorization[j] = new_factorization[j].lhs ~ substitute(new_factorization[j].rhs, Dict(new_factorization[i].lhs => factored[ID_1].lhs))
                                    end
                                    last_match = factored[ID_1].lhs
                                    deleteat!(new_factorization, i)
                                    done = false
                                    break
                                end
                            end
                        end
                        # Add remaining elements of the factorization to the location of index0
                        if length(new_factorization) >= 1
                            for j = length(new_factorization)-1:-1:1
                                insert!(factored, index0, new_factorization[j])
                            end
                            @eval $factored[$index0+length($new_factorization)-1] = $factored[$index0+length($new_factorization)-1].lhs ~ $new_factorization[end].rhs
                        else
                            # Or, if no elements exist, replace all instances of this new auxiliary variable with the last known match
                            for i in eachindex(factored)
                                @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $last_match))
                            end
                            deleteat!(factored, index0)
                        end
                        break
                    end
    
                end
    
                # 9) 1 / (1 + exp(-x)) = Sigmoid(x)
                if exprtype(factored[index0].rhs)==DIV
                    # Get index of the denominator
                    index1 = findfirst(x -> isequal(x.lhs, factored[index0].rhs.den), factored)
                    if !isnothing(index1) && exprtype(factored[index1].rhs)==ADD
                        # Make sure denominator is (1 + [...])
                        if factored[index1].rhs.coeff==1
                            for key in keys(factored[index1].rhs.dict)
                                # Get the index of the other term
                                index2 = findfirst(x -> isequal(x.lhs, key), factored)
                                if !isnothing(index2) && exprtype(factored[index2].rhs)==TERM
                                    # If the other term is a "TERM" type, make sure it's "exp"
                                    if factored[index2].rhs.f==exp
                                        scan_flag = true
                                        # Get the index of the argument of "exp()"
                                        index3 = findfirst(x -> isequal(x.lhs, arguments(factored[index2].rhs)[]), factored)
                                        if isnothing(index3)
                                            # exp's argument is a base-level variable. Add a negation of that variable
                                            # and take the sigmoid of that.
                                            # Check first to make sure these terms aren't already being used
                                            ID_1 = findfirst(x -> isequal(x.rhs, -1*arguments(factored[index2].rhs)[]), factored)
                                            if isnothing(ID_1)
                                                newsym = gensym(:aux)
                                                newsym = Symbol(string(newsym)[3:5] * string(newsym)[7:end])
                                                newvar = genvar(newsym)
                                                insert!(factored, index0, Equation(Symbolics.value(newvar), -1*arguments(factored[index2].rhs)[]))
                                                @eval $factored[$index0+1] = $factored[$index0+1].lhs ~ SCMC_sigmoid($newvar)
                                            else
                                                ID_2 = findfirst(x -> isequal(x.rhs, SCMC_sigmoid(factored[ID_1].lhs)), factored)
                                                if isnothing(ID_2)
                                                    @eval $factored[$index0] = $factored[$index0].lhs ~ SCMC_sigmoid($factored[$ID_1].lhs)
                                                else
                                                    for i in eachindex(factored)
                                                        @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $factored[$ID_2].lhs))
                                                    end
                                                    deleteat!(factored, index0)
                                                end
                                            end
                                        elseif (exprtype(factored[index3].rhs)==MUL) && (factored[index3].rhs.coeff==-1)
                                            # Ideal case. The argument of exp() is already negated. Simply make
                                            # the RHS "SCMC_sigmoid(arg)" (if it doesn't already exist)
                                            arg = arguments(factored[index3].rhs)[2]
                                            arg = isequal(arg,-1) ? arguments(factored[index3].rhs)[1] : arg
                                            ID_1 = findfirst(x -> isequal(x.rhs, SCMC_sigmoid(arg)), factored)
                                            if isnothing(ID_1)
                                                @eval $factored[$index0] = $factored[$index0].lhs ~ SCMC_sigmoid($arg)
                                            else
                                                for i in eachindex(factored)
                                                    @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $factored[$ID_1].lhs))
                                                end
                                                deleteat!(factored, index0)
                                            end
                                        else
                                            # exp's argument is not yet negated. Add a negation of that variable and 
                                            # take the sigmoid of that, if needed.
                                            ID_1 = findfirst(x -> isequal(x.rhs, -1*factored[index3].lhs), factored)
                                            if isnothing(ID_1)
                                                newsym = gensym(:aux)
                                                newsym = Symbol(string(newsym)[3:5] * string(newsym)[7:end])
                                                newvar = genvar(newsym)
                                                insert!(factored, index0, Equation(Symbolics.value(newvar), -1*factored[index3].lhs))
                                                @eval $factored[$index0+1] = $factored[$index0+1].lhs ~ SCMC_sigmoid($newvar)
                                            else
                                                ID_2 = findfirst(x -> isequal(x.rhs, SCMC_sigmoid(factored[ID_1].lhs)), factored)
                                                if isnothing(ID_2)
                                                    @eval $factored[$index0] = $factored[$index0].lhs ~ SCMC_sigmoid($factored[$ID_1].lhs)
                                                else
                                                    for i in eachindex(factored)
                                                        @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $factored[$ID_2].lhs))
                                                    end
                                                    deleteat!(factored, index0)
                                                end
                                            end
                                        end
                                        break
                                    end
                                end
                            end
                        end
                    end
                end
                
                # 10) sin(x) = cos(x - pi/2)
                if exprtype(factored[index0].rhs)==TERM
                    if factored[index0].rhs.f==sin
                        # We found sin(arg). Check if (arg - pi/2) exists,
                        # and if so, also check if cos(arg - pi/2) exists.
                        scan_flag = true
                        index1 = findfirst(x -> isequal(x.rhs, arguments(factored[index0].rhs)[] - pi/2), factored)
                        if !isnothing(index1)
                            index2 = findfirst(x -> isequal(x.rhs, cos(factored[index1].lhs)), factored)
                            if !isnothing(index2)
                                # cos(arg - pi/2) exists already (index2). Remove all reference to index0 and replace with index2
                                for i in eachindex(factored)
                                    @eval $factored[$i] = $factored[$i].lhs ~ substitute($factored[$i].rhs, Dict($factored[$index0].lhs => $factored[$index2].lhs))
                                end
                                deleteat!(factored, index0)
                            else
                                # arg - pi/2 exists already (index1), but not cos(arg - pi/2). Change
                                # index0 to be cos of index1.lhs instead of sin of arg
                                @eval $factored[$index0] = $factored[$index0].lhs ~ cos($factored[$index1].lhs)
                            end
                        else
                            # (arg - pi/2) doesn't exist, so we need to create it
                            newsym = gensym(:aux)
                            newsym = Symbol(string(newsym)[3:5] * string(newsym)[7:end])
                            newvar = genvar(newsym)
                            insert!(factored, index0, Equation(Symbolics.value(newvar), arguments(factored[index0].rhs)[] - pi/2))
                            @eval $factored[$index0+1] = $factored[$index0+1].lhs ~ cos($newvar)
                        end
                        break
                    end
                end
            end
        end

        # Perform a cleaning step to eliminate any expressions that are no
        # longer used in the final result
        factored = cleaner(factored)
    end

    return factored
end

# This function "cleans" a factored expression by removing any terms
# that do not participate in the final expression.
function cleaner(list::Vector{Equation})
    IDs = zeros(Bool, length(list))
    # Start with the final term, and begin looking for participating
    # terms.
    cleaner(list[end].rhs, list, IDs)
    return list[IDs]
end
function cleaner(ex::BasicSymbolic, list::Vector{Equation}, IDs::Vector{Bool})
    # Add this expression to the list of IDs
    ID = findfirst(x -> isequal(x.rhs, ex), list)
    isnothing(ID) && return nothing
    IDs[ID] = true

    # Find all variables in "ex", add them to the list, and check those
    # expressions as well.
    if exprtype(ex)==SYM
        # No need to add anything else
        nothing
    elseif exprtype(ex)==TERM
        new_ID = findfirst(x -> isequal(x.lhs, ex.arguments[]), list)
        !isnothing(new_ID) && cleaner(list[new_ID].rhs, list, IDs)
    elseif exprtype(ex)==ADD
        for (key, val) in ex.dict
            new_ID = findfirst(x -> isequal(x.lhs, key), list)
            !isnothing(new_ID) && cleaner(list[new_ID].rhs, list, IDs)
        end
    elseif exprtype(ex)==MUL
        for (key, val) in ex.dict
            new_ID = findfirst(x -> isequal(x.lhs, key), list)
            !isnothing(new_ID) && cleaner(list[new_ID].rhs, list, IDs)
        end
    elseif exprtype(ex)==DIV
        new_ID = findfirst(x -> isequal(x.lhs, ex.den), list)
        isnothing(new_ID) && return nothing
        cleaner(list[new_ID].rhs, list, IDs) # Assuming no numerator since we're splitting division
    elseif exprtype(ex)==POW
        new_ID = findfirst(x -> isequal(x.lhs, ex.base), list)
        !isnothing(new_ID) && cleaner(list[new_ID].rhs, list, IDs)
        new_ID = findfirst(x -> isequal(x.lhs, ex.exp), list)
        !isnothing(new_ID) && cleaner(list[new_ID].rhs, list, IDs)
    else
        nothing
    end
    return nothing
end
function cleaner(ex::Real, list::Vector{Equation}, IDs::Vector{Bool})
    error("This probably shouldn't be reached. Submit an issue with the expression that brought you here.")
    # Add this expression to the list of IDs
    ID = findfirst(x -> isequal(x.rhs, ex), list)
    isnothing(ID) && return nothing
    IDs[ID] = true
    return nothing
end
cleaner(ex::Any, list::Vector{Equation}, IDs::Vector{Bool}) = nothing

# A helper function to pull all the multiplications connected to one
# factor.
function pull_mult(factors, index; args=[])
    if exprtype(factors[index].rhs)==MUL
        # Identify the new args
        for key in keys(factors[index].rhs.dict)
            new_index = findfirst(x -> isequal(x.lhs, key), factors)
            if !isnothing(new_index)
                args = pull_mult(factors, new_index, args=args)
            else
                push!(args, key)
            end
        end
        if !isone(factors[index].rhs.coeff)
            push!(args, factors[index].rhs.coeff)
        end
    else
        # Add the index of the factor to the list of args
        push!(args, factors[index].rhs)
    end
    return args
end


# A helper function to calculate the sparsity of a factorization, given
# the full gradlist
function detect_sparsity(factored, gradlist)
    # Idea is to check every element of "factored" and pull out a list of indices
    # in gradlist that the variables depend on.
    sparsity = Vector{Int}[]
    string_gradlist = string.(gradlist)

    for i in eachindex(factored)
        vars = string.(pull_vars(extract(factored, i)))
        indices = zeros(Int, length(vars))
        for j in eachindex(indices)
            indices[j] = findfirst(==(string(vars[j])), string_gradlist)
        end
        push!(sparsity, indices)
    end

    return sparsity
end

# This function writes the preamble of a CUDA kernel to the specified file.
function preamble_string(file_hash::String, inputs::Vector{String}, kernel_ID::Int, n_temp_vars::Int, dims::Int)
    kernel_inputs = String(inputs[1])
    for i = 2:length(inputs)
        kernel_inputs *= ", $(String(inputs[i]))"
    end

    buffer = Base.IOBuffer()

    write(buffer, "function f_$(file_hash)_$kernel_ID($(kernel_inputs))\n")
    write(buffer, "    idx = threadIdx().x + (blockIdx().x - Int32(1)) * blockDim().x\n")
    write(buffer, "    stride = blockDim().x * gridDim().x\n")
    write(buffer, "    col = Int32(1)\n")
    write(buffer, "    colmax = Int32($dims)\n")

    for i = 1:n_temp_vars
        write(buffer, "    temp$(i)_lo = 0.0\n")
        write(buffer, "    temp$(i)_hi = 0.0\n")
        write(buffer, "    temp$(i)_cv = 0.0\n")
        write(buffer, "    temp$(i)_cc = 0.0\n")
        write(buffer, "    temp$(i)_cvgrad = @MVector zeros(Float64, $dims)\n")
        write(buffer, "    temp$(i)_ccgrad = @MVector zeros(Float64, $dims)\n")
    end

    # A development flag for enabling @inbounds in the generated kernels. Although @inbounds
    # should speed up the kernels due to the number of indexed memory accesses, in my testing
    # it actually made the kernels SLOWER. This could be due to GPU-side optimizations that
    # expect bounds-checking, and perhaps the compiler can't do those optimizations if I get
    # rid of the bounds check... Might be worth investigating at some point, but leave as
    # "false" for now.
    enable_inbounds = false
    if enable_inbounds
        write(buffer, "    @inbounds while idx <= Int32(size($(inputs[1]),1))\n")
    else
        write(buffer, "    while idx <= Int32(size($(inputs[1]),1))\n")
    end

    return String(take!(buffer))
    
end

# This function deciphers the type of RHS and selects
# the correct operation to add to the kernel in `file`
function write_operation(file::IOStream, RHS::BasicSymbolic{Real}, inputs::Vector{Any}, gradlist::Vector{String}, sparsity::Vector{Int})
    if exprtype(RHS) == ADD
        if ~iszero(RHS.coeff)
            write(file, SCMC_cadd_kernel(inputs..., RHS.coeff, gradlist, sparsity))
        elseif inputs[1]==inputs[2]
            write(file, SCMC_add_to_kernel(inputs[1], inputs[3], gradlist, sparsity))
        elseif inputs[1]==inputs[3]
            write(file, SCMC_add_to_kernel(inputs[1], inputs[2], gradlist, sparsity))
        else
            write(file, SCMC_add_kernel(inputs..., gradlist, sparsity))
        end
    elseif exprtype(RHS) == MUL
        if isone(RHS.coeff)
            write(file, SCMC_mult_kernel(inputs..., gradlist, sparsity))
        elseif RHS.coeff==-1
            write(file, SCMC_negate_kernel(inputs..., gradlist, sparsity))
        else
            write(file, SCMC_cmul_kernel(inputs..., RHS.coeff, gradlist, sparsity))
        end
    elseif exprtype(RHS) == DIV
        if isone(RHS.num)
            write(file, SCMC_inv_kernel(inputs..., gradlist, sparsity))
        else
            close(file)
            error("Division shouldn't have a numerator at this point ($RHS)")
        end
    elseif exprtype(RHS) == POW
        if isinteger(RHS.exp)
            write(file, SCMC_int_power_kernel(inputs..., Int(RHS.exp), gradlist, sparsity))
        else
            write(file, SCMC_float_power_kernel(inputs..., RHS.exp, gradlist, sparsity))
        end
    elseif exprtype(RHS) == TERM
        if RHS.f==exp
            write(file, SCMC_exp_kernel(inputs..., gradlist, sparsity))
        elseif RHS.f==log
            write(file, SCMC_log_kernel(inputs..., gradlist, sparsity))
        elseif RHS.f==SCMC_sigmoid
            write(file, SCMC_sigmoid_kernel(inputs..., gradlist, sparsity))
        elseif RHS.f==sqrt
            write(file, SCMC_float_power_kernel(inputs..., 0.5, gradlist, sparsity))
        elseif RHS.f==cos
            write(file, SCMC_cos_kernel(inputs..., gradlist, sparsity))
        elseif RHS.f==abs
            write(file, SCMC_abs_kernel(inputs..., gradlist, sparsity))
        else
            close(file)
            error("Some function was used that we can't handle yet ($RHS)")
        end
    else
        close(file)
        error("Something's wrong? $RHS")
    end
    return nothing
end

# Once the kernel is almost finished, the last-used temporary variable needs
# to be written to the output matrix. This adds text to accomplish that task,
# and then concludes the kernel.
function postamble(outvar::String)
    buffer = Base.IOBuffer()
    write(buffer, "\n")
    write(buffer, "        OUT[idx,1] = $(outvar)_cv\n")
    write(buffer, "        OUT[idx,2] = $(outvar)_cc\n")
    write(buffer, "        OUT[idx,3] = $(outvar)_lo\n")
    write(buffer, "        OUT[idx,4] = $(outvar)_hi\n")
    write(buffer, "        col = Int32(1)\n")
    write(buffer, "        while col <= colmax\n")
    write(buffer, "            OUT[idx,end-2*colmax+col] = $(outvar)_cvgrad[col]\n")
    write(buffer, "            OUT[idx,end-1*colmax+col] = $(outvar)_ccgrad[col]\n")
    write(buffer, "            col += Int32(1)\n")
    write(buffer, "        end\n")
    write(buffer, "        idx += stride\n")
    write(buffer, "    end\n")
    write(buffer, "    return nothing\n")
    write(buffer, "end\n")
    return String(take!(buffer))
end

# A special case, for the quadratic kernel
function postamble_quadaff(vars::Vector{String}, gradlist::Vector{String})
    # Some quick calculation is needed to figure out convex/concave relaxations
    # and bounds. Subgradients are already set, however.

    # Basic introduction
    buffer = Base.IOBuffer()
    write(buffer, "        #####################################\n")
    write(buffer, "        ## Finalize Relaxations and Bounds ##\n")
    write(buffer, "        #####################################\n")
    write(buffer, "\n")
    write(buffer, "        temp1_cv += intercept_cv\n")
    write(buffer, "        temp1_cc += intercept_cc\n")
    write(buffer, "        temp1_lo += intercept_cv\n")
    write(buffer, "        temp1_hi += intercept_cc\n")
    write(buffer, "\n")

    # Adjust lower/upper bounds to account for every participating variable
    for var in vars
        ID = findfirst(==(var), gradlist)
        write(buffer, "        temp1_cv += temp1_cvgrad[$ID] * $(var)[idx,1]\n")
        write(buffer, "        temp1_cc += temp1_ccgrad[$ID] * $(var)[idx,1]\n")
        write(buffer, "        if temp1_cvgrad[$ID] >= 0.0\n")
        write(buffer, "            temp1_lo += temp1_cvgrad[$ID] * $(var)[idx,2]\n")
        write(buffer, "        else\n")
        write(buffer, "            temp1_lo += temp1_cvgrad[$ID] * $(var)[idx,3]\n")
        write(buffer, "        end\n")
        write(buffer, "        if temp1_ccgrad[$ID] >= 0.0\n")
        write(buffer, "            temp1_hi += temp1_ccgrad[$ID] * $(var)[idx,3]\n")
        write(buffer, "        else\n")
        write(buffer, "            temp1_hi += temp1_ccgrad[$ID] * $(var)[idx,2]\n")
        write(buffer, "        end\n")
        write(buffer, "\n")
    end

    # Write bounds, relaxations, and subgradients to the output
    write(buffer, "        OUT[idx,1] = temp1_cv\n")
    write(buffer, "        OUT[idx,2] = temp1_cc\n")
    write(buffer, "        OUT[idx,3] = temp1_lo\n")
    write(buffer, "        OUT[idx,4] = temp1_hi\n")
    write(buffer, "        col = Int32(1)\n")
    write(buffer, "        while col <= colmax\n")
    write(buffer, "            OUT[idx,end-2*colmax+col] = temp1_cvgrad[col]\n")
    write(buffer, "            OUT[idx,end-1*colmax+col] = temp1_ccgrad[col]\n")
    write(buffer, "            col += Int32(1)\n")
    write(buffer, "        end\n")
    write(buffer, "        idx += stride\n")
    write(buffer, "    end\n")
    write(buffer, "    return nothing\n")
    write(buffer, "end\n")
    return String(take!(buffer))
end

# This function creates the final CPU function that calls the 
# CUDA kernels previously generated.
function outro(expr_hash::String, kernel_nums::Vector{Int}, inputs::Vector{Vector{String}}, outputs::Vector{String}, blocks::Int32, indep_vars::Vector{Symbol})
    # Start by getting the inputs to this final CPU function
    final_inputs = "OUT"
    for i in eachindex(indep_vars)
        final_inputs *= ", $(string(indep_vars[i]))"
    end

    # First, count how many temporary variables we need. We can
    # re-use a temporary variable if it's no longer needed after
    # it's used. While doing this, modify the inputs/outputs
    # to reflect the new variable names
    temp_vars = 0
    if length(inputs) > 1
        temp_tracker = Int[]
        for i in eachindex(inputs)
            # We don't need a temporary variable for the final kernel;
            # it's just "OUT"
            if i==length(inputs)
                continue
            end

            # We need a temporary variable. Put one in the right place.
            if isempty(temp_tracker)
                # If the list is empty, add this index to the list
                push!(temp_tracker, i)
                temp_vars += 1

                # Change names as appropriate
                starting_name = outputs[i]
                final_name = "tempvar_$(length(temp_tracker))"
                for j in eachindex(outputs)
                    if isequal(outputs[j], starting_name)
                        outputs[j] = final_name
                    end
                end
                for j in eachindex(inputs)
                    for k in eachindex(inputs[j])
                        if isequal(inputs[j][k], starting_name)
                            inputs[j][k] = final_name
                        end
                    end
                end
            else
                # If the list has elements, check if any are unused
                # (i.e., they've been switched to 0)
                unused_ID = findfirst(x -> iszero(x), temp_tracker)
                if isnothing(unused_ID)
                    push!(temp_tracker, i)
                    temp_vars += 1

                    # Change names as appropriate
                    starting_name = outputs[i]
                    final_name = "tempvar_$(length(temp_tracker))"
                    for j in eachindex(outputs)
                        if isequal(outputs[j], starting_name)
                            outputs[j] = final_name
                        end
                    end
                    for j in eachindex(inputs)
                        for k in eachindex(inputs[j])
                            if isequal(inputs[j][k], starting_name)
                                inputs[j][k] = final_name
                            end
                        end
                    end
                else
                    temp_tracker[unused_ID] = i

                    # Change names as appropriate
                    starting_name = outputs[i]
                    final_name = "tempvar_$(unused_ID)"
                    for j in eachindex(outputs)
                        if isequal(outputs[j], starting_name)
                            outputs[j] = final_name
                        end
                    end
                    for j in eachindex(inputs)
                        for k in eachindex(inputs[j])
                            if isequal(inputs[j][k], starting_name)
                                inputs[j][k] = final_name
                            end
                        end
                    end
                end
            end

            # Check if we can get rid of any previous temporary variables
            for j in eachindex(temp_tracker)
                iszero(temp_tracker[j]) && continue
                used = false
                for k in i+1:length(inputs)
                    if outputs[temp_tracker[j]] in inputs[k]
                        used = true
                        break
                    end
                end
                if !used
                    temp_tracker[j] = 0
                end
            end
        end
    end

    # Start writing to the buffer
    buffer = Base.IOBuffer()
    write(buffer, "\n")
    write(buffer, "function f_$(expr_hash)($(final_inputs))\n")


    # A development flag for enabling @inbounds in the generated kernels. Although @inbounds
    # should speed up the kernels due to the number of indexed memory accesses, in my testing
    # it actually made the kernels SLOWER. This could be due to GPU-side optimizations that
    # expect bounds-checking, and perhaps the compiler can't do those optimizations if I get
    # rid of the bounds check... Might be worth investigating at some point, but leave as
    # "false" for now. This particular `enable_inbounds` section turns on checks for correct
    # matrix sizes in the main kernel. If these checks pass, it should be fine to use @inbounds
    # in the lower-level kernels without memory access errors.
    enable_inbounds = false
    if enable_inbounds
        # Verify that the sizes of the input matrices are all the same. Doing this check here
        # allows us to use @inbounds later.
        write(buffer, "    # Verify that the sizes of the input matrices are all the same.\n")
        for i in eachindex(indep_vars)
            if i==1
                write(buffer, "    @assert size($(string(indep_vars[i])))")
            else
                write(buffer, " == size($(string(indep_vars[i])))")
            end
        end
        write(buffer, "\n\n")

        # Verify that the output matrix is the right size
        write(buffer, "    # Verify that the output matrix is the right size.\n")
        write(buffer, "    @assert size(OUT) == (size($(string(indep_vars[1])),1), $(4+2*length(indep_vars)))\n\n")
    end

    # Check if we need temporary variables, and if so, add them now.
    if temp_vars > 0
        write(buffer, "    # Create necessary temporary variables.\n")
        for i = 1:temp_vars
            write(buffer, "    tempvar_$i = similar(OUT)\n")
        end
        write(buffer, "\n")
    end

    # Now, call the kernels
    temp_tracker = Int[]
    write(buffer, "    # Call the kernels in order.\n")
    for i in eachindex(inputs)
        # Write the inputs for this kernel (first argument is the mutated output)
        kernel_inputs = string(outputs[i])
        for j in eachindex(inputs[i])
            kernel_inputs *= ", $(string(inputs[i][j]))"
        end

        # Create the kernel call
        # Note: The number of threads should never be lower than 256.
        # Max registers per SM is 65536, max registers per thread is 255,
        # so with 256 threads per SM, each can use all 255 registers (256^2 = 65536)
        write(buffer, "    CUDA.@sync @cuda always_inline=true threads=256 blocks=$blocks f_$(expr_hash)_$(kernel_nums[i])($kernel_inputs)\n")
    end

    # If we wrote temporary variables, clear them from memory
    if temp_vars > 0
        write(buffer, "\n")
        write(buffer, "    # Clear temporary variables from memory.\n")
        for i = 1:temp_vars
            write(buffer, "    CUDA.unsafe_free!(tempvar_$i)\n")
        end
    end

    # And wrap up the string
    write(buffer, "end\n")
    return String(take!(buffer))
end

complexity(num::Num) = complexity(perform_substitutions(factor(num, split_div=true)))
function complexity(factorized::Vector{Equation})
    # Set up the complexity list
    complexity = zeros(Int, length(factorized))

    # If we need the factorized list to be sorted topologically, we can do:
    # edgelist, varids = eqn_edges(factorized)
    # g = SimpleDiGraph(edgelist)

    # We shouldn't need the list to be sorted topologically, since the final
    # term should always correspond to the full expression. So, we can start
    # at the final term and go backwards
    complexity[end] = _complexity(complexity, factorized, length(factorized))
    return complexity
end

function _complexity(complexity::Vector{Int}, factorized::Vector{Equation}, start::Int)
    # Given the term of index [start], return the complexity (# lines)
    # of that expression and its subexpressions.
    # Note: The range for each total_lines row is because the expressions
    #       are longer if we need to assume subgradients (for base-level
    #       variables). In a large majority of the expression, unless it's
    #       very simple, it should usually NOT be assuming subgradients,
    #       so it'll be the shorter version. Still, using this lower value
    #       will make complexity() slightly optimistic.
    total_lines = 0
    RHS = factorized[start].rhs
    if exprtype(RHS) == ADD
        total_lines += 34 # Ranges from 34--49
        for (key,val) in RHS.dict
            new_ID = findfirst(x -> isequal(x.lhs, key), factorized)
            if ~isnothing(new_ID)
                total_lines += _complexity(complexity, factorized, new_ID)
            end
        end
    elseif exprtype(RHS) == MUL
        if ~isone(RHS.coeff)
            total_lines += 36 # Ranges from 36--43
            for (key,val) in RHS.dict
                new_ID = findfirst(x -> isequal(x.lhs, key), factorized)
                if ~isnothing(new_ID)
                    total_lines += _complexity(complexity, factorized, new_ID)
                end
            end
        else
            total_lines += 340 # Ranges from 340--561
            for (key,val) in RHS.dict
                new_ID = findfirst(x -> isequal(x.lhs, key), factorized)
                if ~isnothing(new_ID)
                    total_lines += _complexity(complexity, factorized, new_ID)
                end
            end
        end
    elseif exprtype(RHS) == DIV
        total_lines += 688 # Ranges from 688-960
        new_ID = findfirst(x -> isequal(x.lhs, RHS.den), factorized)
        if !isnothing(new_ID)
            total_lines += _complexity(complexity, factorized, new_ID)
        end
    elseif exprtype(RHS) == POW
        if iseven(RHS.exp)
            total_lines += 200
        elseif RHS.exp==3
            total_lines += 300
        elseif isinteger(RHS.exp)
            total_lines += 450
        else
            total_lines += 190
        end
        new_ID = findfirst(x -> isequal(x.lhs, RHS.base), factorized)
        if !isnothing(new_ID)
            total_lines += _complexity(complexity, factorized, new_ID)
        end
    elseif exprtype(RHS) == TERM
        if RHS.f==exp
            total_lines += 212 # Ranges from 212--310
            new_ID = findfirst(x -> isequal(x.lhs, RHS.arguments[1]), factorized)
            if !isnothing(new_ID)
                total_lines += _complexity(complexity, factorized, new_ID)
            end
        elseif RHS.f==log
            total_lines += 378 # Ranges from 378--520
            new_ID = findfirst(x -> isequal(x.lhs, RHS.arguments[1]), factorized)
            if !isnothing(new_ID)
                total_lines += _complexity(complexity, factorized, new_ID)
            end
        elseif RHS.f==SCMC_sigmoid
            total_lines += 1579 # Ranges from 1579--2008
            new_ID = findfirst(x -> isequal(x.lhs, RHS.arguments[1]), factorized)
            if !isnothing(new_ID)
                total_lines += _complexity(complexity, factorized, new_ID)
            end
        elseif RHS.f==sqrt
            total_lines += 190
            new_ID = findfirst(x -> isequal(x.lhs, RHS.arguments[1]), factorized)
            if !isnothing(new_ID)
                total_lines += _complexity(complexity, factorized, new_ID)
            end
        elseif RHS.f==cos || RHS.f==sin
            total_lines += 300
            new_ID = findfirst(x -> isequal(x.lhs, RHS.arguments[1]), factorized)
            if !isnothing(new_ID)
                total_lines += _complexity(complexity, factorized, new_ID)
            end
        elseif RHS.f==abs
            total_lines += 280
            new_ID = findfirst(x -> isequal(x.lhs, RHS.arguments[1]), factorized)
            if !isnothing(new_ID)
                total_lines += _complexity(complexity, factorized, new_ID)
            end
        else
            error("Some function was used that we can't handle yet ($RHS)")
        end
    elseif exprtype(RHS) == SYM
        nothing
    else
        error("Unknown type")
    end
    complexity[start] = total_lines
    return total_lines
end

# Use extract() to get the fully substituted form of each index
# of `factorized`, and then call pull_vars() on those forms to
# get the counts of unique variables. 
function var_counts(factorized::Vector{Equation})
    vars = zeros(Int, length(factorized))
    for i in eachindex(factorized)
        vars[i] = length(pull_vars(extract(factorized, i)))
    end
    return vars
end


# A test function to identify structures within an expression. This may be useful]
# in the future to, e.g., re-use kernels for multiple operations (note that we would
# need to disable sparsity in re-used kernels, or have it so that all potential re-uses
# are covered by the sparsity list).
function structure(expr::BasicSymbolic{Real})
    # We want to uniquely identify a given expression, but in
    # a way that would allow any particular set of variables
    # to be used in place of the existing variables without
    # modifying the overall structure of the expression.

    vars = pull_vars(expr)
    reordered_vars = Num[]
    string_expr = string(expr)
    locs = [findfirst(string(vars[i]), string_expr) for i in eachindex(vars)]
    cleaned_locs = zeros(length(locs))
    cleaned_locs[(!).(isnothing.(locs))] .= getfield.(locs[(!).(isnothing.(locs))], 1)
    cleaned_locs[isnothing.(locs)] .= Inf
    complete = false
    count = 1
    while !complete
        complete = true
        min_info = findmin(cleaned_locs)
        if isfinite(min_info[1])
            complete = false
            push!(reordered_vars, vars[min_info[2]])
            regex = "$(string(vars[min_info[2]]))(?!\\d)"
            regex = replace(regex, r"\[" => "\\[")
            regex = replace(regex, r"\]" => "\\]")
            regex = replace(regex, r"\s" => "\\s")
            string_expr = replace(string_expr, Regex(regex) => "unique_var_$count")
            count += 1
            cleaned_locs[min_info[2]] = Inf
        end
    end
    return string_expr, reordered_vars
end

# A test function that can be used to genericize an expression.
# Base-level variables get called BASE_###, and auxiliary
# variables get called AUX_###. This might be useful for checking
# the structure, if kernels might be used for multiple subexpressions
generic(expr::Num) = generic(expr.val)
function generic(expr::BasicSymbolic{Real})
    output = string(expr)
    all_vars = pull_vars(expr)
    base_count = 1
    aux_count = 1
    for var in all_vars
        if startswith(string(var), "aux")
            if aux_count < 10
                output = replace(output, "$var" => "AUX_0$aux_count")
            else
                output = replace(output, "$var" => "AUX_$aux_count")
            end
            aux_count += 1
        else
            if base_count < 10
                output = replace(output, "$var" => "BASE_0$base_count")
            else
                output = replace(output, "$var" => "BASE_$base_count")
            end
            base_count += 1
        end
    end
    return output
end