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)
    if vors == "volume"
        # args = (var, eid, fid, grid_data, geo_factors, fv_info, t, dt)
        code =
"var =       args[1]; # variable list
eid =       args[2]; # element ID
fid =       args[3]; # NA
grid =      args[4]; # grid data struct
geo_facs =  args[5]; # geometric factors
fv_data =   args[6]; # FV specific info
refel =     args[7]; # reference element
time =      args[8]; # time
dt =        args[9]; # dt
"
    else # surface
        # args = (var, eid, fid, els, grid_data, geo_factors, fv_info, t, dt)
        code =
"var =       args[1]; # variable list
eid =       args[2]; # element this flux is applied to
fid =       args[3]; # face ID
els =       args[4]; # neighborhood element IDs [[left], [right]]
grid =      args[5]; # grid data struct
geo_facs =  args[6]; # geometric factors
fv_data =   args[7]; # FV specific info
refel =     args[8]; # reference element
time =      args[9]; # time
dt =        args[10]; # dt
"
    end
    return code;
end

# Allocate, compute, or fetch all needed values
function prepare_needed_values_fv_julia(entities, var, lorr, vors)
    # If there are no entities, do nothing
    if length(entities) < 1
        return "";
    end

    # Options
    average_coefficients = false; # true to integrate and average coefficients over cells and faces.

    # 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;
    fv_order = fv_info.fluxOrder;

    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_entity_name(entities[i]);
        if is_unknown_var(entities[i], var) && lorr == LHS
            # The unknown variable will not be included as it is being solved for,
            # but it should be checked to make sure it is just a linear dependence. 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;
                elseif 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] = (Finch.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"
                else
                    nodesymbol = "nodex"
                end
                coef_index = get_coef_index(entities[i]);

                if vors == "volume"
                    if average_coefficients
                        code *= cname * " = zeros(refel.Np);\n";
                        # code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Finch.genfunctions["*string(cval)*"]).func" * cargs * " end\n";
                        code *= "for coefi = 1:refel.Np " * cname * "[coefi] = evaluate_coefficient(Finch.coefficients["*string(coef_index)*"], "*string(entities[i].index)*", "*nodesymbol*"[:,coefi], time, eid, fid) 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;
                        piece_needed[2] = true;
                    else # just evaluate coefficient at center
                        code *= cname * " = evaluate_coefficient(Finch.coefficients["*string(coef_index)*"], "*string(entities[i].index)*", fv_data.cellCenters[:,eid], time, eid, fid);\n";
                    end

                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] = (Finch.genfunctions["*string(cval)*"]).func" * cargs * " end\n";
                        code *= "for coefi = 1:refel.Np " * cname * "[coefi] = evaluate_coefficient(Finch.coefficients["*string(coef_index)*"], "*string(entities[i].index)*", "*nodesymbol*"[:,coefi], time, eid, fid) 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.";
                        # 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
                        need_deriv_matrix = true;
                        piece_needed[[4,5,8]] .= true;

                    else # no derivatives, only need surface nodes
                        if average_coefficients
                            code *= cname * " = zeros(refel.Nfp[frefelind[1]]);\n";
                            # code *= "for coefi = 1:refel.Nfp[frefelind[1]] " * cname * "[coefi] = (Finch.genfunctions["*string(cval)*"]).func" * cargs * " end\n";
                            code *= "for coefi = 1:refel.Nfp[frefelind[1]] " * cname * "[coefi] = evaluate_coefficient(Finch.coefficients["*string(coef_index)*"], "*string(entities[i].index)*", "*nodesymbol*"[:,coefi], time, eid, fid) end\n";
                            piece_needed[5] = true;
                            piece_needed[6] = true;

                            # 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
                        else # just evaluate coefficient at center
                            code *= cname * " = evaluate_coefficient(Finch.coefficients["*string(coef_index)*"], "*string(entities[i].index)*", fv_data.faceCenters[:,fid], time, eid, fid);\n";
                        end
                    end
                end

            elseif ctype == 3 # a known variable value
                piece_needed[1] = true;
                # This generates something like: coef_u_1 = copy((Finch.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;
                #     elseif occursin("CENTRAL", entities[i].flags[flagi])
                #         cellside = 3;
                #     elseif occursin("NEIGHBORHOOD", entities[i].flags[flagi])
                #         cellside = 4;
                #     end
                # end
                cellside = get_face_side_info(entities[i]); # in code_generator_utils.jl
                if vors == "surface"
                    l2gsymbol = "els1"
                else
                    l2gsymbol = "el"
                end
                if cellside == 1
                    l2gsymbol = "els1"
                elseif cellside == 2
                    l2gsymbol = "els2"
                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
                        if variables[cval].discretization == FV
                            code *= cname * " = Finch.variables["*string(cval)*"].values["*indstr*", "*l2gsymbol*"];\n";
                        else
                            code *= cname * " = sum(Finch.variables["*string(cval)*"].values["*indstr*", Finch.grid_data.loc2glb[:,eid]]) / size(Finch.grid_data.loc2glb,1);\n";
                        end
                    end

                else # surface
                    if variables[cval].discretization == FV
                        # Need to reconstruct it using neighboring cells
                        if fv_order == 1
                            if length(entities[i].derivs) > 0
                                code *= cname * " = Finch.variables["*string(cval)*"].values["*indstr*", els2] - Finch.variables["*string(cval)*"].values["*indstr*", els1];\n";
                                code *= cname * " = (els1 != els2 && 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 * (Finch.variables["*string(cval)*"].values["*indstr*", els1] + Finch.variables["*string(cval)*"].values["*indstr*", els2]);\n";
                                elseif cellside < 3
                                    code *= cname * " = Finch.variables["*string(cval)*"].values["*indstr*", "*l2gsymbol*"];\n";
                                elseif cellside == 3 # central
                                    # same as 0 for this case
                                    code *= cname * " = 0.5 * (Finch.variables["*string(cval)*"].values["*indstr*", els1] + Finch.variables["*string(cval)*"].values["*indstr*", els2]);\n";
                                elseif cellside == 4 # neighborhood
                                    # This is a special case only for callback functions.
                                    # Rather than representing a value, it constructs a Neighborhood object to be passed.
                                    code *= cname * " = Finch.Neighborhood(els, cellx, [Finch.variables["*string(cval)*"].values["*indstr*", els[1]], Finch.variables["*string(cval)*"].values["*indstr*", els[2]]]);\n";
                                    piece_needed[4] = true; # cellx
                                end
                            end

                        else # order > 1
                            # collect the cell info
                            piece_needed[6] = true; # face center coords
                            piece_needed[4] = true; # cell centers
                            if length(entities[i].derivs) > 0
                                code *= cname * " = Finch.variables["*string(cval)*"].values["*indstr*", els[2][1]] - Finch.variables["*string(cval)*"].values["*indstr*", els[1][1]];\n";
                                code *= cname * " = (els[1][1] != els[2][1] && 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 || cellside == 3
                                    # Use the full neighborhood
                                    code *= cname * " = Finch.FV_reconstruct_value([cellx[1] cellx[2]], [Finch.variables["*string(cval)*"].values["*indstr*", els[1]] ; Finch.variables["*string(cval)*"].values["*indstr*", els[2]]], facex);\n";
                                elseif cellside == 1 # left
                                    code *= cname * " = Finch.FV_reconstruct_value_left_right(cellx[1], cellx[2], Finch.variables["*string(cval)*"].values["*indstr*", els[1]], Finch.variables["*string(cval)*"].values["*indstr*", els[2]], facex, limiter=\"vanleer\")[1];\n";
                                elseif cellside == 2 # right
                                    code *= cname * " = Finch.FV_reconstruct_value_left_right(cellx[1], cellx[2], Finch.variables["*string(cval)*"].values["*indstr*", els[1]], Finch.variables["*string(cval)*"].values["*indstr*", els[2]], facex, limiter=\"vanleer\")[2];\n";
                                elseif cellside == 2 # right
                                    code *= cname * " = Finch.Neighborhood(els, cellx, [Finch.variables["*string(cval)*"].values["*indstr*", els[1]], Finch.variables["*string(cval)*"].values["*indstr*", els[2]]]);\n";
                                    piece_needed[4] = true; # cellx
                                end
                            end
                        end
                    else # A nodal variable that has a value specified at the surface
                        if length(entities[i].derivs) > 0
                            code *= cname * "_TMP1 = sum(Finch.variables["*string(cval)*"].values["*indstr*", Finch.grid_data.loc2glb[:,els1]]) / size(Finch.grid_data.loc2glb,1);\n";
                            code *= cname * "_TMP2 = sum(Finch.variables["*string(cval)*"].values["*indstr*", Finch.grid_data.loc2glb[:,els2]]) / size(Finch.grid_data.loc2glb,1);\n";
                            code *= cname * " = "*cname*"_TMP2 - "*cname*"_TMP1;\n";
                            code *= cname * " = (els1 != els2 && 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 < 4
                                code *= cname * " = sum(Finch.variables["*string(cval)*"].values["*indstr*", Finch.grid_data.face2glb[:,1,fid]]) / size(Finch.grid_data.face2glb, 1);\n";
                            elseif cellside == 4 # neighborhood
                                # This is a special case only for callback functions.
                                # Rather than representing a value, it constructs a Neighborhood object to be passed.
                                code *= cname * " = sum(Finch.variables["*string(cval)*"].values["*indstr*", Finch.grid_data.face2glb[:,1,fid]]) / size(Finch.grid_data.face2glb, 1);\n";
                                code *= cname * " = Finch.Neighborhood(els, cellx, ["*cname*", "*cname*"]);\n";
                                piece_needed[4] = true; # cellx
                            end
                        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] = (Finch.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"
                else
                    nodesymbol = "nodex"
                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"
                    if average_coefficients
                        code *= cname * " = zeros(refel.Np);\n";
                        # code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Finch.coefficients["*string(cval)*"]).value["*indstr*"].func" * cargs * " end\n";
                        code *= "for coefi = 1:refel.Np " * cname * "[coefi] = evaluate_coefficient(Finch.coefficients["*string(cval)*"], ["*indstr*"], "*nodesymbol*"[:,coefi], time, eid, fid) 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;
                        piece_needed[2] = true;
                    else # just evaluate coefficient at center
                        code *= cname * " = evaluate_coefficient(Finch.coefficients["*string(cval)*"], ["*indstr*"], fv_data.cellCenters[:,eid], time, eid, fid);\n";
                    end

                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] = (Finch.coefficients["*string(cval)*"]).value["*indstr*"].func" * cargs * " end\n";
                        code *= "for coefi = 1:refel.Np " * cname * "[coefi] = evaluate_coefficient(Finch.coefficients["*string(cval)*"], ["*indstr*"], "*nodesymbol*"[:,coefi], time, eid, fid) 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.";
                        # 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
                        need_deriv_matrix = true;
                        piece_needed[[4,5,8]] .= true;

                    else # no derivatives, only need surface nodes
                        if average_coefficients
                            code *= cname * " = zeros(refel.Nfp[frefelind[1]]);\n";
                            # code *= "for coefi = 1:refel.Nfp[frefelind[1]] " * cname * "[coefi] = (Finch.coefficients["*string(cval)*"]).value["*indstr*"].func" * cargs * " end\n";
                            code *= "for coefi = 1:refel.Nfp[frefelind[1]] " * cname * "[coefi] = evaluate_coefficient(Finch.coefficients["*string(cval)*"], ["*indstr*"], "*nodesymbol*"[:,coefi], time, eid, fid) end\n";
                            # 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
                            piece_needed[5] = true;
                            piece_needed[6] = true;
                        else # just evaluate coefficient at center
                            code *= cname * " = evaluate_coefficient(Finch.coefficients["*string(cval)*"], ["*indstr*"], fv_data.faceCenters[:,fid], time, eid, fid);\n";
                        end
                    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_facs.detJ[eid]; # geometric factors\n"
        end
        if piece_needed[5]
            needed_pieces *= "J = geo_facs.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 length(els[2]) == 0 # This is a boundary face. For now, just compute as if neighbor is identical. BCs handled later.
    els[2] = [eid];
