<jupyter_start><jupyter_text>The ProblemSuppose you have a large number of items, each of which has a size, in a particular order. For an example, consider alist of 3114 U.S. Counties sorted by median income, with county population as the size.-Fig bar chart-Now, suppose you want to divide that list into a fixed number of partitions, as evenly as possible. Here, "evenly" means"into a set of bins with the lowest standard deviation in bin size".-Fig bar chart with colored bands left to right-This problem is the subject of this series of articles. It has the advantage of being a very simple problem,at its root: The only way to find the solution is to try every valid combination.That property, in theory, can make finding a solution impractically slow for large numbers of items.However, through the use of dynamic programming I built an initial Python implementation that solved the problem ofputting all U.S. counties into 15 buckets in 84 seconds. Through simple use of Numba, an LLVM-based JIT compiler, Ibrought this down to less than a second. By using Numba's CUDA extensions, I was able to solve it in 28 millisecondson a GeForce 1080Ti.The simplicity of the problem makes for a good case study. I will use it to explain the basics of using dynamicprogramming to speed up whole classes of problems. I will also show how to use Numba to greatly accelerate manycompute-intensive algorithms, and how to use NVidia's CUDA technology to go even faster.The algorithm is simple in another sense, as well: because it checks every possible combination, its performance inspace and time does not depend on the actual data, per se, only the number of items and the number of partitions.  Give Me Solutions, Not ProblemsYou can find my library, Histoptimizer, at link. It provides multiple interfaces to the solver for use with Numpy orPandas, as well as a CLI that can be used to decorate CSVs with partitioning information. Motivating Use CaseAs part of my Covid hobby I wanted to show a 3-dimensional surface graph, with date on the X dimension, and countiesdistributed along the Y dimension in order of population density.The Z measurement would be the trailing 7-day-average of new cases.Covid is stochastic. I felt by aggregating counties into bins I might be able to see if there were population density-based effects on Covid infectivity.I wanted each span on the Y axis to represent approximately the same number of people, so that the audience for thegraph could reason volumetrically about the impact on the U.S. Population.I wanted to partition the ordered series of counties into bins such that the number of people represented by each of theribbons in the graph below were as equal as possible.You can visualize this in a couple of ways, like these bar charts of an ideal sum of items per partition or a cumulativesum. We'll be working with the former.This turns out to be a pretty hard goal to achieve!U.S. County population is very unequally distributed.I tried my standard practice of picking partition boundaries. In a lot of cases, where the absolute values of the sortmetric are important--e.g., some measurements are in tolerance, and some are out--then the boundary values of the binsare pre-determined.Log-based bin boundaries improve things, but not enough.In this case, however, I don't really care whether a bin boundary is at 962 people per km ... Minimizing Variance (and Standard Deviation)My first thought was that I wanted a distribution that minimized the standard deviation of the population of each of thebins. And this is correct, but standard deviation is the square root of the variance, so by minimizing variance we arealso minimizing standard deviation. We can therefore save a needless square root instruction when evaluating cost. When Would You Use This SolutionAny time you have an ordered population of discrete items with an associated set of "sizes", and you would like topartition them into a fixed number of bins or buckets such that the boundaries of the bins are not fixed, and the totalsize of items in each bucket is as equal as possible, as measured by variance or standard deviation.By looping over a range of bin counts, it is also possible to identify a number of bins over a finite range that yields the lowest total variance. What is this solution based onThe basic algorithm is lifted from Skiena's Algorithm Design Handbook chapter on Dynamic Programming. A variance costfunction is substituted for the max sum cost function used in that work.Two accelerated parallel implementations of the algorithm are provided, one CPU-based and one GPU-based:One depends on the Numba JIT compilation package to compile Python to statically-typed native code.The second also depends on Numba, but relies on the CUDA development kit and compatible NVidia hardware to achievegreater acceleration. What the solution doesConsists of three APIs that provide convenient ways to access the functionality using Numpy or Pandas.Partitioners : Take a Numpy array and a number of buckets, returns an array of divider locations.Series: Given a Pandas Series and a number of buckets, return a  series that maps each row in the input series to a bucket numberDataFrame: Given a Pandas Dataframe, add one or more columns that partition the dataframe into various numbers of buckets.It also provides an CLI that can be used to add columns to CSV data, as well as a simple benchmarking tool. Recursive Solutions Traversing all possible combinations Building a Naive SolutionTo be able to test for correctness. Recurrence RelationBase Case. A simple recursive solution. Dynamic Programming Adding a simple cache to our recursive solution That's all Dynamic Programming IsCaching solutions to sub-problems when they will be needed over and over again. Eliminating Recursion -- Skiena's AlgorithmDifference from the maximum sum version.This is our first practical solution. Improving Performance at the EdgesEliminating some unnecessary comparisons due to invalid edges. Accelerating Numeric Computation with NumbaWe're going to take our solution from XX seconds down to 1 using a sophisticated tool that is simple to use. Intro to NumbaWhat it do. Adapting the algorithm to Numba ExperimentsThey weren't very successful. Using guvectorize Moving more loop into the kernel. CUDAWhat it do. Links to Numba tutorial and programming guide. A first attempt -- parallelizing buckets Improving Occupancy -- parallelizing items. Better utilization -- sharding comparisonsBe sure to talk about memory locality, and also perform a test that doesn't have it.Limits of improvement. Minimum Variance Linear Partition ProblemProblem: Integer Partition to minimize σ2 without RearrangementInput: An arrangement S of nonnegative numbers {s1, . . . , sn} and an integer k.Output: Partition S into k ranges to minimize the variance (σ2), without reordering any of the numbers.Old:Given an ordered series of items of various sizes, insert *n* dividers between items to create n+1 buckets insuch a way that the variance (σ2) of the sum of the sizes of the items in each bucket is as low as possible. Let: * *k* be a number of partitions * j be the number of dividers necessary to create *k* partitions, k - 1 * *n* be a number of items * *s1..n* be an ordered series of item sizes. Note that this is 1-based. * *c1..k* be a vector giving k divider locations that come *after* the corresponding (1-based) item locations in *z*,given in increasing order s.t. ci-1 i i+1 for 2 k-1, and ck-1 k * *c0 = 0* and *ck+1 = n*. These locations represent the implicit dividers placed at the beginning and endof the series of items. * *mean* be the sum of the sizes in *s* divided by k * *p1..n* be the prefix sums of items in s,   such that *p0* = 0 and *pi* = *pi-1* + *si* for 1 < i < n. * *size(a,b)* be a function giving the sum of the item sizes in a bucket containing items from *s* on a half-closed interval(containing items *sa* to *sb-1*). This can be quickly given by *pb* - *pa-1* * *variance(a,b)* be the square of the difference between *size(*a,b)* and the value *mean*. This measuresfor a bucket on the half-closed interval the square of the difference of the bucket size and the *overall* mean bucket size.It can be computed from s: *variance(a,b)* = (*size(*a,b)* - *mean*)2* * *buckets(c)* be a function giving a vector *b1..k* of sums of items in each bucket for a given divider vector *c*:   b1 = size(*1,c1), b2 = *size(*c1,c2)*, . . . ,   bk = size(*cj,n) ProblemFind a vector c that minimizes min σ2(b(c)) for given s1s2...snFind a vector of divider locations c that minimizes the statistical variance of the sum of the items in each bucket. Algorithm Hypothesis:To find the vector of divider locations that minimizes the variance, if we can find the optimal locationof the last divider, we will learn the optimal location of the previous divider, and this is true forevery divider, such that we can find the vector that optimizes this value recursively. Proof:1. For 1 cost(c,x)* be the sum of *v* for items the first n divider locations in c,   cost(c,x) = variance(1,c1) + variance(c1,c2) + ... + variance(cx-1, cx)2. *cost(c,j)* is equal to the variance of the size of the buckets defined by *c*.   By finding *c* that minimizes *cost(c,j)* we can solve the problem.3. If:   cost(c,x)    = variance(1,c1) + variance(c1,c2) + ... + variance(cx-2, cx-1) + variance(cx-1, cx)   Then:   cost(c,x-1) = variance(1,c1) + variance(c1,c2) + ... + variance(cx-2, cx-1)   All terms except the last are shared.4. Substitute for the common terms to obtain a recurrence relation:   For *x* = 1:   cost(c, 1) = variance(1,c1)   and for 1 < *x* <= *k*   cost(c,x) = cost(c,x-1) + variance(cx-1 + 1, cx)   If we wish to determine the variance for a vector c (*cost(c,j)*), this becomes:   cost(c,j) = cost(c,j-1) + variance(cj-1 + 1, cj)   The increase in variance created by the last divider is not just dependent on the location of the last divider.   The variance is a function of the size of the *buckets*, and the last divider sits between two buckets. For a given   *cj*, the size of the last bucket is fixed, as it contains items from *z[cj]* to *sn*. The next-to-last bucket,   however, contains items from *z[cj-1]* to *z[cj-1]*, and there are usually many valid values of *cj-1* for any   given *cj*. In order to minimize the added variance we need to find a *combination* of *cj-1* and *cj* that   minimizes the overall variance.   A collection of *n* items has *(n-1)* spaces between items in which to place a divider, and at least the first *(j-2)*   locations must contain previous dividers. Therefore there are only *(n-1)-(j-2)* locations we can place the last two   dividers *cj-1* and *cj*. Therefore there are there are *comb((n-1)-(j-1), 2)* combinations   of divider locations. This is a very manageable number for many uses cases on modern hardware, and by trying all   possible values of cj - 1 for each feasible value of cj we can find the optimal value for both.   Note that in order to find the value of cj that minimizes the variance, we have to recursively determine   the value of element cj-1 that minimizes *cost(c, j - 1)*, *for each possible value of *cj**. But this is   a smaller version of the same problem, and does not depend on the value of cj!   Minimizing this function requires computing the minimum value of cost(c,j-1), but note that this is a smaller version of the   same problem, and the function does not depend on the value of cj! Because the value of *cost(c,j-1)* will be   evaluated multiple times, once for each value of *cj*, we can save considerable time by caching the result.   This means that we can minimize variance(c,j) by minimizing variance(c,j-1), variance(c,j-2), ..., variance(c,1) in turn.   We can derive the vector c by recording the values of cn and cn-1 that minimize the variance at each step when   we determine them.   Notes:   The lowest variance partition of S is a list of divider locations, and that list has a last member. We will find the   optimal location for that last member, and in the process discover the optimal locations of each of the previous members.Will finding the value that minimies variance(c, n-1)Example:Consider s = [7, 41, 97, 53, 67, 24]and j = 2Then p = [0, 7, 41, 138, 191, 258, 282, 0]<jupyter_code>import numpy as np

