text:
#=
Code generation functions for FV for Julia.
=#

# Extract the input args. This must match the arguments passed by the solver.
function handle_input_args_fv_julia(lorr, vors)
    # args = (var, eid, fid, grid_data, geo_factors, fv_info, t, dt)
    code =
"var =       args[1];
eid =       args[2];
fid =       args[3];
grid =      args[4];
geo_facs =  args[5];
fv_data =   args[6];
refel =     args[7];
time =      args[8];
dt =        args[9];
"
    return code;
end

# Allocate, compute, or fetch all needed values
function prepare_needed_values_fv_julia(entities, var, lorr, vors)
    # Only gather the pieces that are used. See the end of this function.
    if vors == "volume"
        # el, nodex, loc2glb, detJ, J
        piece_needed = zeros(Bool, 5);
    else # surface
        # els, nodex, loc2glb, cellx, frefelind, facex, face2glb, normal, face_detJ, area, vol_J
        piece_needed = zeros(Bool, 11);
    end
    need_deriv_matrix = false;
    need_deriv_dist = false;

    code = "";

    # First label the index values
    for i=1:length(indexers)
        if indexers[i].symbol === :elements || indexers[i].symbol === :cells
            # not needed
        else
            code *= "INDEX_VAL_"*string(indexers[i].symbol)*" = kwargs["*string(i)*"]\n";
        end
    end

    for i=1:length(entities)
        cname = make_coef_name(entities[i]);
        if is_unknown_var(entities[i], var) && lorr == LHS
            # TODO

        else
            # Is coefficient(number or function) or variable(array)?
            (ctype, cval) = get_coef_val(entities[i]);
            if ctype == -1
                # It was a special symbol like dt or FACENORMAL
                if entities[i].name == "FACENORMAL1"
                    code *= cname * " = normal["*string(entities[i].index)*"]; # normal vector component\n"
                    piece_needed[8] = true;
                else entities[i].name == "FACENORMAL2"
                    code *= cname * " = -normal["*string(entities[i].index)*"]; # reverse normal vector component\n"
                    piece_needed[8] = true;
                end
            elseif ctype == 0
                # It was a number, do nothing?
            elseif ctype == 1 # a constant wrapped in a coefficient
                # This generates something like: coef_k_i = 4;
                if length(entities[i].derivs) > 0
                    code *= cname * " = 0; # NOTE: derivative applied to constant coefficient = 0\n";
                else
                    code *= cname * " = " * string(cval) * ";\n";
                end

            elseif ctype == 2 # a coefficient function
                # This generates something like:
                ######################################
                # coef_n_i = zeros(refel.Np); # allocate
                # for coefi = 1:refel.Np
                #     coef_k_i[coefi] = (Femshop.genfunctions[cval]).func(x[1,coefi], x[2,coefi],x[3,coefi],time); # evaluate at nodes
                # end
                ######################################
                if vors == "surface" && length(entities[i].derivs) == 0
                    nodesymbol = "facex"
                    piece_needed[6] = true;
                else
                    nodesymbol = "nodex"
                    piece_needed[2] = true;
                end
                cargs = "("*nodesymbol*"[coefi], 0, 0, time)";
                if config.dimension == 2
                    cargs = "("*nodesymbol*"[1, coefi], "*nodesymbol*"[2, coefi], 0, time)";
                elseif config.dimension == 3
                    cargs = "("*nodesymbol*"[1, coefi], "*nodesymbol*"[2, coefi], "*nodesymbol*"[3, coefi], time)";
                end
                if vors == "volume"
                    code *= cname * " = zeros(refel.Np);\n";
                    code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Femshop.genfunctions["*string(cval)*"]).func" * cargs * " end\n";
                    # Apply any needed derivative operators. Interpolate at quadrature points.
                    if length(entities[i].derivs) > 0
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            code *= cname * " = RD"*string(entities[i].derivs[di])*" * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                        end
                        need_deriv_matrix = true;
                        piece_needed[5] = true;
                    end
                    # integrate over cell
                    code *= cname * " = (refel.wg .* detj)' * refel.Q * " * cname * "; # integrate over cell\n";
                    piece_needed[4] = true;

                else # surface
                    # Apply any needed derivative operators. Interpolate at quadrature points.
                    if length(entities[i].derivs) > 0
                        code *= cname * " = zeros(refel.Np);\n";
                        code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Femshop.genfunctions["*string(cval)*"]).func" * cargs * " end\n";
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            code *= cname * " = RD"*string(entities[i].derivs[di])*" * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                        end
                        code *= cname * " = " * cname * "[refel.face2local[frefelind[1]]]; # extract face values only.";
                        need_deriv_matrix = true;
                        piece_needed[[4,5,8]] .= true;

                    else # no derivatives, only need surface nodes
                        code *= cname * " = zeros(refel.Nfp[frefelind[1]]);\n";
                        code *= "for coefi = 1:refel.Nfp[frefelind[1]] " * cname * "[coefi] = (Femshop.genfunctions["*string(cval)*"]).func" * cargs * " end\n";
                        piece_needed[5] = true;
                    end
                    # integrate over face
                    if config.dimension == 1
                        # in 1d there is only one face node
                        code *= cname * " = " * cname * "[1]\n";
                    else
                        code *= cname * " = (refel.surf_wg[frefelind[1]] .* face_detJ)' * refel.surf_Q[frefelind[1]][:, refel.face2local[frefelind[1]]] * " * cname * " / area; # integrate over face\n";
                        piece_needed[9] = true;
                        piece_needed[10] = true;
                    end
                end

            elseif ctype == 3 # a known variable value
                piece_needed[1] = true;
                # This generates something like: coef_u_1 = copy((Femshop.variables[1]).values[1, gbl])
                cellside = 0; # 0 means no side flag
                for flagi=1:length(entities[i].flags)
                    if occursin("DGSIDE1", entities[i].flags[flagi]) || occursin("CELL1", entities[i].flags[flagi])
                        cellside = 1;
                    elseif occursin("DGSIDE2", entities[i].flags[flagi]) || occursin("CELL2", entities[i].flags[flagi])
                        cellside = 2;
                    end
                end
                if vors == "surface"
                    l2gsymbol = "els[1]"
                else
                    l2gsymbol = "el"
                end
                if cellside == 1
                    l2gsymbol = "els[1]"
                elseif cellside == 2
                    l2gsymbol = "els[2]"
                end

                if typeof(entities[i].index) <: Array
                    # It is an indexed variable
                    if length(entities[i].index) == 1
                        indstr = "INDEX_VAL_"*entities[i].index[1];
                    else
                        # There is more than one index. Need to form an expression for it.
                        indstr = "(INDEX_VAL_"*entities[i].index[1];
                        indices = variables[cval].indexer;
                        for indi=2:length(entities[i].index)
                            indstr *= " + ("*string(length(indices[indi-1].range))*"*(INDEX_VAL_"*entities[i].index[indi]*"-1)";
                        end
                        for indi=1:length(entities[i].index)
                            indstr *= ")";
                        end
                    end

                else
                    indstr = string(entities[i].index);
                end

                if vors == "volume"
                    # Apply any needed derivative operators.
                    if length(entities[i].derivs) > 0
                        # Need a derivative here...
                        # TODO
                    else
                        code *= cname * " = Femshop.variables["*string(cval)*"].values["*indstr*", "*l2gsymbol*"];\n";
                    end
                else
                    if length(entities[i].derivs) > 0
                        code *= cname * " = Femshop.variables["*string(cval)*"].values["*indstr*", els[2]] - Femshop.variables["*string(cval)*"].values["*indstr*", els[1]];\n";
                        code *= cname * " = (els[1] != els[2] && abs(normal["*string(entities[i].derivs[1])*"]) > 1e-10) ? "*cname*" ./ dxyz["*string(entities[i].derivs[1])*"]  : 0\n"
                        need_deriv_dist = true;
                        piece_needed[4] = true; # cellx
                        piece_needed[8] = true; # normal
                    else
                        if cellside == 0
                            # No side was specified, so use the average
                            code *= cname * " = 0.5 * (Femshop.variables["*string(cval)*"].values["*indstr*", els[1]] + Femshop.variables["*string(cval)*"].values["*indstr*", els[2]]);\n";
                        else
                            code *= cname * " = Femshop.variables["*string(cval)*"].values["*indstr*", "*l2gsymbol*"];\n";
                        end
                    end
                end

            elseif ctype == 4 # an indexed coefficient
                # This generates something like:
                ######################################
                # coef_n_1 = zeros(refel.Np); # allocate
                # for coefi = 1:refel.Np
                #     coef_k_1[coefi] = (Femshop.coefficients[cval]).value[INDEX_VAL_i].func(x[1,coefi], x[2,coefi],x[3,coefi],time); # evaluate at nodes
                # end
                ######################################
                if vors == "surface" && length(entities[i].derivs) == 0
                    nodesymbol = "facex"
                    piece_needed[6] = true;
                else
                    nodesymbol = "nodex"
                    piece_needed[2] = true;
                end
                cargs = "("*nodesymbol*"[coefi], 0, 0, time)";
                if config.dimension == 2
                    cargs = "("*nodesymbol*"[1, coefi], "*nodesymbol*"[2, coefi], 0, time)";
                elseif config.dimension == 3
                    cargs = "("*nodesymbol*"[1, coefi], "*nodesymbol*"[2, coefi], "*nodesymbol*"[3, coefi], time)";
                end

                indstr = "";
                for indi=1:length(entities[i].index)
                    if indi>1
                        indstr *= ",";
                    end
                    indstr *= "INDEX_VAL_"*entities[i].index[indi];
                end

                if vors == "volume"
                    code *= cname * " = zeros(refel.Np);\n";
                    code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Femshop.coefficients["*string(cval)*"]).value["*indstr*"].func" * cargs * " end\n";
                    # Apply any needed derivative operators. Interpolate at quadrature points.
                    if length(entities[i].derivs) > 0
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            code *= cname * " = RD"*string(entities[i].derivs[di])*" * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                        end
                        need_deriv_matrix = true;
                        piece_needed[5] = true;
                    end
                    # integrate over cell
                    code *= cname * " = (refel.wg .* detj)' * refel.Q * " * cname * "; # integrate over cell\n";
                    piece_needed[4] = true;

                else # surface
                    # Apply any needed derivative operators. Interpolate at quadrature points.
                    if length(entities[i].derivs) > 0
                        code *= cname * " = zeros(refel.Np);\n";
                        code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Femshop.coefficients["*string(cval)*"]).value["*indstr*"].func" * cargs * " end\n";
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            code *= cname * " = RD"*string(entities[i].derivs[di])*" * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                        end
                        code *= cname * " = " * cname * "[refel.face2local[frefelind[1]]]; # extract face values only.";
                        need_deriv_matrix = true;
                        piece_needed[[4,5,8]] .= true;

                    else # no derivatives, only need surface nodes
                        code *= cname * " = zeros(refel.Nfp[frefelind[1]]);\n";
                        code *= "for coefi = 1:refel.Nfp[frefelind[1]] " * cname * "[coefi] = (Femshop.coefficients["*string(cval)*"]).value["*indstr*"].func" * cargs * " end\n";
                        piece_needed[5] = true;
                    end
                    # integrate over face
                    if config.dimension == 1
                        # in 1d there is only one face node
                        code *= cname * " = " * cname * "[1]\n";
                    else
                        code *= cname * " = (refel.surf_wg[frefelind[1]] .* face_detJ)' * refel.surf_Q[frefelind[1]][:, refel.face2local[frefelind[1]]] * " * cname * " / area; # integrate over face\n";
                        piece_needed[9] = true;
                        piece_needed[10] = true;
                    end
                end

            end
        end # if coefficient
    end # entity loop

    if vors == "volume"
        # el, nodex, loc2glb, detJ, J
        needed_pieces = "";
        if piece_needed[1]
            needed_pieces *= "el = eid;\n"
        end
        if piece_needed[2] || piece_needed[3]
            needed_pieces *= "loc2glb = grid.loc2glb[:,eid]; # local to global map\n"
        end
        if piece_needed[2]
            needed_pieces *= "nodex = grid.allnodes[:,loc2glb]; # node coordinates\n"
        end
        if piece_needed[4]
            needed_pieces *= "detj = geo_factors.detJ[eid]; # geometric factors\n"
        end
        if piece_needed[5]
            needed_pieces *= "J = geo_factors.J[eid]; # geometric factors\n"
        end

        if need_deriv_matrix
            needed_pieces *=