end
";
        if fv_order == 1
            needed_pieces *= "els1 = els[1][1]; els2 = els[2][1]; # For first order only.\n";
        end
        # els, nodex, loc2glb, cellx, frefelind, facex, face2glb, normal, face_detJ, area, vol_J
        if piece_needed[5]
            needed_pieces *= "frefelind = [grid.faceRefelInd[1,fid], grid.faceRefelInd[2,fid]]; # refel based index of face in both elements\n"
        end
        if piece_needed[8] || true # probably always need
            needed_pieces *= "normal = grid.facenormals[:, fid];\n"
        end
        if piece_needed[2] || piece_needed[3]
            needed_pieces *= "loc2glb = [grid.loc2glb[:,els[1][1]], grid.loc2glb[:, els[2][1]]]; # 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[:, els[1]], fv_data.cellCenters[:, els[2]]]; # 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 = fv_data.faceCenters[:, fid];  # face node coordinates\n"
            needed_pieces *= "facex = grid.allnodes[:, grid.face2glb[:,1,fid]];  # 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[els[1][1]], geo_facs.J[els[2][1]]);\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][:,1] - cellx[1][:,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,
    how_many_terms = 0;
    for i=1:length(terms)
        how_many_terms = how_many_terms + length(terms[i]);
    end
    if how_many_terms < 1
        if lorr == LHS && vors == "surface"
            code = "return [0   0] # There were no terms to compute here";
        else
            code = "return [0] # There were no terms to compute here";
        end
        return code;
    end

    # Allocate the vector or matrix to be returned if needed
    if dofsper > 1
        if lorr == RHS
            code *= "result = zeros("*string(dofsper)*"); # Allocate for returned values.\n"
        else
            if vors == "volume"
                code *= "result = zeros("*string(dofsper)*", "*string(dofsper)*"); # Allocate for returned values.\n"
            else
                code *= "result = zeros(2 * "*string(dofsper)*", "*string(dofsper)*"); # Allocate for returned values.\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, vors);

                        # 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, vors);

                    # 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 result 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 *= "result["*string(emi)*"] = " * subvector[emi] * "\n";
                end
            end
            if vors == "surface"
                code *= "
if els[2][1] == eid && els[1][1] != els[2][1]
    result = -result; # Since this flux is applied to element eid, make sure it's going in the right direction.
end
return result;
"
            end

        end

    else # one dof
        terms = terms[1];
        if lorr == LHS
            if vors == "volume"
                result = "0";
            else
                result1 = "0";
                result2 = "0";
            end

        else
            result = "0";
        end

        #process each term
        first_term1 = true;
        first_term2 = true;
        for i=1:length(terms)
            (term_result1, term_result2, var_ind) = generate_term_calculation_fv_julia(terms[i], var, lorr, vors);
            if lorr == LHS && vors == "surface"
                # There are two results, one for each side of the face
                if !(term_result1 == "0" || term_result1 == "")
                    if first_term1
                        result1 = term_result1;
                        first_term1 = false;
                    else
                        result1 *= " .+ " * term_result1;
                    end
                end
                if !(term_result2 == "0" || term_result2 == "")
                    if first_term2
                        result2 = term_result2;
                        first_term2 = false;
                    else
                        result2 *= " .+ " * term_result2;
                    end
                end
            else # only term_result1 is used
                if !(term_result1 == "0" || term_result1 == "")
                    if first_term1
                        result = term_result1;
                        first_term1 = false;
                    else
                        result *= " .+ " * term_result1;
                    end
                end
            end

        end # term loop

        if lorr == LHS && vors == "surface" # LHS surface returns two pieces per dof
            code *= "result = [" * result1 * "    " * result2 * "];\n";

        elseif vors == "surface" # RHS surface
            code *= "result = [" * result * "];\n";
            code *= "
if els[2][1] == eid && els[1][1] != els[2][1]
    result = -result; # Since this flux is applied to element eid, make sure it's going in the right direction.
end
return result;
"
        else # volume (LHS and RHS)
            code *= "result = [" * result * "];\n";
        end
    end

    return code;
end

function generate_term_calculation_fv_julia(term, var, lorr, vors)
    # This will return 3 things: for volume terms, the first is the result.
    # For surface terms, there is one for each side.
    # The third is the var_ind
    result1 = "";
    result2 = "";
    # Note: separate_factors return test and trial info, but FV will not have any test functions.
    # var_part refers to the unknown variable part.
    # example:
    # 0.1*D1__u_1*_FACENORMAL1_1  ->  0.1 * normal[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);
        which_side = get_face_side_info(var_part); # 0=none, 1/2=side 1/2, 3=average, 4=neighborhood

        if coef_part === nothing # Only an unknown variable
            # handle negative sign
            if typeof(var_part) == Expr && (var_part.args[1] === :.- || var_part.args[1] === :-)
                if vors == "surface"
                    if which_side == 0 || which_side == 3 # average them
                        result1 = "-0.5";
                        result2 = "-0.5";
                    elseif which_side == 1
                        result1 = "-1";
                    elseif which_side == 2
                        result2 = "-1";
                    else
                        # neighborhoods are not supported for implicit steppers
                    end
                else
                    result1 = "-1";
                end

            else
                if vors == "surface"
                    if which_side == 0 || which_side == 3 # average them
                        result1 = "0.5";
                        result2 = "0.5";
                    elseif which_side == 1
                        result1 = "1";
                    elseif which_side == 2
                        result2 = "1";
                    else
                        # neighborhoods are not supported for implicit steppers
                    end
                else
                    result1 = "1";
                end
            end

        else
            if vors == "surface"
                # This is tricky. First make two copies of the coef_part.
                # Since the coef_part could actually contain a variable hidden in an expression like a conditional,
                # Set the correct side to 1 and the other side to 0, or if no side, set both to 0.5
                # Also make side 2 negative... because I need to make upwinding work.
                # This is a bad way to rig this up, but I'm out of ideas here.
                coef_side1 = copy(coef_part);
                coef_side2 = copy(coef_part);
                coef_side1 = replace_lhs_surface_var_entities(coef_side1, var, 1);
                coef_side2 = replace_lhs_surface_var_entities(coef_side2, var, 2);
                if which_side == 0 || which_side == 3 # average them
                    result1 = string(replace_entities_with_symbols(coef_side1));
                    result2 = string(replace_entities_with_symbols(coef_side2));
                elseif which_side == 1
                    result1 = string(replace_entities_with_symbols(coef_side1));
                elseif which_side == 2
                    result2 = string(replace_entities_with_symbols(coef_side2));
                else
                    # neighborhoods are not supported for implicit steppers
                end
            else
                result1 = string(replace_entities_with_symbols(coef_part));
            end
        end

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

    return (result1, result2, var_ind);
end

# A special function for handling LHS surface parts.
# Swaps variable entities with a number depending on side.
# negate will make side 2 negative.
function replace_lhs_surface_var_entities(ex, var, side, negate=true)
    if typeof(ex) == Expr
        for i=1:length(ex.args)
            ex.args[i] = replace_lhs_surface_var_entities(ex.args[i], var, side, negate);
        end
    elseif typeof(ex) == SymEntity
        if is_unknown_var(ex, var)
            which_side = get_face_side_info(ex);
            if which_side == 0 || which_side == 3
                ex = 0.5;
            elseif which_side == side
                ex = 1;
            else
                ex = 0;
            end
            if side == 2 && negate
                ex = -ex;
            end
        end
    else
        # number?
    end
    return ex;
