Suppose you need n_1 in [10,30], n_2 in [20,40], n_3 in [30,50] and n1+n2+n3=90
If you need each possible triplet (n_1, n_2, n_3) to be equally-likely, that's going to be difficult. The number of triples of the form (20, n_2, n_3) is greater than the number of triples (10, n_2, n_3), so you can't just pick n_1 uniformly.
The incredibly slow but accurate way is to generate the all 5 randoms in the correct ranges and reject the whole group if the sum is not correct.
. . . AHA!
I found a way to parametrize the choice effectively. First, though, for simplicity note that the sum of the low bounds is the minimum possible sum. If subtract the sum of the low bounds from the target number and subtract the low bound from each generated number, you get a problem where each number is in the interval [0, max_k-min_k]. That simplifies the math and array (list) handling. Let n_k be the 0-based choice with 0<=n_k<=max_k-min_k.
The order of the sums is lexicographic, with all sums beginning with n_1=0 (if any) first, then n_1==1 sums, etc. Sums are sorted by n_2 in each of those groups, then by n_3, and so on. If you know how many sums add to the target (call that T), and how many sums start with n_1=0, 1, 2, ... then you can find the starting number n1 of sum number S in in that list. Then you can reduce the problem to adding n_2+n_3+... to get T-n_1, finding sum number S - (number original sums starting with number less than n_1).
Let pulse(n) be a list of n+1 ones: (n+1)*[1] in Python terms. Let max_k,min_k be the limits for the k'th choice, and m_k = max_k-min_k be the upper limit for 0-based choices. Then there are 1+m_1 different "sums" from the choice of the first number, and pulse(m_k) gives the distribution: 1 was to make each sum from 0 to m_1. For the first two choices, there are m_1+m_+1 different sums. It turns out that the convolution of pulse(m_1) with pulse(m_2) gives the distribution.
Time to stop for some code:
def pulse(width, value=1):
''' Returns a vector of (width+1) integer ones. '''
return (width+1)*[value]
def stepconv(vector, width):
''' Computes the discrete convolution of vector with a "unit"
pulse of given width.
Formula: result[i] = Sum[j=0 to width] 1*vector[i-j]
Where 0 <= i <= len(vector)+width-1, and the "1*" is the value
of the implied unit pulse function: pulse[j] = 1 for 0<=j<=width.
'''
result = width*[0] + vector;
for i in range(len(vector)):
result[i] = sum(result[i:i+width+1])
for i in range(len(vector), len(result)):
result[i] = sum(result[i:])
return result
That's coded specifically for only doing convolutions with a "pulse" array, so every linear combination in the convolution is just a sum.
Those are used only in the constructor of the final class solution:
class ConstrainedRandom(object):
def __init__(self, ranges=None, target=None, seed=None):
self._rand = random.Random(seed)
if ranges != None: self.setrange(ranges)
if target != None: self.settarget(target)
def setrange(self, ranges):
self._ranges = ranges
self._nranges = len(self._ranges)
self._nmin, self._nmax = zip(*self._ranges)
self._minsum = sum(self._nmin)
self._maxsum = sum(self._nmax)
self._zmax = [y-x for x,y in self._ranges]
self._rconv = self._nranges * [None]
self._rconv[-1] = pulse(self._zmax[-1])
for k in range(self._nranges-1, 0, -1):
self._rconv[k-1] = stepconv(self._rconv[k], self._zmax[k-1])
def settarget(self, target):
self._target = target
def next(self, target=None):
k = target if target != None else self._target
k = k - self._minsum;
N = self._rconv[0][k]
seq = self._rand.randint(0,N-1)
result = self._nranges*[0]
for i in range(len(result)-1):
cv = self._rconv[i+1]
r_i = 0
while k >= len(cv):
r_i += 1
k -= 1
while cv[k] <= seq:
seq -= cv[k]
r_i += 1
k -= 1
result[i] = r_i
result[-1] = k # t
return [x+y for x,y in zip(result, self._nmin)]
# end clss ConstrainedRandom
Use that with:
ranges = [(low, high), (low, high), ...]
cr = ConstrainedRandom(ranges, target)
seq = cr.next();
print(seq)
assert sum(seq)==target
seq = cr.next(); # get then get the next one.
...etc. The class could be trimmed down a bit, but the main space overhead is in the _rconv list, which has the stored convolutions. That's roughly N*T/2, for O(NT) storage.
The convolutions only use the ranges, with a lot of randoms generated with the same constraints, the table construction time "amortizes away" to zero. The time complexity of .next() is roughly T/2 on average and O(T), in terms of the number of indexes into the _rconv lists.
To see how the algorithm works, assume a sequence of 3 zero-based choices, with max values (5,7,3), and a 0-based target T=10. Define or import the pulse and stepconv functions in an Idle session, then:
>>> pulse(5)
[1, 1, 1, 1, 1, 1]
>>> K1 = pulse (5)
>>> K2 = stepconv(K1, 7)
>>> K3 = stepconv(K2, 3)
>>> K1
[1, 1, 1, 1, 1, 1]
>>> K2
[1, 2, 3, 4, 5, 6, 6, 6, 5, 4, 3, 2, 1]
>>> K3
[1, 3, 6, 10, 14, 18, 21, 23, 23, 21, 18, 14, 10, 6, 3, 1]
>>> K3[10]
18
>>> sum(K3)
192
>>> (5+1)*(7+1)*(3+1)
192
K3[i] shows the number of different choice n_1, n_2, n_3 such that 0 <= n_k <= m_k and Σ n_k = i. Letting * mean convolution when applied to two of these lists. Then pulse(m_2)*pulse(m_3) is gives the distribution of sums of n_2 and n_3:
>>> R23 = stepconv(pulse(7),3)
>>> R23
[1, 2, 3, 4, 4, 4, 4, 4, 3, 2, 1]
>>> len(R23)
11
Every value from 0 to T=10 is (barely) possible, so any choice is possible for the first number and there are R23[T-n_1] possible triplets adding to T=10 that start with N1. So, once you've found that there are 18 possible sums adding to 10, generate a random number S = randint(18) and count down through the R23[T:T-m_1-1:-1] array:
>>> R23[10:10-5-1:-1]
[1, 2, 3, 4, 4, 4]
>>> sum(R23[10:10-5-1:-1])
18
Note the sum of that list is the total computed in K3[10] above. A sanity check. Anyway, if S==9 was the random choice, then find how many leading terms of that array can be summed without exceeding S. That's the value of n_1. In this case 1+2+3 <= S but 1+2+3+4 > S, so n_1 is 3.
As described above, you can then reduce the problem to find n_2. The final number (n_3 in this example) will be uniquely determined.