"
# Note on derivative matrices:
# RQn are vandermond matrices for the derivatives of the basis functions
# with Jacobian factors. They are made like this.
# |RQ1|   | rx sx tx || Qx |
# |RQ2| = | ry sy ty || Qy |
# |RQ3|   | rz sz tz || Qz |

"
            if config.dimension == 1
                needed_pieces *= "(RQ1, RD1) = build_deriv_matrix(refel, J);\n";
                needed_pieces *= "TRQ1 = RQ1';\n"

            elseif config.dimension == 2
                needed_pieces *= "(RQ1, RQ2, RD1, RD2) = build_deriv_matrix(refel, J);\n";
                needed_pieces *= "(TRQ1, TRQ2) = (RQ1', RQ2');\n"

            elseif config.dimension == 3
                needed_pieces *= "(RQ1, RQ2, RQ3, RD1, RD2, RD3) = build_deriv_matrix(refel, J);\n";
                needed_pieces *= "(TRQ1, TRQ2, TRQ3) = (RQ1', RQ2', RQ3');\n"
            end
        end
        if need_deriv_dist
            needed_pieces *= "TODO: neighboring cell-based derivative for volume. generate_code_layer_fv_julia.jl - prepare_needed_values_fv_julia()"
        end

    else # surface
        needed_pieces =
"if grid.face2element[1, fid] == eid # The normal points out of e"*(piece_needed[8] ? "\n\tnormal = grid.facenormals[:, fid];" : "")*"
    neighbor = grid.face2element[2, fid];"*(piece_needed[5] ? "\n\tfrefelind = [grid.faceRefelInd[1,fid], grid.faceRefelInd[2,fid]]; # refel based index of face in both elements" : "")*"
else # The normal points into e"*(piece_needed[8] ? "\n\tnormal = -grid.facenormals[:, fid];" : "")*"
    neighbor = grid.face2element[1, fid];"*(piece_needed[5] ? "\n\tfrefelind = [grid.faceRefelInd[2,fid], grid.faceRefelInd[1,fid]]; # refel based index of face in both elements" : "")*"
end
if neighbor == 0 # This is a boundary face. For now, just compute as if neighbor is identical. BCs handled later.
    neighbor = eid;
end\n
";
        # els, nodex, loc2glb, cellx, frefelind, facex, face2glb, normal, face_detJ, area, vol_J
        if piece_needed[1]
            needed_pieces *= "els = (eid, neighbor); # indices of elements on both sides\n"
        end
        if piece_needed[2] || piece_needed[3]
            needed_pieces *= "loc2glb = (grid.loc2glb[:,eid], grid.loc2glb[:, neighbor]); # volume local to global\n"
        end
        if piece_needed[2]
            needed_pieces *= "nodex = (grid.allnodes[:,loc2glb[1][:]], grid.allnodes[:,loc2glb[2][:]]); # volume node coordinates\n"
        end
        if piece_needed[4]
            needed_pieces *= "cellx = (fv_data.cellCenters[:, eid], fv_data.cellCenters[:, neighbor]); # cell center coordinates\n"
        end
        if piece_needed[5] || piece_needed[6]
            needed_pieces *= "face2glb = grid.face2glb[:,:,fid];         # global index for face nodes for each side of each face\n"
        end
        if piece_needed[6]
            needed_pieces *= "facex = grid.allnodes[:, face2glb[:, 1]];  # face node coordinates\n"
        end
        if piece_needed[9]
            needed_pieces *= "face_detJ = geo_facs.face_detJ[fid];       # detJ on face\n"
        end
        if piece_needed[10]
            needed_pieces *= "area = geo_facs.area[fid];                 # area of face\n"
        end
        if piece_needed[11]
            needed_pieces *= "vol_J = (geo_facs.J[eid], geo_facs.J[neighbor]);\n"
        end

        if need_deriv_matrix
            needed_pieces *= "TODO: derivative matrices for surface. generate_code_layer_fv_julia.jl - prepare_needed_values_fv_julia()"
        end
        if need_deriv_dist
            needed_pieces *= "dxyz = norm(cellx[2] - cellx[1]) .* normal; # normal scaled by distance between cell centers\n"
        end
    end

    code = needed_pieces * "\n" * code;

    return code;
end