#items = [7, 41, 97, 53, 67, 24]
items = [1, 3, 5, 7, 9]
num_items = len(items)
num_dividers = 2
num_buckets = num_dividers + 1
prefix_sum = [0] * (num_items + 1)
# Add a 0 at the beginning to make item offsets 1-based.
items = [0] + items

# insert a 0 at the beginning of the partition list so that the first divider
# location is at c[1]. Append a "divider" at the last location as well.
c = [0] + [1, 3] + [num_items]

for item in range(1, num_items + 1):
    prefix_sum[item] = prefix_sum[item - 1] + items[item]

# Last value of the prefix sum is everything added together.
mean = prefix_sum[-1] / num_buckets
def s(a,b):
    return prefix_sum[b] - prefix_sum[a-1]

# Variance is the square of the difference between the mean of the sum
def v(a,b):
    (s(a,b) - mean)**2

def bucket_sums(c, n):
    """
    Assumes c has a 0 prepended, divider locations one-based.
    """
    sums = [0]*(n+1)
    for i in range(0, n+1):
        sums[i] = s(c[i] + 1, c[i+1])
    return sums

def variance(c, n):
    if n == 1:
        return v(0, c[1])
    else:
        return variance(c, n-1) + v(c[n-1], c[n])
bucket_sums(c, 1)
import matplotlib.pyplot as plt
import numpy as np


# prepare some coordinates
x, y, z = np.indices((8, 8, 8))

# draw cuboids in the top left and bottom right corners, and a link between them
cube1 = (x < 3) & (y < 3) & (z < 3)
cube2 = (x >= 5) & (y >= 5) & (z >= 5)
link = abs(x - y) + abs(y - z) + abs(z - x) <= 2

# combine the objects into a single boolean array
voxels = cube1 | cube2 | link

# set the colors of each object
colors = np.empty(voxels.shape, dtype=object)
colors[link] = 'red'
colors[cube1] = 'blue'
colors[cube2] = 'green'

# and plot everything
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.voxels(voxels, facecolors=colors, edgecolor='k')

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.06981997882103777, 'fraction_numerical': 0.014542887398517473, 'mean_word_length': 4.136330674401741, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '538445', 'license': '0BSD', 'n_tokens_mistral': 4075, 'n_tokens_neox': 3824, 'n_words': 2218, 'path': 'MSDiPoOS.ipynb', 'repo_name': 'delusionary/histoptimizer', 'size': '14165'}