end

# Generate the assembly loop structures and insert the content
function generate_assembly_loop_fv_julia(var, indices, parallel_type)
    # 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];
        parallel_type = ["none", parallel_type];
    end

    code =
"# Label things that were allocated externally
if is_explicit
    sourcevec = allocated_vecs[1];
    fluxvec = allocated_vecs[2];
    facefluxvec = allocated_vecs[3];
    face_done = allocated_vecs[4];
else
    lhsmatV = allocated_vecs[3];
    sourcevec = allocated_vecs[4]; # In the implicit case, these vectors are for the RHS
    fluxvec = allocated_vecs[5];
    facefluxvec = allocated_vecs[6];
    face_done = allocated_vecs[7];

    for i=1:length(lhsmatV)
        lhsmatV[i] = 0;
    end
end

nel = fv_grid.nel_owned;
dofs_squared = dofs_per_node*dofs_per_node;
"
    # generate the loop structures
    code *= "# Loops\n"
    loop_start = ""
    loop_end = ""
    for i=1:length(indices)
        if indices[i] == "elements" || indices[i] == "cells"
            loop_start *= "for eid in elemental_order\n";
            loop_start *= "
for i=1:length(face_done)
    face_done[i] = 0; # Reset face_done in elemental loop
end\n";
            loop_end *= "end # loop for elements\n";
        else
            indexer_index = 0;
            for j=1:length(indexers)
                if indices[i].symbol === indexers[j].symbol
                    indexer_index += j;
                    break;
                end
            end
            loop_start *= "for indexing_variable_"*string(indices[i].symbol)*" in Finch.indexers["*string(indexer_index)*"].range\n";
            loop_start *= "    Finch.indexers["*string(indexer_index)*"].value = indexing_variable_"*string(indices[i].symbol)*";\n";
            loop_end *= "end # loop for "*string(indices[i].symbol)*"\n";
        end
    end

    # The index offset depends on the order in which indexers were used during variable creation.
    # They must all be the same for now.
    ind_offset = ""
    prev_ind = ""
    if typeof(var) <: Array
        var_indices = var[1].indexer;
    else
        var_indices = var.indexer;
    end
    if !(typeof(var_indices) <: Array)
        var_indices = [var_indices];
    end
    for i=1:length(var_indices)
        if length(prev_ind) == 0
            ind_offset *= "(indexing_variable_"*string(var_indices[i].symbol);
            prev_ind = string(length(var_indices[i].range));
        else
            ind_offset *= " + "*prev_ind*" * (indexing_variable_"*string(var_indices[i].symbol)*" - 1";
            prev_ind = string(length(var_indices[i].range));
        end
    end
    for i=1:length(var_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"
    index_offset = "dofs_per_node * (eid - 1) + ("*ind_offset*" - 1) * dofs_per_loop + 1";
    face_index_offset = "dofs_per_node * (fid - 1) + ("*ind_offset*" - 1) * dofs_per_loop + 1";

    code *= loop_start;

    # insert the content
    code *= "
    # determine index in global vector
    index_offset = "*index_offset*";
    first_ind = index_offset; # First global index for this loop
    last_ind = index_offset + dofs_per_loop-1; # last global index

    # Zero the result vectors for this element
    sourcevec[first_ind:last_ind] .= 0;
    fluxvec[first_ind:last_ind] .= 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*");

        # Add to global source vector
        sourcevec[first_ind:last_ind] .= source;
    end

    # Compute LHS volume integral = (block)diagonal components
    if !(source_lhs === nothing) && !is_explicit
        sourceargs = (var, eid, 0, fv_grid, fv_geo_factors, fv_info, refel, t, dt);
        source = source_lhs.func(sourceargs; "*index_args*");
        # Add to global matrix
        for di = 1:dofs_per_loop
            first = (eid-1)*dofs_squared + "*ind_offset*"; # This needs to change for multi dof TODO
            last = first + dofs_per_loop - 1;
            lhsmatV[first:last] = source[((di-1)*dofs_per_loop + 1):(di*dofs_per_loop)];
        end
    end

    ##### Flux integrated over the faces #####
    # Loop over this element's faces.
    for i=1:refel.Nfaces
        fid = grid_data.element2face[i, eid];
        do_face_here = (face_done[fid] < dofs_per_node);
        if do_face_here
            face_done[fid] += dofs_per_loop;
        end
        # Only one element on either side is available here. For more use parent/child version.
        (leftel, rightel) = grid_data.face2element[:,fid];
        if rightel == 0
            neighborhood = [[leftel],[]];
        else
            neighborhood = [[leftel],[rightel]];
        end

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

        if !(flux_rhs === nothing)
            if do_face_here

                fluxargs = (var, eid, fid, neighborhood, 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

        # LHS surface integral
        if !(flux_lhs === nothing) && !is_explicit
            # These update the (block)diagonal for eid and one (block of)off diagonal for neighborID
            if do_face_here
                face_done[fid] = 1; # Possible race condition.
                neighborID = (rightel==eid ? leftel : rightel);
                eid_flux_index = (leftel==eid ? 1 : 2);
                neighbor_flux_index = (leftel==eid ? 2 : 1);
                nfirst_ind = (neighborID-1)*dofs_per_node + "*ind_offset*";
                nlast_ind = neighborID*dofs_per_node + "*ind_offset*" - 1;

                fluxargs = (var, eid, fid, neighborhood, fv_grid, fv_geo_factors, fv_info, refel, t, dt);
                flux = flux_lhs.func(fluxargs; "*index_args*"); println(flux);
                flux = flux_lhs.func(fluxargs; "*index_args*") .* fv_geo_factors.area[fid];
                # insert the matrix elements
                for di = 1:dofs_per_loop
                    # The eid components
                    first = nel * dofs_squared + (fid-1)*dofs_squared*4 + "*ind_offset*";
                    last = first + dofs_per_loop - 1;
                    lhsmatV[first:last] = flux[((di-1)*dofs_per_loop + 1):(di*dofs_per_loop), eid_flux_index] ./ fv_geo_factors.volume[eid];

                    # The neighborID components
                    if neighborID > 0
                        # first = nel * dofs_squared + (fid-1)*dofs_squared*4 + dofs_squared + "*ind_offset*";
                        # last = first + dofs_per_loop - 1;
                        first += dofs_squared;
                        last += dofs_squared;
                        lhsmatV[first:last] = flux[((di-1)*dofs_per_loop + 1):(di*dofs_per_loop), neighbor_flux_index] ./ fv_geo_factors.volume[eid];
                    end
                end

                # While we're here, set the components for the neighborID rows as well.
                if neighborID > 0
                    for di = 1:dofs_per_loop
                        # The eid components
                        first = nel * dofs_squared + (fid-1)*dofs_squared*4 + dofs_squared*2 + "*ind_offset*";
                        last = first + dofs_per_loop - 1;
                        lhsmatV[first:last] = -flux[((di-1)*dofs_per_loop + 1):(di*dofs_per_loop), eid_flux_index] ./ fv_geo_factors.volume[neighborID];

                        # The neighborID components
                        first += dofs_squared;
                        last += dofs_squared;

                        lhsmatV[first:last] = -flux[((di-1)*dofs_per_loop + 1):(di*dofs_per_loop), neighbor_flux_index] ./ fv_geo_factors.volume[neighborID];
                    end
                end

            else
                # This flux has either been computed or is being computed by another thread.
                # Everything was set there.
            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 = "*ind_offset*";
                    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], fid, t) .* geo_factors.area[fid];

                        fluxvec[index_offset + dofind-1] += (bflux - facefluxvec[face_index_offset + dofind-1]) ./ geo_factors.volume[eid];
                        facefluxvec[face_index_offset + dofind-1] = bflux;
                    elseif prob.bc_type[var[vi].index, fbid] == DIRICHLET
                        # Set variable array and handle after the face loop
                        var[vi].values[compo,eid] = FVSolver.evaluate_bc(prob.bc_func[var[vi].index, fbid][compo], eid, fid, t);
                    else
                        printerr(\"Unsupported boundary condition type: \"*prob.bc_type[var[vi].index, fbid]);
                    end
                end

            else # one variable
                d="*ind_offset*";
                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], fid, t) .* geo_factors.area[fid];

                    fluxvec[index_offset] += (bflux - facefluxvec[face_index_offset]) ./ geo_factors.volume[eid];
                    facefluxvec[face_index_offset] = bflux;
                elseif prob.bc_type[var.index, fbid] == DIRICHLET
                    # Set variable array and handle after the face loop
                    var.values[d,eid] = FVSolver.evaluate_bc(prob.bc_func[var.index, fbid][d], eid, fid, t);
                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

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using LinearAlgebra

# If i, i+1 are conjugate pair, then they need to be either both in I or both not in I.
# If one of them is in I and the other is not then LAPACK will consider that both of them are in I.
function conditionnumber(sf::Schur, I)
    n = length(sf.values)
    select = zeros(LinearAlgebra.BlasInt, n)
    for i in I
        select[i] = 1
    end
    LinearAlgebra.LAPACK.trsen!('E', 'N', select, copy(sf.T), copy(sf.Z))[4]
end

# Manocha, D. & Demmel, J. Algorithms for intersecting parametric and algebraic curves II: multiple intersections
# Graphical Models and Image Processing, Elsevier, 1995, 57, 81-100
function clusterordschur(M::AbstractMatrix{<:Real}, ɛ)
    if isempty(M)
        # See bug JuliaLang/julia#...
        return Matrix{float(eltype(M))}(undef, 0, 0), Vector{Int}[]
    else
        _clusterordschur(M, ɛ)
    end
