AggarwalParkMatrixSearch.cs

Go to the documentation of this file.

namespace ULSAlgorithms.Exact.AggarwalPark.Internal;
 
/// <summary>
/// SMAWK row-minimum search for the implicit Monge matrices generated by the
/// Aggarwal-Park divide-and-conquer lot-sizing recursion.
/// </summary>
/// <remarks>
/// <para>
/// Rows are cumulative-demand query points in nondecreasing order and columns
/// are predecessor lines in nonincreasing transformed marginal-cost order.
/// The implicit entry is
/// <c>intercept[column] + slope[column] * prefixDemand[row]</c>.
/// This matrix is Monge, so all row minima are obtained in linear time in the
/// number of rows plus columns.
/// </para>
/// <para>
/// The implementation stores only reduced column indices. Row subsets created
/// by the SMAWK recursion are represented arithmetically by
/// <c>(rowStart, rowStep, rowCount)</c>, avoiding temporary row arrays.
/// </para>
/// </remarks>
internal static class AggarwalParkMatrixSearch
{
    public static void FindRowMinima(
        int rowStart,
        int rowCount,
        ReadOnlySpan<int> columns,
        ReadOnlySpan<double> prefixDemand,
        ReadOnlySpan<double> slope,
        ReadOnlySpan<double> intercept,
        Span<int> argMin,
        Span<int> workspace)
    {
        if (rowCount <= 0)
        {
            return;
        }
 
        if (columns.IsEmpty)
        {
            throw new ArgumentException(
                "At least one predecessor column is required.",
                nameof(columns));
        }
 
        var requiredWorkspace = checked((2 * rowCount) + 2);
 
        if (workspace.Length < requiredWorkspace)
        {
            throw new ArgumentException(
                $"SMAWK workspace must contain at least {requiredWorkspace} entries.",
                nameof(workspace));
        }
 
        Search(
            rowStart,
            rowStep: 1,
            rowCount,
            columns,
            prefixDemand,
            slope,
            intercept,
            argMin,
            workspace);
    }
 
    private static void Search(
        int rowStart,
        int rowStep,
        int rowCount,
        ReadOnlySpan<int> columns,
        ReadOnlySpan<double> prefixDemand,
        ReadOnlySpan<double> slope,
        ReadOnlySpan<double> intercept,
        Span<int> argMin,
        Span<int> workspace)
    {
        if (rowCount == 0)
        {
            return;
        }
 
        var reducedCount = 0;
 
        for (var columnPosition = 0;
             columnPosition < columns.Length;
             columnPosition++)
        {
            var column = columns[columnPosition];
 
            while (reducedCount > 0)
            {
                var comparisonRowPosition = reducedCount - 1;
 
                if (comparisonRowPosition >= rowCount)
                {
                    break;
                }
 
                var row =
                    rowStart +
                    (comparisonRowPosition * rowStep);
 
                var oldColumn =
                    workspace[reducedCount - 1];
 
                if (!IsStrictlyBetter(
                        row,
                        column,
                        oldColumn,
                        prefixDemand,
                        slope,
                        intercept))
                {
                    break;
                }
 
                reducedCount--;
            }
 
            if (reducedCount < rowCount)
            {
                workspace[reducedCount] = column;
                reducedCount++;
            }
        }
 
        var reducedColumns =
            workspace[..reducedCount];
 
        var oddRowCount = rowCount / 2;
 
        if (oddRowCount > 0)
        {
            Search(
                rowStart + rowStep,
                checked(rowStep * 2),
                oddRowCount,
                reducedColumns,
                prefixDemand,
                slope,
                intercept,
                argMin,
                workspace[reducedCount..]);
        }
 
        var lowerPosition = 0;
 
        for (var rowPosition = 0;
             rowPosition < rowCount;
             rowPosition += 2)
        {
            var row =
                rowStart +
                (rowPosition * rowStep);
 
            if (rowPosition > 0)
            {
                var lowerColumn =
                    argMin[row - rowStep];
 
                while (lowerPosition < reducedCount &&
                       reducedColumns[lowerPosition] != lowerColumn)
                {
                    lowerPosition++;
                }
 
                if (lowerPosition >= reducedCount)
                {
                    throw new InvalidOperationException(
                        "SMAWK lower-bound column was not found.");
                }
            }
 
            var upperPosition = reducedCount - 1;
 
            if (rowPosition + 1 < rowCount)
            {
                var upperRow = row + rowStep;
                var upperColumn = argMin[upperRow];
 
                upperPosition = lowerPosition;
 
                while (upperPosition < reducedCount &&
                       reducedColumns[upperPosition] != upperColumn)
                {
                    upperPosition++;
                }
 
                if (upperPosition >= reducedCount)
                {
                    throw new InvalidOperationException(
                        "SMAWK upper-bound column was not found.");
                }
            }
 
            var bestPosition = lowerPosition;
 
            for (var position = lowerPosition + 1;
                 position <= upperPosition;
                 position++)
            {
                if (IsStrictlyBetter(
                        row,
                        reducedColumns[position],
                        reducedColumns[bestPosition],
                        prefixDemand,
                        slope,
                        intercept))
                {
                    bestPosition = position;
                }
            }
 
            argMin[row] = reducedColumns[bestPosition];
 
            if (rowPosition + 1 < rowCount)
            {
                lowerPosition = upperPosition;
            }
        }
    }
 
    private static bool IsStrictlyBetter(
        int row,
        int candidateColumn,
        int currentColumn,
        ReadOnlySpan<double> prefixDemand,
        ReadOnlySpan<double> slope,
        ReadOnlySpan<double> intercept)
    {
        var candidate = Evaluate(
            row,
            candidateColumn,
            prefixDemand,
            slope,
            intercept);
 
        var current = Evaluate(
            row,
            currentColumn,
            prefixDemand,
            slope,
            intercept);
 
        // Strict comparison keeps the left-most column when values tie, which
        // is the tie convention required by monotone matrix searching.
        return candidate < current;
    }
 
    public static double Evaluate(
        int row,
        int column,
        ReadOnlySpan<double> prefixDemand,
        ReadOnlySpan<double> slope,
        ReadOnlySpan<double> intercept)
    {
        var product =
            slope[column] *
            prefixDemand[row];
 
        var value =
            intercept[column] +
            product;
 
        if (!double.IsFinite(product) ||
            !double.IsFinite(value))
        {
            throw new ArithmeticException(
                "Numerical overflow while evaluating an Aggarwal-Park Monge-matrix entry.");
        }
 
        return value;
    }
}