function make_elemental_computation_fv_julia(terms, var, dofsper, offset_ind, lorr, vors)
    # Here is where I make some assumption about the form of the expression.
    # Since it was expanded by the parser it should look like a series of terms: t1 + t2 + t3...
    # Where each term, if LHS, involves one unknown component and possibly some coefficients.
    code = "";

    # If there were no terms,
    if length(terms) < 1
        code = "return nothing # There were no terms to compute here";
    end

    # Allocate the vector or matrix to be returned if needed
    if dofsper > 1
        if lorr == RHS
            code *= "cell_average = zeros("*string(dofsper)*"); # Allocate for returned values.\n"
        else
            if vors == "volume"
                code *= "cell_matrix = zeros(refel.Np * "*string(dofsper)*", "*string(dofsper)*"); # Allocate for returned matrix.\n"
            else
                code *= "cell_matrix = zeros(refel.Nfp[frefelind[1]] * "*string(dofsper)*", "*string(dofsper)*"); # Allocate for returned matrix.\n"
            end
        end
    end

    # Separate the factors of each term and form the calculation
    if dofsper > 1
        # Subvectors for each component
        if lorr == LHS
            subvector = Array{String, 2}(undef, dofsper, dofsper);
        else # RHS
            subvector = Array{String, 1}(undef, dofsper);
        end
        for smi=1:length(subvector)
            subvector[smi] = "";
        end

        if typeof(var) <: Array
            for vi=1:length(var) # variables
                # Process the terms for this variable
                for ci=1:length(terms[vi]) # components
                    for i=1:length(terms[vi][ci])
                        (term_result, var_ind) = generate_term_calculation_fv_julia(terms[vi][ci][i], var, lorr);

                        # println(terms)
                        # println(terms[vi])
                        # println(terms[vi][ci])
                        # println(terms[vi][ci][i])
                        # println(term_result * " : "*string(test_ind)*", "*string(trial_ind))

                        # Find the appropriate subvector for this term
                        subveci = offset_ind[vi] + ci;
                        subvecj = var_ind;
                        if lorr == LHS
                            subvec_ind = subveci + dofsper * (subvecj-1);
                        else
                            subvec_ind = subveci;
                        end

                        if length(subvector[subvec_ind]) > 1
                            subvector[subvec_ind] *= " .+ " * term_result;
                        else
                            subvector[subvec_ind] = term_result;
                        end
                    end
                end

            end # vi

        else # only one variable
            # Process the terms for this variable
            for ci=1:length(terms) # components
                for i=1:length(terms[ci])
                    (term_result, var_ind) = generate_term_calculation_fv_julia(terms[ci][i], var, lorr);

                    # Find the appropriate submatrix for this term
                    subvec_ind = ci;

                    if length(subvector[subvec_ind]) > 1
                        subvector[subvec_ind] *= " .+ " * term_result;
                    else
                        subvector[subvec_ind] = term_result;
                    end
                end
            end

        end

        # Put the subvector together into cell_average or cell_matrix
        if lorr == LHS
            for emi=1:dofsper
                for emj=1:dofsper
                    if length(subvector[emi]) > 1
                        if vors == "volume"
                            rangei = "("*string(emi-1)*"*refel.Np + 1):("*string(emi)*"*refel.Np)";
                        else
                            rangei = "("*string(emi-1)*"*refel.Nfp[frefelind[1]] + 1):("*string(emi)*"*refel.Nfp[frefelind[1]])";
                        end

                        code *= "cell_matrix["*rangei*"] = " * subvector[emi, emj] * "\n";
                    end
                end
            end
            code *= "return cell_matrix;\n"

        else # RHS
            for emi=1:dofsper
                if length(subvector[emi]) > 1
                    code *= "cell_average["*string(emi)*"] = " * subvector[emi] * "\n";
                end
            end
            code *= "return cell_average;\n"
        end

    else # one dof
        terms = terms[1];
        if lorr == LHS
            if vors == "volume"
                result = "zeros(refel.Np)";
            else
                result = "zeros(refel.Nfp[frefelind[1]])";
            end

        else
            result = "[0]";
        end

        #process each term
        first_term = true;
        for i=1:length(terms)
            (term_result, var_ind) = generate_term_calculation_fv_julia(terms[i], var, lorr);

            if !(term_result == "0")
                if first_term
                    result = term_result;
                    first_term = false;
                else
                    result *= " .+ " * term_result;
                end
            end
        end
        code *= "cell_average = [" * result * "];\n";
        code *= "return cell_average;\n"
    end

    return code;
end

function generate_term_calculation_fv_julia(term, var, lorr)
    result = "";
    # Note: separate_factors return test and trial info, but FV will not have any test functions.
    # trial_part refers to the unknown variable part.
    # example:
    # 0.1*D1__u_1*_FACENORMAL1_1    ->  ???               on LHS
    #                             ->  0.1 * (coef_D1xu_1 .* normal[1])     on RHS
    if lorr == LHS
        (test_part, var_part, coef_part, test_ind, var_ind) = separate_factors(term, var);
        # # LHS: ??
        # TODO
        #
        # Seriously, need to do

    else
        (test_part, var_part, coef_part, test_ind, var_ind) = separate_factors(term);
        # RHS: coef_part
        if !(coef_part === nothing)
            result = string(replace_entities_with_symbols(coef_part));
        else
            result = "0";
        end
    end

    return (result, var_ind);
end

# Generate the assembly loop structures and insert the content
function generate_assembly_loop_fv_julia(indices)
    # Each of the indices must be passed to the functions in a named tuple.
    # Pass all defined indexers.
    index_args = "";
    for i=1:length(indexers)
        if i > 1
            index_args *= ", ";
        end
        index_args *= "index_val_"*string(indexers[i].symbol)*"=indexing_variable_"*string(indexers[i].symbol);
    end
    if length(indexers) == 0
        index_args = "empty_index = 1";
    end

    # Make sure the elemental loop is included
    elements_included = false;
    for i=1:length(indices)
        if indices[i] == "elements" || indices[i] == "cells"
            elements_included = true;
        end
    end
    # If elements were not included, make them the outermost loop
    if !elements_included
        indices = ["elements"; indices];
    end

    code =
"# Label things that were allocated externally
sourcevec = allocated_vecs[1];
fluxvec = allocated_vecs[2];
facefluxvec = allocated_vecs[3];
face_done = allocated_vecs[4];

