If you want to keep all results, this can be done more efficiently using meshgrid. You can create N-dimensional grids of values, flatten them, then take the mean of the result.
So for the 2D case, meshgrid is pretty straightforward:
>>> a = np.arange(5)
>>> [a1, a2] = np.meshgrid(a, a)
>>> print(a1)
[[0 1 2 3 4]
[0 1 2 3 4]
[0 1 2 3 4]
[0 1 2 3 4]
[0 1 2 3 4]]
>>> print(a2)
[[0 0 0 0 0]
[1 1 1 1 1]
[2 2 2 2 2]
[3 3 3 3 3]
[4 4 4 4 4]]
>>> print(a1.ravel())
[0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4]
>>> print(a2.ravel())
[0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4]
As you can see, the flattened versions of a1 and a2 combined provide every possible combination. This can be extended to higher dimensions just by putting more a values:
>>> a = np.arange(3)
>>> [a1, a2, a3] = np.meshgrid(a, a, a)
>>> b = np.vstack([a1.ravel(), a2.ravel(), a3.ravel()])
>>> print(b)
[[0 0 0 1 1 1 2 2 2 0 0 0 1 1 1 2 2 2 0 0 0 1 1 1 2 2 2]
[0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2]
[0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2]]
This can be extended to any number of dimensions, allowing of course for memory constraints.
You can generalize the meshgrid function call by using unpacking. These two function calls are equivalent:
np.meshgrid(a, b, c)
np.meshgrid(*[a, b, c])
The * in this case means "unpack the list as function arguments" (also works with tuples, iteratirs, and many other things). You can combine this with the [x]*n syntax, where n is a number. This means "duplicate the contents of list [x] n times". So [5]*4 is similar to [5, 5, 5, 5] ("similar", not "identical", but the differences are not relevant to what we are doing).
So you can combine this to make calling meshgrid with multiple values easier. These three lines of code are the same:
np.meshgrid(a, a, a)
np.meshgrid([a, a, a])
np.meshgrid([a]*3)
You can exclude combinations where all the values are equal using logical indexing. First, you determine whether each element is equal to the value of the first element in that combination. Then you find combinations where that is True for all values. Then you get the columns for which that isn't the case. So something like:
>>> c = b[:, ~(b[0] == b).all(0)]
>>> print(c)
[[0 0 1 1 1 2 2 2 0 0 0 1 1 2 2 2 0 0 0 1 1 1 2 2]
[0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2]
[1 2 0 1 2 0 1 2 0 1 2 0 2 0 1 2 0 1 2 0 1 2 0 1]]
If you want to avoid memory issues you can make use of some internal knowledge of how numpy handles arrays. Arrays are simply continguous (or strided, i.e. continous in steps) blocks of memory. By working with this sort of memory, numpy can make use of very fast libraries for doing calculations on a block of memory at a time. But since this is simply a sequence of bytes, there are different ways to look at this block of memory. The same block could be seen as a 1D array of floats, a 2D array of integers, etc. Numpy can make use of this to avoid making unnecessary copies. So arr.ravel(), for example, "flattens" an array by treating the same underlying data as a 1D array. The catch (or benefit, depending on your use-case) is that changes to a raveled array will affect the original.
meshgrid can make use of this to avoid making unnecessary copies. Consider this:
>>> a = np.arange(10)
>>> b = np.meshgrid(*[a, a, a, a])
>>> c = np.meshgrid(*[a, a, a, a], copy=False)
>>> b[0].flags['OWNDATA']
True
>>> c[0].flags['OWNDATA']
False
flags['OWNDATA'] tells you whether the array is a view into the data stored in another array or not. In the case of b, it is, which means each b array takes up memory equal to the size of the array. In the case of c, however, it isn't, which means c takes almost no additional memory, it is still using the same block of memory as a. The catch is that any change to any element to any of the arrays in c will change the value in a and any of the values corresponding to that value in any of the c arrays.
If you combine this with ravel, you can get a view as well, but reshape(-1) is more reliable. Consider this:
>>> c[0].ravel().flags['OWNDATA']
True
>>> c[0].reshape(-1).flags['OWNDATA']
False
Slice of arrays, such as a[5:1000] also return views. "Fancy" indexing, where you get multiple elements of an array by index or boolean value, return new arrays. So you can slice out as many values as you want.
Getting the histogram of rows or columns of an array is a bit more difficult. The np.histogram can't do it. What you can do is follow the instructions here to get the unique rows (you need unique columns so just transpose), but when they use the unique function add the return_counts=True argument to also get the counts of each unique row. You will also need to sort each column using d.sort(axis=0) at the beginning so that columns with different orders of the same value are not counted separately.
Note that this doesn't use the highly-optimized numpy array processing libraries, so it will be very slow for large data set.
You can also read all the lines into a numpy array at once using loadtxt
So to implement what you want you can do:
from __future__ import print_function
import numpy as np
n = 4 # This is the number of repetitions
maxn = 10000
data = np.loadtxt('testdata.txt')
vals = np.meshgrid(*[data]*n, copy=False)
vals = np.vstack([val.reshape(-1)[:maxn] for val in vals])
vals = vals[:, ~(vals[0] == vals).all(0)]
p_avg = vals.mean(0)
print(p_avg)
Or as a function:
from __future__ import print_function
import numpy as np
def allmeans(data, n):
vals = np.meshgrid(*[data]*n)
vals = np.vstack([val.ravel() for val in vals])
vals = vals[:, ~(vals[0] == vals).all(0)]
return vals.mean(0)
data = np.loadtxt('testdata.txt')
print(allmeans(data, 4))
Doing variance or standard deviation just requires replacing mean with var and std, respectively. You could even pass it as a function:
from __future__ import print_function
import numpy as np
def allagg(data, n, func):
vals = np.meshgrid(*[data]*n)
vals = np.vstack([val.ravel() for val in vals])
vals = vals[:, ~(vals[0] == vals).all(0)]
return func(vals, 0)
data = np.loadtxt('testdata.txt')
print('Mean: \n', allagg(data, 4, np.mean))
print('Variance: \n', allagg(data, 4, np.var))
print('Standard Deviation: \n', allagg(data, 4, np.std))
NOTE: The differences between [x, x, x] and [x]*3 come into play when x is itself mutable (a data type that can change, like a list or numpy array). Consider this example:
>>> a = [[1], [1], [1]]
>>> b = [[1]]*3
>>> a[0][0] = 5
>>> b[0][0] = 5
>>> print(a)
[[5], [1], [1]]
>>> print(b)
[[5], [5], [5]]
This difference happens because in a each [1] is an independent list, while in b each [1] is the same list. Since you can change lists in-place, changes to a list in one part of your program affect that list everywhere else in the program. So these are more equivalent:
>>> x = [1]
>>> a = [x, x, x]
>>> b = [x]*3
>>> a[0][0] = 5
>>> print(a)
[[5], [5], [5]]
>>> print(b)
[[5], [5], [5]]
Notice I didn't even have top change b here, changes to a change every x used everywhere. This can be used to do funny things like create nested infinitely-nested lists:
>>> x = [1, 1]
>>> x[0] = x
>>> x
[[...], 1]
This issue affects lists of numpy arrays as well, but it doesn't matter to use because we aren't changing any of those arrays in-place.
EDIT 1: Change code to exclude duplicated
EDIT 2: Explain unpacking
EDIT 3: Warning about [x]*n syntax
EDIT 4: Added optimization tips and histogram