end
function _clusterordschur(M::AbstractMatrix{<:Real}, ɛ)
    # M = Z * T * Z' and "values" gives the eigenvalues
    if VERSION >= v"0.7-"
        sf = schur(M)
        Z = sf.Z
        v = sf.values
    else
        sf = schurfact(M)
        Z = sf[:Z]
        v = sf[:values]
    end
    # documentation says that the error on the eigenvalues is ɛ * norm(T) / conditionnumber
    nT = norm(sf.T)

    _atol(I) = ɛ * nT / conditionnumber(sf, I)

    A = typeof(_atol(1))
    V = real(eltype(v))

    ONE = one(one(V) / one(A))

    clusters = Vector{Int}[]
    λ = V[]
    atol = A[]
    # conditionnumber requires that conjugate pair need to be treated together so we first need to handle them
    # If they are in the same cluster then pair them, otherwise it is complex solution so we reject them
    i = 1
    while i <= lastindex(v)
        if isreal(v[i])
            push!(clusters, [i])
            push!(λ, v[i])
            push!(atol, _atol(i))
            i += 1
        else
            @assert i < lastindex(v) && !isreal(v[i+1])
            pairatol = _atol([i, i+1])
            if abs(v[i] - v[i+1]) / pairatol < ONE
                # Pair conjugate pairs into a real eigenvalue
                push!(clusters, [i, i+1])
                push!(λ, real((v[i] + v[i+1]) / 2)) # The imaginary part should be zero anyway
                push!(atol, pairatol)
            end
            i += 2
        end
    end
    σ = sortperm(λ)

    # For eigenvalues not clustered yet, their eigenvalues is quite large.
    # Therefore, if we cluster all i, j close enough at once we migth cluster too much
    # The technique used here is to cluster only the closest pair.
    # Once they are matched, a new atol is computed and if the cluster is complete,
    # this atol will be small which will avoid addition of new eigenvalues.
    while true
        I = 0
        J = 0
        best = ONE
        for i in eachindex(clusters)
            for j in 1:(i-1)
                d = abs(λ[i] - λ[j]) / min(atol[i], atol[j])
                if d < best
                    I = i
                    J = j
                    best = d
                end
            end
        end
        if best < ONE
            # merge I with J
            nI = length(clusters[I])
            nJ = length(clusters[J])
            λ[I] = (λ[I] * nI + λ[J] * nJ) / (nI + nJ)
            append!(clusters[I], clusters[J])
            atol[I] = _atol(clusters[I])
            deleteat!(λ, J)
            deleteat!(clusters, J)
            deleteat!(atol, J)
        else
            break
        end
    end
    Z, clusters
end

ext: jl


text:
using BenchmarkTools

function day13(printoutput=true)
    dots = Set{Tuple{Int,Int}}()
    folds = Tuple{Char,Int}[]
    for line in filter(!isempty, strip.(readlines("AoCdata/AoC_2021_day13.txt")))
        left, right = split(line, r"[ =,]")[(end - 1):end]
        if left in ["x", "y"]
            push!(folds, (first(left), parse(Int, right)))
        else
            push!(dots, (parse(Int, left) + 1, parse(Int, right) + 1))
        end
    end

    for (i, fold) in enumerate(folds)
        if first(fold) == 'y'  # fold in half with a horzontal fold line, lifting bottom edge up and folding down toward top
            for d in dots
                if last(d) > last(fold)
                    delete!(dots, d)
                    push!(dots, (first(d), 2 * last(fold) - last(d) + 2))
                end
            end
        else           # fold in half with a vertical fold line, lifting right edge up and to left
            for d in dots
                if first(d) > last(fold)
                    delete!(dots, d)
                    push!(dots, (2 * last(fold) - first(d) + 2, last(d)))
                end
            end
        end
        printoutput && i == 1 && println("Part 1: ", length(dots))
    end
    mat = fill(". ", maximum([first(d) for d in dots]), maximum([last(d) for d in dots]))
    for d in CartesianIndices(mat)
        Tuple(d) in dots && (mat[d] = "# ")
    end
    return printoutput &&
           println("Part 2 graphic:\n", join([prod(line) for line in eachcol(mat)], "\n"))
end

day13()

@btime day13(false)

#=
Part 1: 724
Part 2 graphic:
. # # . . # # # . . . . # # . # # # . . # # # # . # # # . . # . . # . # . . .
# . . # . # . . # . . . . # . # . . # . # . . . . # . . # . # . . # . # . . .
# . . . . # . . # . . . . # . # # # . . # # # . . # . . # . # . . # . # . . .
# . . . . # # # . . . . . # . # . . # . # . . . . # # # . . # . . # . # . . .
# . . # . # . . . . # . . # . # . . # . # . . . . # . # . . # . . # . # . . .
. # # . . # . . . . . # # . . # # # . . # # # # . # . . # . . # # . . # # # #
  1.034 ms (5295 allocations: 583.92 KiB)
=#


text:
module SparseMatrices

# package code goes here
include("csr.jl")

end # module


text:
include("constants.jl")
include("VelocityScale.jl")
include("holtslag_psim.jl")
include("holtslag_psis.jl")
include("holtslag_local.jl")

function holtslag(θv,U,V,LH,SH,ustar,ρA,ZF)
# a vertical diffusion scheme based on Holtslag and Boville (1993)
# Includes enhenced diffusion below the PBL
#θv-virtual potential temperature

Rcrit=0.5;
b=8.5
eps=0.1;
minshear=0.003
d=10

Z=ZF[2:end];
DZ=diff(ZF,dims=1);
ZM=Z-DZ./2;
SZM=length(ZM);

ZD=[Z[1]:0.01*(Z[1]):Z[end];]

#wθv-surface virtual heat flux, positive upward
wθv=SH./ρA./cpa.+humid_fac.*θv[1].*LH./ρA./Av;

#TF=(SH+LH)./ρA./cpa;
ustar=abs(ustar)
L=-θv[1].*ustar.^3 ./(karman.*gravity_mks.*wθv);
HD=ZD./L;
H=ZF[2:end]./L;

# function to calculate velocity VelocityScale

# surface velocity scale
out=VelocityScale.(0.1*ZD,ZD,L,ustar,θv[1],wθv)
ws=map(x->x[1], out)
wm=map(x->x[2], out)

θs=θv[1].+max.(b.*wθv./ws,0);

# interpolation to high resolution
spl = Spline1D(ZM, U); UD=spl(ZD);
spl = Spline1D(ZM, V);  VD=spl(ZD);
spl = Spline1D(ZM, θv); TD=spl(ZD);
#calculate PBLH at high resolution
ha=(Rcrit.*θs.*(UD.^2+VD.^2)./(gravity_mks.*(TD-θs)))-ZD;
ha[ha.>0].=Inf;
if all(isinf.(ha))
    if wθv>0
        id=SZM;
    else
        id=1;
    end
else
    ~,id=findmin(abs.(ha));
end
id=max(id,1);
h=abs(ZD[id[1]]);
# n - levels inside PBL
n=(h.-ZF).*((h.-ZF).>0);
n[n.==0].=Inf;
~,i=findmin(n);
# frac - fraction of the upper PBL level that is indide the PBL
frac=(h-ZF[i])./DZ[i];

# nsl - levels inside surface layer
nsl=(0.1 .*h.-ZF).*((0.1 .*h.-ZF).>0);
nsl[nsl.==0].=Inf;
~,isl=findmin(nsl);

#ipbl - levels inside the PBL
#npbl - levels outside the PBL
ipbl=[1:i;];
if i==SZM
    npbl=[];
else
    npbl=i+1:SZM;
end

kM=zeros(SZM); kS=zeros(SZM);
wm=zeros(SZM); wt=zeros(SZM);

out=VelocityScale.(ZM[ipbl],h,L,ustar,θv[1],wθv)

wt[ipbl]=map(x->x[1], out)
wm[ipbl]=map(x->x[2], out)
wstar=((gravity_mks/θv[1])*wθv*(wθv>0)*h)^(1/3)

kS[ipbl]=karman.*wt[ipbl].*min.(ZM[ipbl],h).*(1 .- min.(ZM[ipbl]/h,1)).^2
kM[ipbl]=karman.*wm[ipbl].*min.(ZM[ipbl],h).*(1 .- min.(ZM[ipbl]/h,1)).^2

kM[i]=kM[i].*frac
kS[i]=kS[i].*frac

#vertical diffusion above the PBL
S=max.(abs.(diff([0; U])./diff(ZF)),minshear);
GRi=(gravity_mks./θv[1:end-1]).*diff(θv)./diff(ZM)./S[1:end-1].^2;

fm=[0; holtslag_local.(GRi)]

lc=30

i2=min(i+1,SZM)
kM[npbl]=kM[npbl]+lc.^2 .*S[npbl].*fm[npbl];
kS[npbl]=kM[npbl]+lc.^2 .*S[npbl].*fm[npbl];

γcq=a*(LH./ρA./Av)/wm.^2*wstar/h
γct=a*(SH./ρA./cpa)/wm.^2*wstar/h
γcm=a*(ustar^2)/wstar*wm.^2/h

println(h)
return kM,kS,γcq,γct,γcm
#[i max(kM) max(kS) wθv max(Ts) U(1)]

end


