"""
Five algorithmic solutions to the maximum subarray problem (see Bentley's Programming Pearls).

Author: Michael Goldwasser

Module can be executed as a test harness.

Use -h flag for documentation on usage.
"""

#-----------------------------------------------------------------------
def algorithm1(A):
    """A brute-force algorithm using cubic-time.

    This algorithm tries every possible i < j pair, and computes the
    sum of entries in A[i:j].

    Returns tuple designating optimal (maxsum, start, stop).
    """
    n = len(A)
    maxsofar = (0, 0, 0)                    # (t, start, stop) with sum(A[start:stop]) = t
    for i in range(n):                      # potential start index
        for j in range(i+1, n+1):           # potential stop index
            total = 0
            for k in range(i,j):            # compute sum(A[i:j])
                total += A[k]
            if total > maxsofar[0]:         # if new best, record parameters
                maxsofar = (total, i, j)
    return maxsofar

#-----------------------------------------------------------------------
def algorithm2a(A):
    """A quadratic time approach to brute force by using previous sums.

    Uses fact that sum(A[i:j]) = sum(A[i:j-1]) + A[j-1] to avoid a third nested loop.

    Returns tuple designating optimal (maxsum, start, stop).
    """
    n = len(A)
    maxsofar = (0, 0, 0)                    # (t, start, stop) with sum(A[start:stop]) = t
    for i in range(n):                      # potential start index
        total = 0                           # maintain sum(A[i:j]) as j increases
        for j in range(i+1, n+1):
            total += A[j-1]                 # thus total is now sum(A[i:j])
            if total > maxsofar[0]:
                maxsofar = (total, i, j)
    return maxsofar

#-----------------------------------------------------------------------
def algorithm2b(A):
    """A quadratic time approach to brute force using prefix sums.

    precomputes sum(A[0:j]) for each j.

    That allows for O(1) time computation of sum(A[i:j]) = sum(A[0:j]) - sum(A[0:i])

    Returns tuple designating optimal (maxsum, start, stop).
    """
    n = len(A)

    cumulative = [0] * (n+1)                 # cumulative[i] = sum(A[0:i])
    for i in range(n):
        cumulative[i+1] = cumulative[i] + A[i]
    
    maxsofar = (0, 0, 0)                    # (t, start, stop) with sum(A[start:stop]) = t
    for i in range(n):                      # potential start index
        for j in range(i+1, n+1):
            total = cumulative[j] - cumulative[i]
            if total > maxsofar[0]:
                maxsofar = (total, i, j)
    return maxsofar

#-----------------------------------------------------------------------
def recurse(A, start, stop):
    """Recursion for algorithm3 that only considers implicit slice A[start:stop])."""
    if stop == start:                         # zero elements
        return (0, 0, 0)               
    elif stop == start + 1:                   # one element
        if A[start] > 0:
            return (A[start], start, stop)
        else:
            return (0, start, start)
    else:                                     # two or more elements
        mid = (start + stop) // 2

        # find maximum sum(A[i:mid]) for i < mid
        total = 0
        lmax = (0, mid)                       # (t,i) such that sum(A[i:mid]) = t
        for i in range(mid-1, start-1, -1):
            total += A[i]
            if total > lmax[0]:
                lmax = (total,i)
                
        # find maximum sum(A[mid:j]) for j > mid
        total = 0
        rmax = (0, mid)                       # (t, j) such that sum(A[mid:j]) = t
        for j in range(mid+1, stop+1):
            total += A[j-1]
            if total > rmax[0]:
                rmax = (total,j)

        overlay = (lmax[0]+rmax[0], lmax[1], rmax[1])

        return max(recurse(A, start, mid),
                   recurse(A, mid, stop),
                   overlay)

def algorithm3(A):
    """A divide-and-conquer approach achieving O(n log n) time.

    We find the maximum solution from the left half, the maximum from the right,
    and the maximum solution that straddles the middle.  One of those three is
    the true optimal solution.

    Returns tuple designating optimal (maxsum, start, stop).
    """
    return recurse(A, 0, len(A))

#-----------------------------------------------------------------------
def algorithm4(A):
    """A linear algoirthm based on a simple greedy strategy.

    As we let index i increase, we keep track of two pieces of information:
      * maxhere     which is maximum subarray ending specifically at index i
      * maxsofar    which is maximum subarray for anything from A[0:i]

    Returns tuple designating optimal (maxsum, start, stop).
    """
    n = len(A)
    maxsofar = (0, 0, 0)           # (t, start, stop) such that sum(A[start,stop]) = t
    maxhere = (0, 0)               # (t, start) such that sum(A[start:i]) = t
    for i in range(1, n+1):
        if maxhere[0] + A[i-1] > 0:
            maxhere = (maxhere[0] + A[i-1], maxhere[1])
        else:
            maxhere = (0, i)
        if maxhere[0] > maxsofar[0]:
            maxsofar = (maxhere[0], maxhere[1], i)
    return maxsofar


#-----------------------------------------------------------------------
if __name__ == '__main__':
    from optparse import OptionParser
    import sys
    import time
    import random

    def quit(message=None):
        if message: print(message)
        sys.exit(1)

    parser = OptionParser(usage='usage: %prog [options]')

    parser.add_option('-n', dest='size', type='int', default=512,
                      help='number of elements in data set [default: %default]')
    
    parser.add_option('-a', dest='num_alg', metavar='NUM', type='int', default=5,
                      help='number of algorithms to test [default: %default]')
    
    parser.add_option('-s', dest='seed', type='int', default=None,
                      help='random seed for data set [default: %default]')   

    (options,args) = parser.parse_args()

    if options.size <= 0:
        quit('n must be positive')
    
    if not 1 <= options.num_alg <= 5:
        quit('invalid number of algorithms to test')
    
    if options.seed is None:
        options.seed = random.randrange(1000000)
#        print('Using seed: {0}'.format(options.seed))
    random.seed(options.seed)

    data = [random.uniform(-100,100) for _ in range(options.size)]
          
    algorithms = (algorithm1, algorithm2a, algorithm2b, algorithm3, algorithm4)

    print('Running tests for {0} algorithms using n={1} and seed={2}'.format(options.num_alg,options.size,options.seed))
    for alg in reversed(algorithms[len(algorithms)-options.num_alg:]):
        start = time.time()
        answer = alg(data)
        end = time.time()
        print("{0:<11} found sum(data[{1}:{2}])={3:.2f} using {4:.3f} seconds of computation".format(
            alg.__name__, answer[1], answer[2], answer[0], (end-start)))