"
    # generate the loop structures
    code *= "# Loops\n"
    loop_start = ""
    loop_end = ""
    ind_offset = ""
    prev_ind = ""
    for i=1:length(indices)
        if indices[i] == "elements" || indices[i] == "cells"
            loop_start *= "for eid in elemental_order\n";
            loop_start *= "face_done .= 0; # Reset face_done in elemental loop\n";
            loop_end *= "end # loop for elements\n";
        else
            loop_start *= "for indexing_variable_"*string(indices[i].symbol)*" in "*string(indices[i].range)*"\n";
            loop_end *= "end # loop for "*string(indices[i].symbol)*"\n";
            if length(prev_ind) == 0
                ind_offset *= "(indexing_variable_"*string(indices[i].symbol);
                prev_ind = string(length(indices[i].range));
            else
                ind_offset *= " + "*prev_ind*" * (indexing_variable_"*string(indices[i].symbol)*" - 1";
                prev_ind = string(length(indices[i].range));
            end
        end
    end
    for i=2:length(indices)
        ind_offset *= ")";
    end
    index_offset = "dofs_per_loop * ("*ind_offset*" - 1 + dofs_per_node * (eid - 1)) + 1"
    face_index_offset = "dofs_per_loop * ("*ind_offset*" - 1 + dofs_per_node * (fid - 1)) + 1"

    code *= loop_start;

    # insert the content
    code *= "
    # determine index in global vector
    index_offset = "*index_offset*";

    # Zero the result vectors for this element
    sourcevec[index_offset:(index_offset + dofs_per_loop-1)] .= 0;
    fluxvec[index_offset:(index_offset + dofs_per_loop-1)] .= 0;

    ##### Source integrated over the cell #####
    # Compute RHS volume integral
    if !(source_rhs === nothing)

        sourceargs = (var, eid, 0, grid_data, geo_factors, fv_info, refel, t, dt);
        source = source_rhs.func(sourceargs; "*index_args*") ./ geo_factors.volume[eid];

        # Add to global source vector
        sourcevec[index_offset:(index_offset + dofs_per_loop-1)] .= source;
    end

    ##### Flux integrated over the faces #####
    # Loop over this element's faces.
    for i=1:refel.Nfaces
        fid = grid_data.element2face[i, eid];

        # determine index in global face vector
        face_index_offset = "*face_index_offset*";

        if !(flux_rhs === nothing)
            if face_done[fid] < dofs_per_node
                face_done[fid] += dofs_per_loop;

                fluxargs = (var, eid, fid, grid_data, geo_factors, fv_info, refel, t, dt);
                flux = flux_rhs.func(fluxargs; "*index_args*") .* geo_factors.area[fid];

                # Add to global flux vector for faces
                facefluxvec[face_index_offset:(face_index_offset + dofs_per_loop-1)] .= flux;
                # Combine all flux for this element
                fluxvec[index_offset:(index_offset + dofs_per_loop-1)] .+= flux ./ geo_factors.volume[eid];

            else
                # This flux has either been computed or is being computed by another thread.
                # The state will need to be known before paralellizing, but for now assume it's complete.
                flux = -facefluxvec[face_index_offset:(face_index_offset + dofs_per_loop-1)];
                fluxvec[index_offset:(index_offset + dofs_per_loop-1)] .-= facefluxvec[face_index_offset:(face_index_offset + dofs_per_loop-1)] ./ geo_factors.volume[eid];
            end
        end

        # Boundary conditions are applied to flux
        fbid = grid_data.facebid[fid]; # BID of this face
        if fbid > 0
            facex = grid_data.allnodes[:, grid_data.face2glb[:,1,fid]];  # face node coordinates

            if typeof(var) <: Array
                dofind = 0;
                for vi=1:length(var)
                    compo = indexing_variable_speed #1:var[vi].total_components
                    dofind = dofind + 1;
                    if prob.bc_type[var[vi].index, fbid] == NO_BC
                        # do nothing
                    elseif prob.bc_type[var[vi].index, fbid] == FLUX
                        # compute the value and add it to the flux directly
                        # Qvec = (refel.surf_wg[grid_data.faceRefelInd[1,fid]] .* geo_factors.face_detJ[fid])' * (refel.surf_Q[grid_data.faceRefelInd[1,fid]])[:, refel.face2local[grid_data.faceRefelInd[1,fid]]]
                        # Qvec = Qvec ./ geo_factors.area[fid];
                        # bflux = FV_flux_bc_rhs_only(prob.bc_func[var[vi].index, fbid][compo], facex, Qvec, t, dofind, dofs_per_node) .* geo_factors.area[fid];

                        bflux = FVSolver.FV_flux_bc_rhs_only_simple(prob.bc_func[var[vi].index, fbid][compo], fv_info.faceCenters[:,fid], t);

                        fluxvec[index_offset + dofind-1] += (bflux - facefluxvec[face_index_offset + dofind-1]) ./ geo_factors.volume[eid];
                        facefluxvec[face_index_offset + dofind-1] = bflux;
                    else
                        printerr(\"Unsupported boundary condition type: \"*prob.bc_type[var[vi].index, fbid]);
                    end
                end

            else # one variable
                d=indexing_variable_speed;
                if prob.bc_type[var.index, fbid] == NO_BC
                    # do nothing
                elseif prob.bc_type[var.index, fbid] == FLUX
                    # compute the value and add it to the flux directly
                    # Qvec = (refel.surf_wg[grid_data.faceRefelInd[1,fid]] .* geo_factors.face_detJ[fid])' * (refel.surf_Q[grid_data.faceRefelInd[1,fid]])[:, refel.face2local[grid_data.faceRefelInd[1,fid]]]
                    # Qvec = Qvec ./ geo_factors.area[fid];
                    # bflux = FV_flux_bc_rhs_only(prob.bc_func[var.index, fbid][d], facex, Qvec, t, dofind, dofs_per_node) .* geo_factors.area[fid];

                    bflux = FVSolver.FV_flux_bc_rhs_only_simple(prob.bc_func[var.index, fbid][d], fv_info.faceCenters[:,fid], t);

                    fluxvec[index_offset] += (bflux - facefluxvec[face_index_offset]) ./ geo_factors.volume[eid];
                    facefluxvec[face_index_offset] = bflux;
                else
                    printerr(\"Unsupported boundary condition type: \"*prob.bc_type[var.index, fbid]);
                end
            end
        end# BCs

    end# face loop
"
    # close loops
    code *= loop_end;

    # finish
    code *=
"
if !(source_rhs === nothing)
    return sourcevec + fluxvec;
else
    return fluxvec;
end
"
    # println(code)
    return code;
end


text:
# import functions to be extended
import FoodCachingExperiments: add!, remove!, wait!, countfooditems
# import some types and functions
import FoodCachingExperiments: Tray, Food, FoodItem, InspectionObserver, MaintenanceDiet, observe!

# Define a DummyModel that can keep track of available caching trays, food items
# and inspection observers that record the number of inspections a bird makes to
# a certain tray.
struct DummyModel
    trays::Vector{Tray}
    foodtypes::Vector{FoodItem}
    observers::Vector{InspectionObserver}
end
DummyModel() = DummyModel([], [], [])

# For the DummyModel, adding different objects appends them to the appropriate list.
add!(m::DummyModel, f::FoodItem) = push!(m.foodtypes, f)
add!(m::DummyModel, t::Tray) = push!(m.trays, t)
add!(m::DummyModel, o::InspectionObserver) = push!(m.observers, o)
add!(::DummyModel, ::Type{MaintenanceDiet}) = nothing

# For the DummyModel, removing objects pops them from the respective lists.
remove!(m::DummyModel, f::FoodItem) = removefromlist!(m.foodtypes, f)
remove!(m::DummyModel, f::Food) = removefromlist!(m.foodtypes, x -> x.id == f)
remove!(m::DummyModel, t::Tray) = removefromlist!(m.trays, t)
remove!(m::DummyModel, o::InspectionObserver) = removefromlist!(m.observers, o)
remove!(::DummyModel, ::Type{MaintenanceDiet}) = nothing
function remove!(m::DummyModel, ::Type{Any})
    for list in (:trays, :foodtypes, :observers)
        empty!(getproperty(m, list))
    end
end
# helper functions for removal
removefromlist!(list, obj) = removefromlist!(list, x -> x == obj)
removefromlist!(list, f::Function) = removefromlist!(list, findfirst(f, list))
removefromlist!(list, i::Int) = splice!(list, i)
removefromlist!(::Any, ::Nothing) = nothing

# Count all food items that are freely available and cached in trays.
function countfooditems(m::DummyModel, typ)
	n = countfooditems(m.foodtypes, typ)
	for tray in m.trays
		n += countfooditems(tray, typ)
	end
	n
end

# The DummyModel removes random numbers of items from the freely available ones
# and caches random numbers of items in the available caching trays. The
# difference between removed and cached items gets interpreted as eaten items.
# Also, when trays are available the DummyModel inspects randomly one of them.
function wait!(m::DummyModel, delta, extracondition = nothing)
    for foodtype in m.foodtypes
        n = rand(1:20) # number of items to remove
        if n > foodtype.n  # trying to remove more items than available
            n = foodtype.n
            removefromlist!(m.foodtypes, foodtype)
        else
            foodtype.n -= n # decrement the number of avaiable items
        end
        for tray in m.trays
            n == 0 && break
            n_cache = rand(1:n)
            add!(tray, foodtype.id, n_cache) # cache some items
            n -= n_cache
        end
    end
    for o in m.observers
        if length(m.trays) > 0
            observe!(m.observers[1], rand(m.trays))
        end
    end
end

ext: jl


text:
"""
    Druid

A Druid SQL querying library.

See also: [Client](@ref), [execute](@ref)
"""
module Druid
using URIs
using HTTP: request, IOError, StatusError
using JSON
using Tables
using CSV

# Druid connection
export Client

# Queries
export execute, Timeseries, TopN, GroupBy, Scan, Search, TimeBoundary, SegmentMetadata, DatasourceMetadata, Sql

# DataSources
export Table, Lookup, Unioned, Inline, QuerySource, Join

# Granularities and Interval
export SimpleGranularity, DurationGranularity, PeriodGranularity, Interval

# Aggregatioons
export Count, SingleField, StringAgg, Grouping, Filtered

# PostAggregations
export Arithmetic, FieldAccess, FinalizingFieldAccess, ConstantPA, Greatest, Least, JavaScriptPA, HyperUniqueCardinality

# Dimension specs
export DefaultDS, ExtractionDS, ListFiltered, RegexFiltered, PrefixFiltered, LookupDS, Map, MapLookupDS

# ExtractionFunctions
export RegexEF, PartialEF, SearchQueryEF, SubstringEF, StrlenEF,
    TimeFormatEF, TimeParseEF, JavaScriptEF, RegisteredLookupEF, InlineLookupEF,
    CascadeEF, StringFormatEF, UpperEF, LowerEF, BucketEF

# Filters
export Selector, ColumnComparison, RegexF, AndF, OrF, NotF,
    JavaScriptF, Contains, InsensitiveContains, Fragment, SearchF, InF, Like,
    Bound, IntervalF, TrueF, ExtractionFilter

# Having Filters
export EqualTo, GreaterThan, LessThan, DimSelector, AndH, OrH, NotH

# TopN Metric Specs
export Numeric, Dimension, Inverted

# LimitSpec
export DefaultLS, OrderByColumn

# SQL Parameter
export Parameter

# Virtual columns
export Expression

# Export abstract types
export Query, Granularity, Aggregator, PostAggregator, Filter, HavingSpec, DataSource,
    DimensionSpec, ExtractionFunction, SearchQuerySpec, TopNMetricSpec, LimitSpec,
    VirtualColumn

"""
    Granularity

Supertype of all granularities.
"""
abstract type Granularity end

"""
    Aggregator

Supertype of all aggregators.
"""
abstract type Aggregator end

"""
    PostAggregator

Supertype of all post aggregators.
"""
abstract type PostAggregator end

"""
    ExtractionFunction

Supertype of all extraction functions.
"""
abstract type ExtractionFunction end

"""
    DimensionSpec

Supertype of all dimension specs.
"""
abstract type DimensionSpec end

"""
    SearchQuerySpec

Supertype of all search query specs.
"""
abstract type SearchQuerySpec end

"""
    Filter

Supertype of all filters.
"""
abstract type Filter end

"""
    HavingSpec

Supertype of all having filters.
"""
abstract type HavingSpec end

"""
    DataSource

Supertype of all data sources.
"""
abstract type DataSource end

"""
    TopNMetricSpec

Supertype of all topN metric specs.
"""
abstract type TopNMetricSpec end

"""
    LimitSpec

Supertype of all limit specs.
"""
abstract type LimitSpec end

"""
    VirtualColumn

Supertype of all virtual columns.
"""
abstract type VirtualColumn end

"""
    Query

Supertype of all queries.
"""
abstract type Query end

"""
    QueryResult

Supertype of all query results.
"""
abstract type QueryResult end

include("utils.jl")
include("client.jl")
include("intervals.jl")
include("granularities.jl")
include("dimensionspecs.jl")
include("filters.jl")
include("having.jl")
include("aggregations.jl")
include("post_aggregations.jl")
include("datasources.jl")
include("virtualcolumns.jl")
include("queries.jl")

end # module


text:
using FastGaussQuadrature

"""
    lagrange_1d(points, i, x)

Evaluate the Lagrange interpolation polynomial defined over nodal `points` with index `i`
at point `x`.
"""
function lagrange_1d(points, i, x)
    l = 1.0
    x_i = points[i]
    for x_j in points
        if x_j != x_i
            l *= (x - x_j) / (x_i - x_j)
        end
    end
    return l
end

"""
    lagrange_diff(points, i, x)

Evaluate the derivative of the Lagrange interpolation polynomial defined over nodal `points`
with index `i` at point `x`.
"""
function lagrange_diff(points, i, x)
    dl = 0.0
    x_i = points[i]
    for x_j in points
        if x_j != x_i
            p = 1.0
            for x_m in points
                if (x_m != x_i) & (x_m != x_j)
                    p *= (x - x_m) / (x_i - x_m)
                end
            end
            dl += 1.0/(x_i - x_j) * p
        end
    end
    return dl
end

"""
    get_quadpoints(n)

Compute quadrature points and weights for Gaussian quadrature
of order `n`.
The points (and thus the corresponding weights) are normalized to the range
``[0.0, 1.0]``.

Return a tuple of `points, weights`.
"""
function get_quadpoints(n)
    x, w = gausslegendre(n)
    # They are scaled in [-1.0, 1.0]
    # and we need [0.0, 1.0]
    (x .+ 1)./2, w ./ 2
end

"""
    Basis

A standard 1-dimensional basis of `order::Integer`
with `quadpoints::Array{Float64,1}` and
corresponding `quadweights::Array{Float64,1}`
Is the basis (pun intended) for tensor-product
bases.

    Basis(order::Integer, dimensions)

Initialize a basis of `order::Integer` and
dimensions `dimensions`.
"""
struct Basis
    quadpoints::Array{Float64,1}
    quadweights::Array{Float64,1}
    order::Int64
    dimensions::Int64

    function Basis(order, dimensions)
        quadpoints, quadweights = get_quadpoints(order)
        new(quadpoints, quadweights, order, dimensions)
    end
end

"""
    Base.length(basis::Basis)

Return number of points for `basis` in n-dimensions.
"""
#Base.length(basis::Basis) = length(basis.quadpoints)*basis.dimensions
#Base.length(basis::Basis) = length(basis.quadpoints) ^ basis.dimensions
Base.length(basis::Basis) = basis.order ^ basis.dimensions

"""
    Base.size(basis::Basis)

Return number of points for `basis` for each dimensions as tuple.
"""
function Base.size(basis::Basis)
    ntuple(basis.dimensions) do x
       length(basis.quadpoints)
    end
end

"""
    Base.size(basis::Basis, dim)

Return number of points for `basis` for dimensions `dim`.

function Base.size(basis::Basis, dim::Integer)
    length(basis.quadpoints)
end
"""

function Base.size(basis::Basis, dim::Integer)
    if dim > basis.dimensions
        println("Error: Given Dimension exceeds the Dimension of the problem = ",basis.dimensions)
    else
        length(basis.quadpoints)
    end
end

"""
    evaluate_basis(basis::Basis, coeffs, x)

Evaluate the `basis` with coefficients
`coeffs` at point `x`.

# So the input is coeffs , X point, I need to give basis
#multiplication of the basis is used , indicate tensor product
"""
function evaluate_basis(basis::Basis, coeffs, x)
    #coeffs[1] #code given for order 1

    #To collect result
    basisvalue=0

    # To make linear indexes
    # collect(): makes a array of tuples of indexes, same like meshgrid in matlab
    # vec(): converts matrix into vector
    linearindex = vec(collect(Iterators.product(1:basis.order, 1:basis.order)))

    for i in 1:length(basis)
        basisvalue = basisvalue + ( coeffs[i] * (  lagrange_1d(basis.quadpoints, linearindex[i][1], x[1]) *  lagrange_1d(basis.quadpoints, linearindex[i][2], x[2]) ) )
    end

    return basisvalue
end

"""
    project_to_reference_basis(fun, basis::Basis, ndofs)

Project the result of the function `fun` to coefficients
of the basis built of a tensor-product of `basis`.
The function `fun(x,y)`  takes in the ``x, y``-coordinates
and returns a vector with size `ndofs`.
The corresponding coefficients are returned.

#The reshape() is an inbuilt function in julia which is used
to return an array with the same data as the specified array,
but with different specified dimension sizes.

# It returns the matrix of shape[length(basis),ndof]

"""
function project_to_reference_basis(fun, basis::Basis, ndofs)
    #reshape(fun(0.5, 0.5), (1,ndofs)) #made for constant polynomial of order 1

    # make a matrix of shape [length(basis),ndof]
    coeffarray = zeros( ( length(basis) , ndofs ) )

    # To make linear indexes
    # collect(): makes a array of tuples of indexes, same like meshgrid in matlab
    # vec(): converts matrix into vector

    linearindex = vec(collect(Iterators.product(1:basis.order, 1:basis.order)))

    for i in 1:length(basis)
        coeffarray[i,:] = fun(  basis.quadpoints[linearindex[i][1]]  ,   basis.quadpoints[linearindex[i][2]]  ) # coefficients are nothing but values of function at that point
    end

    return coeffarray # No need to reshaping ,as matrix is made fully and filled it
end

"""
    massmatrix(basis, dimensions)

Return the mass-matrix for a `dimensions`-dimensional
tensor-product basis built up from the 1d-basis `basis`.
"""
function massmatrix(basis, dimensions)
    #ones(1,1)
    mass_matrix = Matrix(1.0I, basis.order^2, basis.order^2)
    ij = 1
    for j = 1:basis.order
        for i = 1:basis.order                                  #not generalising to dimensions>2
            mass_matrix[ij,ij] = basis.quadweights[i] * basis.quadweights[j]                        #eq 4.12 and 4.19
            ij += 1
        end
    end
    mass_matrix

end

"""
    derivativematrix(basis)

Returns the 2-dimensional derivative matrix for `basis`.
Multiplying this with flux-coefficients of shape
`(dimensions * basissize_2d, ndofs)` returns the
coefficients of the corresponding derivative.
"""
function derivativematrix(basis)
    #zeros(1, basis.dimensions * 1)

    # Finding Derivative matrix Try3
    # Finding Dx only, Derivative of the basis function in the direction of x
    Dx = zeros(basis.order^2, basis.order^2)
    Dy = zeros(basis.order^2, basis.order^2)
    linearindex = vec(collect(Iterators.product(1:basis.order, 1:basis.order)))

    for j in 1:basis.order^2
        for p in 1:basis.order^2
            Dx[j,p] = lagrange_diff(basis.quadpoints,linearindex[j][1],basis.quadpoints[linearindex[p][1]]) * lagrange_1d(basis.quadpoints, linearindex[j][2], basis.quadpoints[linearindex[p][2]])
        end
    end

    for j in 1:basis.order^2
        for p in 1:basis.order^2
            Dy[j,p] = lagrange_diff(basis.quadpoints,linearindex[j][2],basis.quadpoints[linearindex[p][2]]) * lagrange_1d(basis.quadpoints, linearindex[j][1], basis.quadpoints[linearindex[p][1]])
        end
    end
    hcat(Dx,Dy)
end

"""
    get_face_quadpoints(basis::Basis, face)

Return the quadrature points at the face `face` for basis `basis`.
"""
function get_face_quadpoints(basis::Basis, face)    # Is there a more efficient/ general way to do this?
    #ones(1.0, 1.0)

    face_quadpoints = zeros(basis.order,2)
    if face == left
        face_quadpoints[:,1] .= 0
        face_quadpoints[:,2] = basis.quadpoints
    elseif face == top
        face_quadpoints[:,1] = basis.quadpoints
        face_quadpoints[:,2] .= 1
    elseif face == right
        face_quadpoints[:,1] .= 1
        face_quadpoints[:,2] = basis.quadpoints
    elseif face == bottom
        face_quadpoints[:,1] = basis.quadpoints
        face_quadpoints[:,2] .= 0
    end
    face_quadpoints
end

function face_projection_matrix(basis, face)
    #ones(1,1)
    face_quadpoints = get_face_quadpoints(basis,face)

    # Matrix of required size
    P_ij = zeros( basis.order  , (basis.order)^2 )

    # Linear indexing
    linearindex = vec(collect(Iterators.product(1:basis.order, 1:basis.order)))

    for j = 1:(basis.order)^2
        for i = 1:basis.order
            P_ij[i,j] = lagrange_1d( basis.quadpoints, linearindex[j][1] , face_quadpoints[i,1] ) *  lagrange_1d( basis.quadpoints, linearindex[j][2], face_quadpoints[i,2] )
        end
    end
    return P_ij
end

"""
    evaluate_m_to_n_vandermonde_basis(basis)

Return the Vandermonde matrix that converts between the 2D-modal
(normalized) Legendre-basis and the 2D-nodal tensor-product basis built
with `basis`.
"""
function evaluate_m_to_n_vandermonde_basis(basis)
    Array{Float64}(I, length(basis), length(basis))
end


text:
using GitHubActions
using Documenter

makedocs(;
    modules=[GitHubActions],
    authors="Chris de Graaf <me@cdg.dev> and contributors",
    repo="https://github.com/julia-actions/GitHubActions.jl/blob/{commit}{path}#L{line}",
    sitename="GitHubActions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://julia-actions.github.io/GitHubActions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/julia-actions/GitHubActions.jl",
)


text:
const RE_PIPE_LAMINAR = 2100
const RE_PIPE_TURBULENT = 4000

"""
    FrictionFactor(Re::ReynoldsNumber,ϵD::RoughnessRatio)

Return the Darcy friction factor for a pipe flow with the given
Reynolds number and roughness ratio. The function chooses the laminar
formula (f = 64/Re) if Reynolds number is smaller than 2100 and
uses the Colebrook equation if Reynolds number is larger than 4000.
For transitional Reynolds numbers it uses the laminar formula, but
warns the user.
"""
function FrictionFactor(Re::ReynoldsNumber,ϵD::RoughnessRatio)
    _check_Re(Re)
    _check_roughness(ϵD)
    if _laminar(Re)
        return _f_laminar(Re)
    elseif _transitional(Re)
        @warn "Reynolds number is transitional. Using laminar relation."
        return _f_laminar(Re)
    elseif _turbulent(Re)
        return _f_turbulent(Re,ϵD)
    end
end

_f_laminar(Re::ReynoldsNumber) = FrictionFactor(64.0/Re)

function _f_turbulent(Re::ReynoldsNumber,ϵD::RoughnessRatio)
    err = 10.0
    fi = _f_haaland(Re,ϵD)
    f = fi
    while err > 1e-10
        f = _f_colebrook(Re,ϵD,fi)
        err = abs(f-fi)
        fi = f
    end
    f
end

function _f_colebrook(Re::ReynoldsNumber,ϵD::RoughnessRatio,f_in::FrictionFactor)
    _check_roughness(ϵD)
    t2 = Re == Inf ? 0.0 : 2.51/Re/sqrt(f_in)
    x = -2.0*log10(ϵD/3.7 + t2)
    return FrictionFactor(1.0/x^2)
end

function _f_haaland(Re::ReynoldsNumber,ϵD::RoughnessRatio)
    _check_roughness(ϵD)
    t2 = Re == Inf ? 0.0 : 6.9/Re
    x = -1.8*log10((ϵD/3.7)^1.11 + t2)
    return FrictionFactor(1.0/x^2)
end

_f_colebrook_fullyrough(ϵD::RoughnessRatio) = _f_colebrook(ReynoldsNumber(Inf),ϵD,FrictionFactor(Inf))

function _check_roughness(ϵD::RoughnessRatio)
    ϵD <= 0.05 || @warn "Roughness ratio is too large for Colebrook equation"
    ϵD >= 0.0 || throw(DomainError(ϵD,"Roughness ratio must be non-negative"))
    nothing
end

function _check_Re(Re::ReynoldsNumber)
    Re >= 0.0 || throw(DomainError(Re,"Reynolds number must be non-negative"))
    nothing
end

_laminar(Re::ReynoldsNumber) = Re <= RE_PIPE_LAMINAR

_transitional(Re::ReynoldsNumber) = (Re > RE_PIPE_LAMINAR && Re < RE_PIPE_TURBULENT)

_turbulent(Re::ReynoldsNumber) = Re >= RE_PIPE_TURBULENT

max_stars_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_issues_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_forks_repo_name: UCLAMAEThreads/ThermofluidQuantities


text:
module NearSurfaceTurbulenceModels

export
    SurfaceEnhancedModelConstant,
    SurfaceFluxDiffusivityBoundaryConditions,
    save_closure_parameters!

using Oceananigans,
      Oceananigans.BoundaryConditions,
      Oceananigans.TurbulenceClosures

using JLD2

struct SurfaceEnhancedModelConstant{T} <: Function
             C₀ :: T
             Δz :: T
             z₀ :: T
    enhancement :: T
    decay_scale :: T
end

"""
    SurfaceEnhancedModelConstant(Δz; FT=Float64, C₀=1/12, z₀=-Δz/2,
                                 enhancement=2, decay_scale=8Δz)

Returns a callable object representing a spatially-variable model constant
for an LES eddy diffusivity model with the surface-enhanced form

    ``C(z) = C₀ * (1 + enhancement * exp((z - z₀) / decay_scale)``

"""
function SurfaceEnhancedModelConstant(Δz; FT=Float64, C₀=1/12, z₀=-Δz/2,
                                      enhancement=2, decay_scale=8Δz)

    return SurfaceEnhancedModelConstant{FT}(C₀, Δz, z₀, enhancement, decay_scale)
end

@inline (C::SurfaceEnhancedModelConstant)(x, y, z) =
    C.C₀ * (1 + C.enhancement * exp((z - C.z₀) / C.decay_scale))

"""
    SurfaceEnhancedModelConstant(filename)

Reconstruct a `SurfaceEnhancedModelConstant` from file.
"""
function SurfaceEnhancedModelConstant(filename::String)
    file = jldopen(filename)

    Lz = file["grid/Lz"]
    Nz = file["grid/Nz"]
    Δz = Lz / Nz

    C₀ = file["closure/C₀"]
    z₀ = file["closure/z₀"]
    enhancement = file["closure/enhancement"]
    decay_scale = file["closure/decay_scale"]

    close(file)

    Cᴬᴹᴰ = SurfaceEnhancedModelConstant(Δz; C₀=1/12, z₀=-Δz/2, enhancement=2, decay_scale=8Δz)

    return Cᴬᴹᴰ
end

#####
##### Saving utilities for surface enhanced model constant
#####

save_closure_parameters!(args...) = nothing # fallback

const EnhancedAMD = AnisotropicMinimumDissipation{FT, PK, <:SurfaceEnhancedModelConstant} where {FT, PK}

function save_closure_parameters!(file, closure::EnhancedAMD)
    file["closure/C₀"] = closure.Cν.C₀
    file["closure/z₀"] = closure.Cν.z₀
    file["closure/enhancement"] = closure.Cν.enhancement
    file["closure/decay_scale"] = closure.Cν.decay_scale
    return nothing
end

#####
##### Utility for obtaining surface flux diffusivity boundary conditions
#####

function SurfaceFluxDiffusivityBoundaryConditions(grid, Qᵇ; Cʷ=0.1)
    w★ = (Qᵇ * grid.Lz)^(1/3) # surface turbulent velocity scaling
    κ₀ = Cʷ * grid.Δz * w★
    return DiffusivityBoundaryConditions(grid, top = BoundaryCondition(Value, κ₀))
end

end # module


text:
module Photometry

using Reexport

include("Aperture/Aperture.jl")

@reexport using .Aperture

end


text:
using HDF5, JLD
using PyCall
@pyimport matplotlib2tikz

# max_immediate_feasibility = 0.476
# mean_immediate_feasibility = 0.322

include("CenterOfMass.jl")
include("setup_parameters.jl")

save_plots = false
do_belief = false

plt[:close]("all")
pygui(!save_plots)

folder = "chance_constrained"

max_robots = 4

type TrialData
  belief::Histogram
  estimate
  error
  applied_f
  f_hat
  csqmi
  robots_measurement
  robots
  current_feasibility
  true_feasibility
  feasibility_no_limit
end

@load "$(folder)/data.fixed_probability.40" data_array csqmi_ratio thetas actuator_limit feasibility_constraint
n_trial = size(data_array, 1)

forces = map(x->x.applied_f, data_array)
num_measurement = map(x->length(x.robots_measurement), data_array)
num_robot = map(x->length(x.robots), data_array)
actuator_limits = actuator_limit * num_measurement
limits = forces - actuator_limits .>-1e-6

current_feasibility = map(x->x.current_feasibility, data_array)
true_feasibility = map(x->x.true_feasibility, data_array)

robots_feasible(robots, theta) =
  lifting_feasibility(attachment_ws[robots],
                      theta[3]*g*point_to_wrench(theta[1:2]),
                      actuator_limit)

current_lifting = Array{Bool}(size(data_array))
for ii = 1:size(data_array, 1)
  for jj = 1:size(data_array, 2)
    current_lifting[ii,jj] = robots_feasible(data_array[ii,jj].robots, thetas[:,ii])
  end
end

lifting_configurations = sum(current_lifting[:,end])

##########################################################################
# compare to behavior of best immediate combination to iterative selection
##########################################################################

best_initial = map(x->robots_feasible([1,3,5,8], thetas[:,x]), 1:size(thetas,2))
mean_best = mean(best_initial)

# print output
println("Found $(lifting_configurations) ($(mean(current_lifting[:,end]))) feasible lifting configurations in $(size(data_array,1)) trials")
println("$(sum(~true_feasibility[:,end])) trials had no feasible lifting configuration remaining")
println("The rate for the best initial configuration is $(mean_best)")

#######################
# plot normalized error
#######################
println("normalized error")

figure()
errors = map(x->x.error, data_array)
mean_errors = mean(errors, 1)

colors = ["b", "g", "r", "m"]

indices = collect(1:length(mean_errors))
plot(indices, mean_errors[:], color = "k", linewidth = 2.0)

cutoff = 0.15
n = 2
stds = std(errors, 1)[:] / sqrt(size(errors,1))
fill([indices;reverse(indices)], [mean_errors[:]+stds;reverse(mean_errors[:]-stds)],
  color = "k", alpha=0.2, linewidth=0.0)

if false
for ii = 1:size(errors, 1)
  plot(indices, errors[ii,:][:], color = "k", alpha = 0.2,
  linewidth = 1.0)

  #robot_added_inds = find(diff(num_robot[jj, ii, :][:]))+1
  #scatter(robot_added_inds, errors[jj, ii, robot_added_inds], color =
  #colors[jj], alpha = 0.2, marker = "1")

  for kk = 2:size(errors, 3)
    scatter(kk, errors[jj, ii, kk], color = colors[jj], alpha = 0.2)
  end
end
end

ylim(0.0, 0.5)
xlim(1.0, size(errors, 2))

xlabel("Iteration")
ylabel("Normalized error")

if save_plots
  matplotlib2tikz.save("$(folder)/convergence.tex", figureheight="\\figureheight",
    figurewidth="\\figurewidth")
end

#######################
# plot number of robots
#######################
println("number robots")

figure()

mean_robots = mean(num_robot, 1)

plot(indices, mean_robots[:], color = "k", linewidth = 2.0)

for ii = 1:size(num_robot, 1)
  plot(indices, num_robot[ii,:][:], color = "k", alpha = 0.2,
    linewidth = 1.0)
end

xlabel("Iteration")
ylabel("Num. robots")

if save_plots
  matplotlib2tikz.save("$(folder)/num_robots.tex", figureheight="\\figureheight",
    figurewidth="\\figurewidth")
end

##################
# plot feasibility
##################
println("feasibility")

figure()

mean_feasibility = mean(current_feasibility, 1)
#plot(indices, mean_feasibility[:], color = "k", linewidth = 2.0)
plot(indices, mean(true_feasibility,1)[:], color = "b")
standard_error(true_feasibility, "b")
plot(indices, mean(current_lifting,1)[:], color = "g")
standard_error(current_lifting,"g")

plot([1;size(errors, 2)], mean_best*ones(2), "--r")
standard_error(repmat(best_initial,1,size(errors,2)),"r")

#stds = std(current_feasibility, 1)[:] / sqrt(size(current_feasibility,1))
#fill([indices;reverse(indices)], [mean_feasibility[:]+stds;reverse(mean_feasibility[:]-stds)],
  #color = "k", alpha=0.2, linewidth=0.0)

if false
for ii = 1:size(current_feasibility, 1)
  plot(indices, current_feasibility[ii,:][:], color = "k", alpha = 0.2,
    linewidth = 1.0)
end
end

legend(["feasibility existence", "current feasibility", "max configuration"], loc="upper left")

xlabel("Iteration")
#ylabel("Feasibility probability")

if save_plots
  matplotlib2tikz.save("$(folder)/feasibility.tex", figureheight="\\figureheight",
    figurewidth="\\figurewidth")
end

#################################
# feasibility by number of robots
#################################

feasibility_by_number = zeros(n_trial,max_robots)
lifting_by_number = zeros(n_trial,max_robots)
for ii = 1:n_trial
  for jj = 1:max_robots
    ind = findfirst(num_robot[ii,:], jj)
    if ind > 0
      feasibility_by_number[ii,jj] = true_feasibility[ii,ind]
      lifting_by_number[ii,jj] = current_lifting[ii,ind]
    else
      # Treat the data pessimistically if not reaching a given number.
      # If not already in a lifting configuration, assume it becomes infeasible,
      # and if not already in a lifting configuration, assume such is not
      # reached
      feasibility_by_number[ii,jj] = lifting_by_number[ii,jj-1]
      lifting_by_number[ii,jj] = lifting_by_number[ii,jj-1]
    end
  end
end

figure()
plot(collect(1:max_robots), mean(feasibility_by_number,1)[:], color = "b")
standard_error(feasibility_by_number,"b")
plot(collect(1:max_robots), mean(lifting_by_number,1)[:], color = "g")
standard_error(lifting_by_number,"g")
plot([1;max_robots], mean_best*ones(2), "--r")
standard_error(repmat(best_initial,1,max_robots),"r")
legend(["feasibility existence", "current feasibility", "max configuration"], loc="upper left")

xlabel("Num. robots")
#ylabel("Feasibility probability")

if save_plots
  matplotlib2tikz.save("$(folder)/lifting.tex", figureheight="\\figureheight",
    figurewidth="\\figurewidth")
end

#############
# plot belief
#############
println("belief")

if save_plots && do_belief
  for ii = 1:size(errors,1)
    trial_folder = "$(folder)/$(ii)"
    mkpath(trial_folder)

    for jj = 1:size(errors,2)
      plt[:close]("all")
      println("belief: trial=$(ii), iter=$(jj)")
      data = data_array[ii,jj]
      robots = data.robots
      belief = data.belief
      robots_measurement = data.robots_measurement
      theta = thetas[:,ii]
      com_p = theta[1:2]
      mass = theta[3]

      figure()
      plot_attachment_points(attachment_ps)
      plot_occupied_points(attachment_ps[setdiff(robots, robots_measurement)])
      cloud = to_cloud(belief, interior_q)
      scaled_ps = (15*cloud[4,:]'/maximum(cloud[4,:]))
      scatter3D(cloud[1,:]', cloud[2,:]', cloud[3,:]', "z", scaled_ps.^2, "purple", alpha = 0.5)

      scatter3D([com_p[1]], [com_p[2]], [mass], "z", 200, "red", marker="*", alpha=1)
      scatter3D([com_p[1]], [com_p[2]], [0], "z", 200, "red", marker=(5,2,0), alpha=1)
      plot3D(com_p[1]*ones(2), com_p[2]*ones(2), [0.0; mass],
        color = "k", alpha = 1, linestyle = "dashdot")

      plot_measurement_points(attachment_ps[robots_measurement])

      robots_polygon = sort(robots)
      push!(robots_polygon, robots_polygon[1])
      polygon = hcat(attachment_ps[robots_polygon]...)
      @show size(polygon)

      plot3D(polygon[1,:][:], polygon[2,:][:], polygon[3,:][:],
        color = "k", alpha = 1, linestyle = "solid")

      if jj > 1
        old_robots = data_array[ii,jj-1].robots
        if length(robots) > length(old_robots)
          plot_new_point(attachment_ps[setdiff(robots, old_robots)[1]])
        end
      end
      fill3d(hcat(circle_ps...), alpha=0.2)

      xlabel("X")
      ylabel("Y")
      zlabel("M")
      #fig[:axes][1][:get_yaxis]()[:set_visible](false)
      #fig[:axes][1][:get_xaxis]()[:set_visible](false)
      #axis("off")
      xlim(-1.0,1.0)
      ylim(-1.0,1.0)
      zlim(0.0,1.405)
      savefig("$(trial_folder)/belief_$(jj).png", pad_inches=0.01, bbox_inches="tight")
    end
  end
end

pygui(true)

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
function timestring(t)
    if t < 1e3
        unit = "ns"
    elseif t < 1e6
        t /= 1e3
        unit = "μs"
    elseif t < 1e9
        t /= 1e6
        unit = "ms"
    else
        t /= 1e9
        unit = "s"
    end
    return string(round(t, digits=3)) * " " * unit
end

function _frame(
    agg=minimum;
    baseline,
    target,
    benchmarks,
    extra_meta=(
        "size" => :size,
        "complex" => :complex,
        "dim" => :dim_max,
    )
)
    names = collect(intersect(keys(baseline), keys(target)))
    targets = [agg(target[name]) for name in names]
    baselines = [agg(baseline[name]) for name in names]
    ratios = [string(round(r.ratio.time, digits=3)) for r in judge.(targets, baselines)]
    bs = [benchmarks[findfirst(b -> b.name == name, benchmarks)] for name in names]
    metas = [key => get_meta(value).(bs)
        for (key, value) in extra_meta
    ]

    return DataFrame([
        "name" => names,
        metas...,
        "target" => [timestring(t.time) for t in targets],
        "baseline" => [timestring(b.time) for b in baselines],
        "ratio" => ratios,
     ])
end

function save_csv()
    result = BenchmarkTools.load(
        joinpath(@__DIR__, "..", "results" , "all_results.json")
    )[1]
    res(f) = joinpath(@__DIR__, "..", "results", "csv", f)
    CSV.write(
        res("involuted-vs-cohomology.csv"),
        _frame(
            baseline=result["involuted"]["cohomology"],
            target=result["involuted"]["involuted"],
            benchmarks=involuted_benchmarks,
            extra_meta=(
                "size" => :size,
                "dim" => :dim_max,
                "threshold" => :threshold,
                "complex" => :complex
            )
        )
    )
    CSV.write(
        res("involuted-vs-homology.csv"),
        _frame(
            baseline=result["involuted"]["homology"],
            target=result["involuted"]["involuted"],
            benchmarks=involuted_benchmarks,
            extra_meta=(
                "size" => :size,
                "dim" => :dim_max,
                "threshold" => :threshold,
                "complex" => :complex
            )
        )
    )
    CSV.write(
        res("involuted-vs-eirene.csv"),
        _frame(
            baseline=result["eirene"]["eirene"],
            target=result["eirene"]["ripserer"],
            benchmarks=eirene_benchmarks,
            extra_meta=(
                "size" => :size,
                "dim" => :dim_max,
                "threshold" => :threshold,
                "complex" => :complex
            )
        )
    )
    CSV.write(
        res("rips-vs-ripser-sparse.csv"),
        _frame(
            baseline=result["ripser"]["rips"]["ripser"],
            target=result["ripser"]["rips"]["ripserer"]["sparse"],
            benchmarks=rips_benchmarks,
            extra_meta=(
                "size" => :size,
                "dim" => :dim_max,
                "threshold" => :threshold,
            )
        )
    )
    CSV.write(
        res("rips-vs-ripser-dense.csv"),
        _frame(
            baseline=result["ripser"]["rips"]["ripser"],
            target=result["ripser"]["rips"]["ripserer"]["dense"],
            benchmarks=rips_benchmarks,
            extra_meta=(
                "size" => :size,
                "dim" => :dim_max,
                "threshold" => :threshold,
            )
        )
    )
    CSV.write(
        res("cubical-vs-ripser-dense.csv"),
        _frame(
            baseline=result["ripser"]["cubical"]["ripser"],
            target=result["ripser"]["cubical"]["ripserer"],
            benchmarks=cubical_benchmarks,
            extra_meta=(
                "size" => :size,
                "dim" => :dim_max,
                "threshold" => :threshold,
            )
        )
    )
    CSV.write(
        res("alpha-rips-vs-ripser.csv"),
        _frame(
            baseline=result["ripser"]["alpha_rips"]["ripser"],
            target=result["ripser"]["alpha_rips"]["ripserer"],
            benchmarks=alpha_rips_benchmarks,
            extra_meta=(
                "size" => :size,
                "dim" => :dim_max,
                "threshold" => :threshold,
            )
        )
    )
end


text:
using ADCME
using PyPlot
db = Database("hessian.db")

close("all")
plot(npoints, npoints, "k--")

for seed in [2, 23, 233]
    c = execute(db, """
    SELECT nlambda, npoints from eigvals where seed=$seed order by npoints
    """)
    vals = collect(c)
    global npoints = [x[2] for x in vals]
    vals = [x[1] for x in vals]
    plot(npoints, vals, ".-", label = "Seed = $seed")
end
legend()
xlabel("\$||\\mathcal{I}||\$")
ylabel("Positive Eigenvalues")
ylim(0, 15)
grid("on")
minorticks_on()
savefig("hessian_eigenvalue_rank.png")


text:
module TestKaleido
using Test

@testset "$file" for file in sort([file for file in readdir(@__DIR__) if
                                   match(r"^test_.*\.jl$", file) !== nothing])
    include(file)
end

using Aqua
using Kaleido
#=
Aqua.test_all(Kaleido)
=#
@testset "Aqua" begin
    @testset "Method ambiguity" begin
        # Not including `Base` due to
        # https://github.com/JuliaCollections/DataStructures.jl/pull/511
        Aqua.test_ambiguities(Kaleido)
    end
    @testset "Unbound type parameters" begin
        Aqua.test_unbound_args(Kaleido)
    end
    @testset "Undefined exports" begin
        Aqua.test_undefined_exports(Kaleido)
    end
end

end  # module


text:
# Now, you need to figure out how to pilot this thing.
#
# It seems like the submarine can take a series of commands like forward 1,
# down 2, or up 3:
#
# forward X increases the horizontal position by X units.
# down X increases the depth by X units.
# up X decreases the depth by X units.
#
# Note that since you're on a submarine, down and up affect your depth, and so
# they have the opposite result of what you might expect.
#
# The submarine seems to already have a planned course (your puzzle input). You
# should probably figure out where it's going. For example:
#
# forward 5
# down 5
# forward 8
# up 3
# down 8
# forward 2
#
# Your horizontal position and depth both start at 0. The steps above would
# then modify them as follows:
#
# forward 5 adds 5 to your horizontal position, a total of 5.
# down 5 adds 5 to your depth, resulting in a value of 5.
# forward 8 adds 8 to your horizontal position, a total of 13.
# up 3 decreases your depth by 3, resulting in a value of 2.
# down 8 adds 8 to your depth, resulting in a value of 10.
# forward 2 adds 2 to your horizontal position, a total of 15.
#
# After following these instructions, you would have a horizontal position of
# 15 and a depth of 10. (Multiplying these together produces 150.)
#
# Calculate the horizontal position and depth you would have after following
# the planned course. What do you get if you multiply your final horizontal
# position by your final depth?

dirs = map(eachline("input.txt")) do l
  (dir, amt) = split(l, ' ')
  amt = parse(Int, amt)
  if dir == "forward"
    (0, amt)
  elseif dir == "up"
    (-amt, 0)
  else
    (amt, 0)
  end
end
(depth, horz) = reduce(.+, dirs)
println("Result: ", depth * horz)


text:
using ARCTR
using LinearAlgebra
using Test
using Printf
using NLPModels
using Stopping

using ADNLPModels
using OptimizationProblems.ADNLPProblems

include("testLDLt.jl")
# include("testAbs.jl") #  MethodError: no method matching ARCSpectral(::ADNLPModel{Float64, Vector{Float64}}; verbose=true)
# include("testMA57_0.jl")
# include("testMA97.jl")

global nbsolver = 0
for solver in ALL_solvers
    global nbsolver += 1
    nlp = extrosnb(n = 2)
    nlpstop = NLPStopping(nlp)
    println(nbsolver, "  ", solver)
    eval(solver)(nlpstop, verbose = true)
    final_nlp_at_x, optimal = nlpstop.current_state, nlpstop.meta.optimal
    @test optimal
    reset!(nlp)
    stats = eval(solver)(nlp, verbose = false)
    @test stats.status == :first_order
    reset!(nlp)
end


text:
module ProtoBufTestApis
using ProtoBuf
using ..Test
import ProtoBuf.meta

print_hdr(tname) = println("testing $tname...")

mutable struct TestType <: ProtoType
    a::AbstractString
    b::Bool
    TestType() = (o=new(); fillunset(o); o)
end #type TestType
meta(t::Type{TestType}) = meta(t, Symbol[:a], Int[], Dict{Symbol,Any}())

function test_apis()
    t = TestType()

    @test !has_field(t, :a)
    @test !has_field(t, :b)

    @test false == try get_field(t, :a); true; catch; false; end

    t.b = true
    @test has_field(t, :b)
    @test (get_field(t, :b) == true)

    @test !isinitialized(t)
    t.a = "hello world"
    @test isinitialized(t)
    @test (get_field(t, :a) ==  "hello world")

    clear(t, :b)
    @test isinitialized(t)
    clear(t)
    @test !isinitialized(t)

    t = protobuild(TestType, Dict(:a => "hello", :b => false))
    @test t.a == "hello"
    @test t.b == false
end

function test_deepcopy()
    ts = ProtoBuf.google.protobuf.Timestamp()
    ts.seconds = 123
    @test !has_field(ts, :nanos)
    ts2 = deepcopy(ts)
    @test !has_field(ts2, :nanos)
end

end # module ProtoBufTestApis

print_hdr("utility api methods")
ProtoBufTestApis.test_apis()
ProtoBufTestApis.test_deepcopy()