text:
"""
    spawnrandompoint(ranfloat, randfloat) -> Array{Float64,1}

Initialize an [x,y] coordinate, uniformly random at water-enz interface, just inside enz.
"""
function spawnrandompoint(randfloat::Float64=rand(), randfloatstub::Float64=rand())
    spawnlinelength = SPAWN_RIGHT_X - SPAWN_LEFT_X + 2 * (SPAWN_ENZYME_MAX_Y - WALL_Y)
    spawncoordinate = spawnlinelength * randfloat - spawnlinelength / 2

    #Top Enzyme Face
    if SPAWN_LEFT_X <= spawncoordinate <= SPAWN_RIGHT_X
        return [spawncoordinate, SPAWN_ENZYME_MAX_Y]
    else
        randomy = SPAWN_ENZYME_MAX_Y + SPAWN_RIGHT_X - abs(spawncoordinate)
        if isapprox(randomy, WALL_Y)        # spawn outside wall
            randomy += BORDER_CORRECTION
        end

        if spawncoordinate < 0              #Left Enzyme Side
            return [SPAWN_LEFT_X, randomy]
        else                                #Top Enzyme Face
            return [SPAWN_RIGHT_X, randomy]
        end
    end
end

"""
    whereoutsidespawn(xy)

Assign a cardinal direction (W, NW, N, NE, or E) when crossing the water-enzyme boundary
"""
function whereoutsidespawn(xy::Array{Float64,1})
    x, y = xy
    if x < ENZYME_LEFT_X
        if y > ENZYME_MAX_Y
            return "NW"
        else
            return "W"
        end
    elseif x > ENZYME_RIGHT_X
        if y > ENZYME_MAX_Y
            return "NE"
        else
            return "E"
        end
    elseif y > ENZYME_MAX_Y
        return "N"
    else
        return "whereoutsidespawn Error"
    end
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
import IterativeSolvers

function lobpcg_itsolve(A, X0; prec=nothing, tol, maxiter, kwargs...)
    @warn "lobpcg_itsolve is not tested. Use at your own risk."
    res = IterativeSolvers.lobpcg(A, false, X0; P=prec, tol=tol, maxiter=maxiter)
    (λ=res.λ, X=res.X,
     residual_norms=res.residual_norms,
     iterations=res.iterations,
     converged=res.converged,
     n_matvec=0)
end


text:
function weight(PPDE, pThresh, alpha, precision)
	((PPDE >= 0) & (PPDE <= 1)) || error("PPDE has to be between 0 and 1")
	if PPDE==1
	        weight = -alpha * 2 * log(1-pThresh)
	else
	        weight = min(2 * (-log(1 - PPDE) + log(1 - pThresh)),
				-alpha * 2 * log(1 - pThresh))
	end

	return round(weight, sigdigits = precision)
end

max_stars_repo_name: htpusa/moominJL

max_issues_repo_name: htpusa/moominJL

max_forks_repo_name: htpusa/moominJL


text:
mktempdir() do tmpdir
    @info("Building 'nongit' in $tmpdir")
    cp(joinpath(@__DIR__, "docs"), joinpath(tmpdir, "docs"))
    include(joinpath(tmpdir, "docs/make.jl"))
    # Copy the build/ directory back so that it would be possible to inspect the output.
    cp(joinpath(tmpdir, "docs/build"), joinpath(@__DIR__, "build"); force=true)
end


text:
@testset "Steiner Tree" begin

    function sum_weight(g::AbstractGraph{<:Integer}, distmx::AbstractMatrix{<:Integer} = weights(g))
        sum_wt = zero(eltype(g))
        for e in edges(g)
            sum_wt += distmx[src(e), dst(e)]
        end
        return sum_wt
    end

    approx_factor(t::Integer) = 2-2/t

    g3 = star_graph(5)
    for g in testgraphs(g3)
        g_st = @inferred(steiner_tree(g, [2, 5]))
        @test ne(g_st) == 2 # [Edge(2, 1), Edge(1, 5)]

        g_copy = SimpleGraph(g)
        LightGraphs.filter_non_term_leaves!(g_copy, [2, 5])
        @test ne(g_copy) == 2 # [Edge(2, 1), Edge(1, 5)]
    end

    g4 = path_graph(11)
    for g in testgraphs(g4)
        g_st = @inferred(steiner_tree(g, [4, 8]))
        @test ne(g_st) == 4

        g_copy = SimpleGraph(g)
        LightGraphs.filter_non_term_leaves!(g_copy, [4, 8])
        @test ne(g_copy) == 4
    end

    g5 = grid([5, 5])
    for g in testgraphs(g5)
        g_st = @inferred(steiner_tree(g, [3, 11, 15, 23]))
        @test sum_weight(g_st) == 8
    end

    d = [0   2   3   4   5
         2   0   60   80  1
         3   60   0  120  150
         4   80  120  0  200
         5  1  150  200  0]

    g6 = complete_graph(5)

    for g in testgraphs(g6)
        g_st = @inferred(steiner_tree(g, [2, 4, 5], d))
        @test sum_weight(g_st, d) <= approx_factor(3)*(1+2+4)
    end

    d[2, 5] = d[5, 2] = 100

    for g in testgraphs(g6)
        g_st = @inferred(steiner_tree(g, [2, 4, 5], d))
        @test sum_weight(g_st, d) <= approx_factor(3)*(2+4+5)
    end

end


text:
using JuMP

function decision_variable(model::Model, S::States, d::Node, I_d::Vector{Node}, base_name::String="")
    # Create decision variables.
    dims = S[[I_d; d]]
    z_d = Array{VariableRef}(undef, dims...)
    for s in paths(dims)
        z_d[s...] = @variable(model, binary=true, base_name=base_name)
    end
    # Constraints to one decision per decision strategy.
    for s_I in paths(S[I_d])
        @constraint(model, sum(z_d[s_I..., s_d] for s_d in 1:S[d]) == 1)
    end
    return z_d
end

struct DecisionVariables
    D::Vector{Node}
    I_d::Vector{Vector{Node}}
    z::Vector{<:Array{VariableRef}}
end

"""
    DecisionVariables(model::Model,  diagram::InfluenceDiagram; names::Bool=false, name::String="z")

Create decision variables and constraints.

# Arguments
- `model::Model`: JuMP model into which variables are added.
- `diagram::InfluenceDiagram`: Influence diagram structure.
- `names::Bool`: Use names or have JuMP variables be anonymous.
- `name::String`: Prefix for predefined decision variable naming convention.

# Examples
```julia
z = DecisionVariables(model, diagram)
```
"""
function DecisionVariables(model::Model, diagram::InfluenceDiagram; names::Bool=false, name::String="z")
    DecisionVariables(diagram.D, diagram.I_j[diagram.D], [decision_variable(model, diagram.S, d, I_d, (names ? "$(name)_$(d.j)$(s)" : "")) for (d, I_d) in zip(diagram.D, diagram.I_j[diagram.D])])
end

function is_forbidden(s::Path, forbidden_paths::Vector{ForbiddenPath})
    return !all(s[k]∉v for (k, v) in forbidden_paths)
end

function path_compatibility_variable(model::Model, base_name::String="")
    # Create a path compatiblity variable
    x = @variable(model, base_name=base_name)

    # Constraint on the lower and upper bounds.
    @constraint(model, 0 ≤ x ≤ 1.0)

    return x
end

struct PathCompatibilityVariables{N} <: AbstractDict{Path{N}, VariableRef}
    data::Dict{Path{N}, VariableRef}
end

Base.length(x_s::PathCompatibilityVariables) = length(x_s.data)
Base.getindex(x_s::PathCompatibilityVariables, key) = getindex(x_s.data, key)
Base.get(x_s::PathCompatibilityVariables, key, default) = get(x_s.data, key, default)
Base.keys(x_s::PathCompatibilityVariables) = keys(x_s.data)
Base.values(x_s::PathCompatibilityVariables) = values(x_s.data)
Base.pairs(x_s::PathCompatibilityVariables) = pairs(x_s.data)
Base.iterate(x_s::PathCompatibilityVariables) = iterate(x_s.data)
Base.iterate(x_s::PathCompatibilityVariables, i) = iterate(x_s.data, i)

function decision_strategy_constraint(model::Model, S::States, d::Node, I_d::Vector{Node}, D::Vector{Node}, z::Array{VariableRef}, x_s::PathCompatibilityVariables)

    # states of nodes in information structure (s_d | s_I(d))
    dims = S[[I_d; d]]

    # Theoretical upper bound based on number of paths with information structure (s_d | s_I(d)) divided by number of possible decision strategies in other decision nodes
    other_decisions = filter(j -> all(j != d_set for d_set in [I_d; d]), D)
    theoretical_ub = prod(S)/prod(dims)/ prod(S[other_decisions])

    # paths that have a corresponding path compatibility variable
    existing_paths = keys(x_s)

    for s_d_s_Id in paths(dims) # iterate through all information states and states of d
        # paths with (s_d | s_I(d)) information structure
        feasible_paths = filter(s -> s[[I_d; d]] == s_d_s_Id, existing_paths)

        @constraint(model, sum(get(x_s, s, 0) for s in feasible_paths) ≤ z[s_d_s_Id...] * min(length(feasible_paths), theoretical_ub))
    end
end

"""
    PathCompatibilityVariables(model::Model,
        diagram::InfluenceDiagram,
        z::DecisionVariables;
        names::Bool=false,
        name::String="x",
        forbidden_paths::Vector{ForbiddenPath}=ForbiddenPath[],
        fixed::FixedPath=Dict{Node, State}(),
        probability_cut::Bool=true,
        probability_scale_factor::Float64=1.0)

Create path compatibility variables and constraints.

# Arguments
- `model::Model`: JuMP model into which variables are added.
- `diagram::InfluenceDiagram`: Influence diagram structure.
- `z::DecisionVariables`: Decision variables from `DecisionVariables` function.
- `names::Bool`: Use names or have JuMP variables be anonymous.
- `name::String`: Prefix for predefined decision variable naming convention.
- `forbidden_paths::Vector{ForbiddenPath}`: The forbidden subpath structures.
    Path compatibility variables will not be generated for paths that include
    forbidden subpaths.
- `fixed::FixedPath`: Path compatibility variable will not be generated
    for paths which do not include these fixed subpaths.
- `probability_cut` Includes probability cut constraint in the optimisation model.
- `probability_scale_factor::Float64`: Adjusts conditional value at risk model to
   be compatible with the expected value expression if the probabilities were scaled there.

# Examples
```julia
x_s = PathCompatibilityVariables(model, diagram; probability_cut = false)
```
"""
function PathCompatibilityVariables(model::Model,
    diagram::InfluenceDiagram,
    z::DecisionVariables;
    names::Bool=false,
    name::String="x",
    forbidden_paths::Vector{ForbiddenPath}=ForbiddenPath[],
    fixed::FixedPath=Dict{Node, State}(),
    probability_cut::Bool=true,
    probability_scale_factor::Float64=1.0)

    if probability_scale_factor ≤ 0
        throw(DomainError("The probability_scale_factor must be greater than 0."))
    end

    if !isempty(forbidden_paths)
        @warn("Forbidden paths is still an experimental feature.")
    end

    # Create path compatibility variable for each effective path.
    N = length(diagram.S)
    variables_x_s = Dict{Path{N}, VariableRef}(
        s => path_compatibility_variable(model, (names ? "$(name)$(s)" : ""))
        for s in paths(diagram.S, fixed)
        if !iszero(diagram.P(s)) && !is_forbidden(s, forbidden_paths)
    )

    x_s = PathCompatibilityVariables{N}(variables_x_s)

    # Add decision strategy constraints for each decision node
    for (d, z_d) in zip(z.D, z.z)
        decision_strategy_constraint(model, diagram.S, d, diagram.I_j[d], z.D, z_d, x_s)
    end

    if probability_cut
        @constraint(model, sum(x * diagram.P(s) * probability_scale_factor for (s, x) in x_s) == 1.0 * probability_scale_factor)
    end

    x_s
end

"""
    lazy_probability_cut(model::Model, diagram::InfluenceDiagram, x_s::PathCompatibilityVariables)

Add a probability cut to the model as a lazy constraint.

# Examples
```julia
lazy_probability_cut(model, diagram, x_s)
```

!!! note
    Remember to set lazy constraints on in the solver parameters, unless your solver does this automatically. Note that Gurobi does this automatically.

"""
function lazy_probability_cut(model::Model, diagram::InfluenceDiagram, x_s::PathCompatibilityVariables)
    # August 2021: The current implementation of JuMP doesn't allow multiple callback functions of the same type (e.g. lazy)
    # (see https://github.com/jump-dev/JuMP.jl/issues/2642)
    # What this means is that if you come up with a new lazy cut, you must replace this
    # function with a more general function (see discussion and solution in https://github.com/gamma-opt/DecisionProgramming.jl/issues/20)

    function probability_cut(cb_data)
        xsum = sum(callback_value(cb_data, x) * diagram.P(s) for (s, x) in x_s)
        if !isapprox(xsum, 1.0)
            con = @build_constraint(sum(x * diagram.P(s) for (s, x) in x_s) == 1.0)
            MOI.submit(model, MOI.LazyConstraint(cb_data), con)
        end
    end
    MOI.set(model, MOI.LazyConstraintCallback(), probability_cut)
end

"""
    expected_value(model::Model,
        diagram::InfluenceDiagram,
        x_s::PathCompatibilityVariables)

Create an expected value objective.

# Arguments
- `model::Model`: JuMP model into which variables are added.
- `diagram::InfluenceDiagram`: Influence diagram structure.
- `x_s::PathCompatibilityVariables`: Path compatibility variables.

# Examples
```julia
EV = expected_value(model, diagram, x_s)
```
"""
function expected_value(model::Model,
    diagram::InfluenceDiagram,
    x_s::PathCompatibilityVariables)

    @expression(model, sum(diagram.P(s) * x * diagram.U(s, diagram.translation) for (s, x) in x_s))
end

"""
    conditional_value_at_risk(model::Model,
        diagram,
        x_s::PathCompatibilityVariables{N},
        α::Float64;
        probability_scale_factor::Float64=1.0) where N

Create a conditional value-at-risk (CVaR) objective.

# Arguments
- `model::Model`: JuMP model into which variables are added.
- `diagram::InfluenceDiagram`: Influence diagram structure.
- `x_s::PathCompatibilityVariables`: Path compatibility variables.
- `α::Float64`: Probability level at which conditional value-at-risk is optimised.
- `probability_scale_factor::Float64`: Adjusts conditional value at risk model to
   be compatible with the expected value expression if the probabilities were scaled there.

# Examples
```julia
α = 0.05  # Parameter such that 0 ≤ α ≤ 1
CVaR = conditional_value_at_risk(model, x_s, U, P, α)
CVaR = conditional_value_at_risk(model, x_s, U, P, α; probability_scale_factor = 10.0)
```
"""
function conditional_value_at_risk(model::Model,
    diagram::InfluenceDiagram,
    x_s::PathCompatibilityVariables{N},
    α::Float64;
    probability_scale_factor::Float64=1.0) where N

    if probability_scale_factor ≤ 0
        throw(DomainError("The probability_scale_factor must be greater than 0."))
    end
    if !(0 < α ≤ 1)
        throw(DomainError("α should be 0 < α ≤ 1"))
    end

    # Pre-computed parameters
    u = collect(Iterators.flatten(diagram.U(s, diagram.translation) for s in keys(x_s)))
    u_sorted = sort(u)
    u_min = u_sorted[1]
    u_max = u_sorted[end]
    M = u_max - u_min
    u_diff = diff(u_sorted)
    if isempty(filter(!iszero, u_diff))
        return u_min    # All utilities are the same, CVaR is equal to that constant utility value
    else
        ϵ = minimum(filter(!iszero, abs.(u_diff))) / 2
    end

    # Variables and constraints
    η = @variable(model)
    @constraint(model, η ≥ u_min)
    @constraint(model, η ≤ u_max)
    ρ′_s = Dict{Path{N}, VariableRef}()
    for (s, x) in x_s
        u_s = diagram.U(s, diagram.translation)
        λ = @variable(model, binary=true)
        λ′ = @variable(model, binary=true)
        ρ = @variable(model)
        ρ′ = @variable(model)
        @constraint(model, η - u_s ≤ M * λ)
        @constraint(model, η - u_s ≥ (M + ϵ) * λ - M)
        @constraint(model, η - u_s ≤ (M + ϵ) * λ′ - ϵ)
        @constraint(model, η - u_s ≥ M * (λ′ - 1))
        @constraint(model, 0 ≤ ρ)
        @constraint(model, 0 ≤ ρ′)
        @constraint(model, ρ ≤ λ * probability_scale_factor)
        @constraint(model, ρ′ ≤ λ′* probability_scale_factor)
        @constraint(model, ρ ≤ ρ′)
        @constraint(model, ρ′ ≤ x * diagram.P(s) * probability_scale_factor)
        @constraint(model, (x * diagram.P(s) - (1 - λ))* probability_scale_factor ≤ ρ)
        ρ′_s[s] = ρ′
    end
    @constraint(model, sum(values(ρ′_s)) == α * probability_scale_factor)

    # Return CVaR as an expression
    CVaR = @expression(model, sum(ρ_bar * diagram.U(s, diagram.translation) for (s, ρ_bar) in ρ′_s) / (α * probability_scale_factor))

    return CVaR
end

# --- Construct decision strategy from JuMP variables ---

"""
    LocalDecisionStrategy(j::Node, z::Array{VariableRef})

Construct decision strategy from variable refs.
"""
function LocalDecisionStrategy(d::Node, z::Array{VariableRef})
    LocalDecisionStrategy(d, @. Int(round(value(z))))
end

"""
    DecisionStrategy(z::DecisionVariables)

Extract values for decision variables from solved decision model.

# Examples
```julia
Z = DecisionStrategy(z)
```
"""
function DecisionStrategy(z::DecisionVariables)
    DecisionStrategy(z.D, z.I_d, [LocalDecisionStrategy(d, z_var) for (d, z_var) in zip(z.D, z.z)])
end


text:
##########################
# FunctionCell Interface #
##########################
"""
	kmedoids(k, d=Distances.Euclidean() [;kwargs])

Constructs an untrained cell that when piped data inside, clusters the data using the
k-medoids algorithm where `k` is the desired number of clusters and `d` is the distance
that will be used to calculate distances from new samples to the clustering generated
from the untrained function cell.

# Keyword arguments (same as in `Clustering`)
  * `init` initialization algorithm; can be `:rand` for random, `:kmpp` for K-medoids++ or `:kmcen` for choosing samples with highest centrality (default `:kmpp`)
or a vector of size 'k' providing the indexes of the initial seeds.
  * `maxiter` maximum number of iterations (default `100`)
  * `tol` tolerable change in objective at convergence (default `1.0e-6`)
  * `display` Level of information; can be `:none`, `:final`, `:iter` (default `:none`)

Read the `Clustering.jl` documentation for more information.
"""
kmedoids(k::Int, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) = FunctionCell(kmedoids, (k,d), ModelProperties(), kwtitle("K-medoids ($k clusters)", kwargs); kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	kmedoids(x, k, d=Distances.Euclidean() [;kwargs])

Performs k-medoids clustering on the dataset `x` (square cost matrix) using `k` number of clusters and returns the clustering; when piping
data to the resulting clustering (e.g. trained cell), the distances from input data to the medoids are
returned, using `d` as distance.
"""
# Training
kmedoids(x::T where T<:CellData, k::Int, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) = kmedoids(getx!(x), k, d; kwargs...)
kmedoids(x::T where T<:AbstractVector, k::Int, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) =
	error("[kmedoids] Vectors are not supported as input data. Try using `v |> dist(v) |> kmedoids(n)`")
kmedoids(x::T where T<:AbstractMatrix, k::Int, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) = begin

	# Perform clustering
	kmeddata = Clustering.kmedoids(x, k; kwargs...)

	@assert length(kmeddata.medoids) == k "[kmedoids] The number of medoids (cluster centers) differes from the desired one k=$k."

	# Build model properties
	modelprops = ModelProperties(length(kmeddata.acosts), k)

	FunctionCell(kmedoids, Model((d, x[:,sort(kmeddata.medoids)], kmeddata), modelprops), kwtitle("K-medoids ($k clusters)", kwargs) );
end

# Execution
kmedoids(x::T where T<:CellData, model::Model{<:Tuple{<:Distances.PreMetric, <:AbstractArray, <:Clustering.KmedoidsResult}}) =
	datacell(kmedoids(getx!(x), model), gety(x))
kmedoids(x::T where T<:AbstractVector, model::Model{<:Tuple{<:Distances.PreMetric, <:AbstractArray, <:Clustering.KmedoidsResult}}) =
	kmedoids(mat(x, LearnBase.ObsDim.Constant{2}()), model)
kmedoids(x::T where T<:AbstractMatrix, model::Model{<:Tuple{<:Distances.PreMetric, <:AbstractArray, <:Clustering.KmedoidsResult}}) =
	# Return distances from each input sample to the medoids
	Distances.pairwise(model.data[1], model.data[2], x)

##########################
# FunctionCell Interface #
##########################
"""
	kmedoids!(medoids, d=Distances.Euclidean() [;kwargs])

Constructs an untrained cell that when piped data inside, clusters the data using the
k-medoids algorithm where `medoids` contains the initial seed observation indices
(will be updated inplace) and `d` is the distance that will be used to calculate distances
from new samples to the clustering generated from the untrained function cell.

# Keyword arguments (same as in `Clustering`)
  * `maxiter` maximum number of iterations (default `100`)
  * `tol` tolerable change in objective at convergence (default `1.0e-6`)
  * `display` Level of information; can be `:none`, `:final`, `:iter` (default `:none`)
"""
kmedoids!(medoids::Vector{Int}, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) =
	FunctionCell(kmedoids!, (medoids,d), ModelProperties(), kwtitle("K-medoids! ($(length(medoids)) clusters)", kwargs); kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	kmedoids!(x, medoids, d=Distances.Euclidean() [;kwargs])

Performs k-medoids clustering on the dataset `x` using `medoids` as initial seed observation indices
returning a clustering; when piping data to the resulting clustering (e.g. trained cell),
the distances from input data to the new medoids are returned, using `d` as distance.
"""
# Training
kmedoids!(x::T where T<:CellData, medoids::Vector{Int}, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) = kmedoids!(getx!(x), medoids, d; kwargs...)
kmedoids!(x::T where T<:AbstractVector, medoids::Vector{Int}, d::Distances.PreMetric=Distances.Euclidean() ; kwargs...) =
	error("[kmedoids!] Vectors are not supported as input data. Try using `v |> dist(v) |> kmedoids!(medoids)`")
kmedoids!(x::T where T<:AbstractMatrix, medoids::Vector{Int}, d::Distances.PreMetric=Distances.Euclidean() ; kwargs...) = begin

	# Perform clustering
	kmeddata = Clustering.kmedoids!(x, medoids; kwargs...)

	# Build model properties
	modelprops = ModelProperties(length(kmeddata.acosts), length(kmeddata.medoids))

	# Return a trained cell that uses `kmedoids` and not `kmedoids!` as execution function
	# since the only difference is in training (e.g. clustering) and not in
	# assigning new data to clusters.
	FunctionCell(kmedoids, Model((d, x[:,sort(kmeddata.medoids)], kmeddata), modelprops), kwtitle("K-medoids ($(length(medoids)) clusters)", kwargs) );
end


text:
 #-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*---------------------------------------------------------------------
#* Author: Gerson Bicca
#*$Id$
#*-------------------------------------------------------------------
type flow_measurement
	flow_measurement()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Type=>"PP"
			)),
			DanaInteger(),
			flow_mol(Dict{Symbol,Any}(
				:Brief=>"Flow Constant",
				:Hidden=>true
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0,
				:PosY=>0.78,
				:Protected=>true
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream",
				:PosX=>1,
				:PosY=>0.78,
				:Protected=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Flow Signal",
				:PosX=>0.5,
				:PosY=>0,
				:Protected=>true
			)),
			[
				:(Outlet.F = Inlet.F),
				:(Outlet.T = Inlet.T),
				:(Outlet.P = Inlet.P),
				:(Outlet.h = Inlet.h),
				:(Outlet.v = Inlet.v),
				:(Outlet.z = Inlet.z),
				:(FlowSignal*FlowConst = Inlet.F),
			],
			[
				"Overall Molar Balance","Temperature","Pressure","Energy Balance","Vapour fraction","Componente Molar Balance","Flow Measurement in kmol/h",
			],
			[:PP,:NComp,:FlowConst,],
			[:Inlet,:Outlet,:FlowSignal,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	FlowConst::flow_mol
	Inlet::stream
	Outlet::stream
	FlowSignal::control_signal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export flow_measurement
function set(in::flow_measurement)
	FlowConst = 1*"kmol/h"

end
function setEquationFlow(in::flow_measurement)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
end
function atributes(in::flow_measurement,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/flow_measurement"
	fields[:Brief]="Model to measure the flow rate in a pipe line. "
	drive!(fields,_)
	return fields
end
flow_measurement(_::Dict{Symbol,Any})=begin
	newModel=flow_measurement()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
### Methods for updating codebooks in ChainQ, LSQ, and LSQ++
### Julieta

# Update a dimension of a codebook using conjugate gradient either LSQR or LSMR
function updatecb!(
  K::SharedMatrix{Float32},
  C::SparseMatrixCSC{Int32,Int32},
  X::Matrix{Float32},
  IDX::UnitRange{Int64},
  codebook_upd_method::AbstractString="lsqr")   # choose the codebook update method out of lsqr or lsmr

  if codebook_upd_method == "lsqr"
    for i = IDX
      K[i,:] = IterativeSolvers.lsqr(C, vec(X[i, :]))
    end
  elseif codebook_upd_method == "lsmr"
    for i = IDX
      K[i,:] = IterativeSolvers.lsmr(C, vec(X[i, :]))
    end
  else
    error("Codebook update method unknown: ", codebook_upd_method)
  end
end

# Same as above but sparse matrix is parameterized differently
function updatecb!(
  K::SharedMatrix{Float32},
  C::SparseMatrixCSC{Float32,Int32},
  X::Matrix{Float32},
  IDX::UnitRange{Int64},
  codebook_upd_method::AbstractString="lsqr")   # choose the codebook update method out of lsqr or lsmr

  if codebook_upd_method == "lsqr"
    for i = IDX
      K[i,:] = IterativeSolvers.lsqr(C, vec(X[i, :]))
    end
  elseif codebook_upd_method == "lsmr"
    for i = IDX
      K[i,:] = IterativeSolvers.lsmr(C, vec(X[i, :]))
    end
  else
    error("Codebook update method unknown: ", codebook_upd_method)
  end
end

# Naive codebook update method. Super slow and nobody should use it
function update_codebooks_naive(
  X::Matrix{Float32}, # d-by-n matrix to update codebooks on.
  B::Matrix{Int16},   # m-by-n matrix. X encoded.
  h::Integer,         # number of entries per codebook.
  V::Bool=false)      # whether to print progress

  if V print("Doing Naive codebook update... "); st=time(); end

  m, _ = size( B )
  B_ok = sparsify_codes(B, h)
  C = convert(Matrix{Float32}, Matrix(B_ok) \ X')
  if V @printf("done in %.3f seconds.\n", time()-st); end
  return K2vec(collect(C'), m, h )
end

# Naive implementation of codebook update based on Cholesky decomposition
function update_codebooks_fast(
  X::Matrix{Float32}, # d-by-n matrix to update codebooks on.
  B::Matrix{Int16},   # m-by-n matrix. X encoded.
  h::Integer,         # number of entries per codebook.
  V::Bool=false,      # whether to print progress
  rho::AbstractFloat=1e-4) # regularization

  if V print("Doing fast codebook update... "); st=time(); end

  m, _ = size(B)

  # size( B_ok ) = n x mh
  B_ok = sparsify_codes(B, h)

  BTB = B_ok' * B_ok
  BTB = convert(Matrix{Float32}, Matrix( BTB ))

  BTXT = (X * B_ok)'

  # Solve sub-problem to solve C
  mh = size(BTB, 1)
  A = BTB + Matrix{Float32}(rho*I, mh, mh)
  b = BTXT

  # size( C ) = mh x d
  C = convert(Matrix{Float32}, A \ b)

  if V @printf("done in %.3f seconds.\n", time()-st); end

  return K2vec(collect(C'), m, h)
end

# Fast matrix multiplications that take advantage of the structure of B
function fast_bin_matmul(
  X::Matrix{T}, # d-by-n matrix to update codebooks on.
  B::Matrix{Int16},   # m-by-n matrix. X encoded.
  h::Integer,         # number of entries per codebook.
  V::Bool=false,      # whether to print progress
  rho::Float64=1e-4) where T <: AbstractFloat # regularization

  d, n = size(X)
  m, n = size(B)

  Bm = Vector{Vector{Int16}}(undef, m)
  for i = 1:m
    Bm[i] = B[i,:]
  end

  BTB = Matrix{Matrix{Float64}}(undef, m, m)

  # Diagonals are easy
  hi = zeros(Float32, h)
  @inbounds for i = 1:m

    # @simd for j = 1:h
    for j = 1:h
      hi[j] = zero(Float32)
    end

    for j = 1:n
      hi[ B[i,j] ] += one(Float32)
    end

    BTB[i,i] = Matrix(Diagonal(hi))
  end

  # Loop on off-diagonal quadrants
  @inbounds for i = 1:m
    Bmi = Bm[i]

    for j = (i+1):m

      Bmj = Bm[j]
      cij = zeros(Float32, h, h)

      for k=1:n
        cij[ Bmj[k], Bmi[k] ] += one(Float32)
      end

      # BTB is symmetric
      BTB[i,j] = cij
      BTB[j,i] = collect(cij')
    end
  end

  # Concatenate into a mh x mh matrix
  BTB = hvcat(m, BTB...)

  BXT = Vector{Matrix{Float64}}(undef, m)
  @inbounds for i=1:m
    BXTi = zeros(Float64, d, h)
    Bmi = Bm[i]
    for j=1:n
      Bmij = Bmi[j]
      @simd for k=1:d
        BXTi[ k, Bmij ] += X[k,j]
      end
    end
    BXT[i] = BXTi
  end

  BXT = hcat( BXT... )'

  # Solve with Cholesky
  A = BTB + rho*I
  b = collect(BXT)

  return A, b
end

# Same as above but matrix multiplications are accelerated by taking advantage
# of the binary structure of B
function update_codebooks_fast_bin(
  X::Matrix{Float32}, # d-by-n matrix to update codebooks on.
  B::Matrix{Int16},   # m-by-n matrix. X encoded.
  h::Integer,         # number of entries per codebook.
  V::Bool=false,      # whether to print progress
  rho::Float64=1e-4) # regularization

  if V print("Doing fast bin codebook update... "); st=time(); end

  m, n = size( B )
  d, n = size( X )

  A, b = fast_bin_matmul(X, B, h, V, rho)

  # Slightly more naive way to do it
  # C = A \ b

  # More low-level (avoids checks) and less readable but more efficient.
  # See https://software.intel.com/en-us/mkl-developer-reference-c-getrf
  lpt = LAPACK.getrf!(A)
  # See https://software.intel.com/en-us/mkl-developer-reference-c-getrs
  C = LAPACK.getrs!('N', lpt[1], lpt[2], b)

  # size( C ) = mh x d
  C = convert(Matrix{Float32}, C)

  if V @printf("done in %.3f seconds.\n", time()-st); end

  return K2vec(collect(C'), m, h)
end

# Same as above but matrix multiplications are accelerated by taking advantage
# of the binary structure of B
function update_codebooks_fast_bin2(
  X::Matrix{Float32}, # d-by-n matrix to update codebooks on.
  B::Matrix{Int16},   # m-by-n matrix. X encoded.
  h::Integer,         # number of entries per codebook.
  V::Bool=false,      # whether to print progress
  rho::Float64=1e-4) # regularization

  if V print("Doing fast bin codebook update... "); st=time(); end

  m, n = size( B )
  d, n = size( X )

  A, b = fast_bin_matmul(X, B, h, V, rho)

  # Slightly more naive way to do it
  # C = A \ b

  Ainv = inv(A)

  C = convert(Matrix{Float32}, Ainv * b)
  return K2vec(collect(C'), m, h)
end

###############################################
## Codebook update for a fully-connected MRF ##
###############################################

function update_codebooks(
  X::Matrix{Float32}, # d-by-n matrix to update codebooks on.
  B::Matrix{Int16},   # m-by-n matrix. X encoded.
  h::Integer,         # number of entries per codebook.
  V::Bool=false,      # whether to print progress
  method::AbstractString="fastbin")   # choose the codebook update method out of lsqr or lsmr

  if !(method in ["fast", "fastbin", "lsmr", "lsqr", "naive"]); error("Codebook update method unknown"); end

  if method == "fast"
    return update_codebooks_fast(X, B, h, V)
  elseif method == "fastbin"
    return update_codebooks_fast_bin(X, B, h, V)
  elseif method == "naive"
    return update_codebooks_naive(X, B, h, V)
  else # (method in ["lsqr", "lsmr"])

    if V print("Doing " * uppercase(method) * " codebook update... "); st=time(); end

    d, n   = size(X)
    m, _   = size(B)
    C = sparsify_codes( B, h )

    K = SharedArray{Float32}(d, size(C,2))
    if nworkers() == 1
      updatecb!( K, C, X, 1:d )
    else
      paridx = splitarray( 1:d, nworkers() )
      @sync begin
        for (i, wpid) in enumerate( workers() )
          @async begin
            remotecall_wait( updatecb!, wpid, K, C, X, paridx[i], method )
          end #@async
        end #for
      end #@sync
    end

    new_C = K2vec( K, m, h )

    if V @printf("done in %.3f seconds.\n", time()-st); end

    return new_C
  end
end # function update_codebooks

function get_cbdims_chain(
  d::Integer, # The number of dimensions
  m::Integer) # The number of codebooks

  subdims = splitarray(1:d, m-1)
  odims = Vector{UnitRange{Integer}}(undef, m)

  odims[1] =  subdims[1]
  for i = 2:m-1
    odims[i] = subdims[i-1][1]:subdims[i][end]
  end
  odims[end] = subdims[end]

  return odims
end

# Update a dimension of a codebook using LSQR
function updatecb_struct!(
  K::SharedMatrix{T},
  C::SparseMatrixCSC{T, Int32},
  X::Matrix{T},
  dim2C, #::Vector{Bool},
  subcbs,
  IDX::UnitRange{Int64}) where T <: AbstractFloat

  for i = IDX
    rcbs      = cat(subcbs[findall(dim2C[i,:])]..., dims=1)
    K[i,rcbs] = IterativeSolvers.lsqr(C[:,rcbs], vec(X[i, :]))
  end
end

function update_codebooks_generic(
  X::Matrix{T},  # d-by-n. The data that was encoded.
  B::Union{Matrix{Int16},SharedArray{Int16,2}}, # d-by-m. X encoded.
  h::Integer,          # number of entries per codebooks
  odimsfunc::Function, # Function that says which dimensions each codebook has
  V::Bool=false) where T <: AbstractFloat      # whether to print progress

  if V print("Doing LSQR codebook update... "); st=time(); end

  d, n   = size( X )
  m, _   = size( B )
  C = sparsify_codes(B, h, T)

  odims = odimsfunc(d, m)

  # Make a map of dimensions to codebooks
  dim2C = zeros(Bool, d, m)
  for i = 1:m; dim2C[odims[i], i] .= true; end

  K = SharedArray{T,2}((d, h*m))
  subcbs = splitarray(1:(h*m), m)

  if nworkers() == 1
    updatecb_struct!(K, C, X, dim2C, subcbs, 1:d)
  else
    paridx = splitarray(1:d, nworkers())
    @sync begin
      for (i, wpid) in enumerate(workers())
        @async begin
          remotecall_wait(updatecb_struct!, wpid, K, C, X, dim2C, subcbs, paridx[i])
        end #@async
      end #for
    end #@sync
  end

  new_C = K2vec( K, m, h )
  if V @printf("done in %.3f seconds.\n", time()-st) end

  return new_C
end

# Update codebooks in a chain
function update_codebooks_chain(
  X::Matrix{T}, # d-by-n. The data that was encoded.
  B::Union{Matrix{Int16},SharedArray{Int16,2}}, # d-by-m. X encoded.
  h::Integer,         # number of entries per codebook.
  V::Bool=false) where T <: AbstractFloat # whether to print progress

  st = time()
  C = update_codebooks_generic(X, B, h, get_cbdims_chain, V)

  return C, (time() - st)
end

# Fast version of the function above
function update_codebooks_chain_bin(
    X::Matrix{T}, # d-by-n. The data that was encoded.
    B::Matrix{Int16}, # d-by-m. X encoded.
    h::Integer,         # number of entries per codebook.
    V::Bool=false,      # whether to print progress
    rho::Float64=1e-4) where T <: AbstractFloat

    start_time = time_ns()

    m, n = size( B )
    d, n = size( X )

    A, b = fast_bin_matmul(X, B, h, V, rho)
    # @show size(A) # mh x mh
    # @show size(b) # mh x d

    # Take subdims that are adjacent and solve each one
    subdims = get_cbdims_chain(d, m)

    C = Vector{Matrix{Float32}}(undef, m)
    for i = 1:m
        C[i] = zeros(Float32, d, h)
    end

    for i = 1:(m-1)
        d1 = (i-1)*h+1:(i+1)*h
        d2 = intersect(subdims[i], subdims[i+1])

        # @show d1, d2

        # Slightly more naive way to do it
        # C = A \ b

        # More low-level (avoids checks) and less readable but more efficient.
        # See https://software.intel.com/en-us/mkl-developer-reference-c-getrf
        lpt = LAPACK.getrf!(A[d1, d1])
        # See https://software.intel.com/en-us/mkl-developer-reference-c-getrs
        Cblock = LAPACK.getrs!('N', lpt[1], lpt[2], b[d1, d2])

        C[i][d2, :] = Cblock[1:h, :]'
        C[i+1][d2, :] = Cblock[h+1:2*h, :]'
    end

    return C, (time_ns() - start_time) / 1e9
end