Memory Used to Subscript Arrays
QUESTION: Whoa! Does it, like, require a LOT of memory to subscript an array in IDL?
ANSWER: Indeed it does. Here is a recent discussion on the IDL newsgroup between Dick Jackson, who brought the topic to our attention, and Reimar Bauer. Dick, first:
This may be old news to some of you, but it surprised me and a couple of colleagues, and I couldn't find any discussion of it on this group, so I'll share it around.
I was surprised to find how much memory is used during access to a subset of an array. I ran this, which makes a 1000x1000 array, and accesses a subset of it using an array of subscripts. Here is a small test program to illustrate the situation:
PRO test a = bindgen(1000, 1000) subscripts = Long(RandomU(seed, 500)*1000) baseMem = (memory())[0] help, a[subscripts, *] highWaterMem = (memory())[3] Print, 'Memory used during access: ', highWaterMem-baseMem END IDL> test BYTE = Array[500, 1000] Memory used during access: 2500076The array being extracted is 0.5 million bytes, but it took 2.5 million bytes to do it! I'm guessing that there's a Long array being made behind the scenes that contains the indices of the elements I'm going to get back.
This came to light when my client was using a couple of 100MB images on a machine with 1GB of RAM, and my program ran out of memory! There are other factors here, but it was unexpected that accessing N bytes from an array requires 5*N bytes for the operation! I know now to be more careful.
This is the reply by Reimar, which has been edited slightly for clarity.
My last posting wasn't detailed enough to describe exactly what is going on. So here is another attempt. Consider a 1000 by 1000 byte array. Then the size of the array is 1MB.
a = bindgen(1000, 1000)If you now want to access this array by subscripts, then the subscripts are never passed by reference, they are passed by value and their size corresponds to their type. For example:
help, a[indgen(50)]This above expression consumes 50*2 + 50 bytes of memory during execution. And the expression below consumes 50*4 + 50 bytes of memory during execution.
help, a[indgen(50,/long)]Now what happens if we set one of the subscript values to an asterisk (*). By default, IDL will assign the index array a type of LONG. For example:
help, a[subscripts,*]is the same as this:
help, a[subscripts[0:499],0:999]This means the subscripts of A are of size 500 x 1000 x 4 and this is 2 MB of memory. Plus, you must add the result after the index operation in bytes of 500 x 1000 which gives 0.5MB. This is the measured 2.5 MB.
During the operation you need this amount of memory. I think what this means is that the array functionality of IDL is purchased at the cost of memory usage.
Error: Array Has Too Many Elements
Recently the amount of memory required for array subscripting was brought home to readers of the IDL newsgroup in a different way.
QUESTION: I have a problem. I have indices that are returned from performing a Where filter on a two-dimensional array. When I use those indices to subscript a similar 2D array, I have no problem. But if I use the same indices to subscript a frame of a 3D array, I get the error message % Array had too many elements. Here is an example of what I mean.
IDL> myArray = FltArr(3660, 1680) IDL> indices = Randomu(seed, 100000) * 3660L * 1680L IDL> myArray[indices, indices] = 5
No problems with this. But what about this?
IDL> DelVar, myArray
IDL> myArray = FltArr(3660, 1680, 8)
IDL> myArray[indices, indices, 4] = 5
% Array had too many elements
% Execution halted at: $MAIN$
Do you have any ideas on what could be going on here.
ANSWER: The answer was provided by JD Smith. He wrote:
I think the answer is pretty simple, if subtle. When IDL encounters a multi-dimensional subscript, it looks to see if all subscript vectors have the same dimensions. If they do, it "threads the list" and constructs indices from them on the fly as:
[vec1,vec2,vec3,...] ==> vec1+vec2*n1+vec3*n1*n2+...
where n1, n2, etc. are the sizes of the 1st, 2nd, etc. dimension of the array being indexed. That is, you have essentially specified a short list of index pairs, triples, etc., of length n_elements(vec1).
If, however, any of the subscripts are unspecified or zero-dimensional, by virtue of using a single index or one of the higher-order range operations (e.g. '*' or '0:5'), a temporary large array of indices has to be pre-created. Why? Because you can no longer "thread the list". For example, if I do this:
a[ [1,2,3] , [4,5,6], 0]
you might think I mean for IDL to generate a list like this:
a[1, 4, 0] a[2, 5, 0] a[3, 6, 0]
but it actually expands the expression to this:
a[1, 4, 0] a[1, 5, 0] a[1, 6, 0] a[2, 4, 0] a[2, 5, 0] a[2, 6, 0] a[3, 4, 0] a[3, 5, 0] a[3, 6, 0]
This can be an extremely important distinction when subscripting with large index vectors.
Here's an example demonstrating this:
IDL> a = RandomU(seed,100,100,100)
IDL> Help, /Memory
heap memory used: 4392171, max: 4392190, gets: 1895, frees: 1475
IDL> a[*,*,*] = 1
IDL> Help, /Memory
heap memory used: 4392203, max: 8392260, gets: 1899, frees: 1477
IDL> Print,(8392260-4392190)/4
1000017
Ah-ha, it seems a temporary index array of 100*100*100 indices was made. Make sense. What if we use three index vectors of the same size?
IDL> a = RandomU(seed,100,100,100)
IDL> r1 = RandomU(seed,100) & r2 = RandomU(seed,100) & r3 = RandomU(seed,100)
IDL> Help, /Memory
heap memory used: 4385413, max: 4385432, gets: 1212, frees: 807
IDL> a[r1,r2,r3] = 1
IDL> Help, /Memory
heap memory used: 4385448, max: 4386374, gets: 1218, frees: 811
In this case, a temporary index array was not needed; the three r vectors were used together directly as a threaded list, and no extra memory was used. How about this:
IDL> a = RandomU(seed,100,100,100)
IDL> Help, /Memory
heap memory used: 4383915, max: 4383934, gets: 1206, frees: 806
IDL> r1 = RandomU(seed,100) & r2 = RandomU(seed,100)
IDL> Help, /Memory
heap memory used: 4384920, max: 4384939, gets: 1211, frees: 807
IDL> a[r1,r2,*] = 1
IDL> Help, /Memory
heap memory used: 4384954, max: 8385481, gets: 1217, frees: 811
It seems 100*100*100 indices where created here too. Looks right. But now, let's stress things a bit:
IDL> a = RandomU(seed,100,100,100)
IDL> r1 = RandomU(seed,1000) & r2 = RandomU(seed,1000) & r3 = RandomU(seed,1000)
IDL> Help, /Memory
heap memory used: 4396199, max: 4396218, gets: 1210, frees: 806
IDL> a[r1,r2,r3] = 1
IDL> Help, /Memory
heap memory used: 4396234, max: 4404360, gets: 1216, frees: 810
Wait a minute, what's happening here? The subscript vectors, despite being larger than the dimensions of the array they're accessing, are still just being used directly, with no additional overhead required for creating a temporary index array. The assignment to 1 occurs 1000 times.
What about this:
IDL> a = RandomU(seed,100,100,100)
IDL> r1 = RandomU(seed,1000) & r2 = RandomU(seed,1000)
IDL> Help, /Memory
heap memory used: 4392105, max: 4392124, gets: 1209, frees: 806
IDL> a[r1,r2,*] = 1
IDL> Help, /Memory
heap memory used: 4400378, max: 404396266, gets: 1902, frees: 1478
IDL> Print, (404396266-4392124)/4
100001035
Uh oh. You can see that IDL had to pre-allocate a temporary index entry on the fly with 1000*1000*100 elements, despite the fact that it was used to index a much smaller array. The assignment to 1 occurs 100,000,000 times! Quite a difference. I can take this to the extreme:
IDL> a = RandomU(seed,1,1,1)
IDL> r1 = RandomU(seed,10000L) & r2 = RandomU(seed,10000L)
IDL> Help, /Memory
heap memory used: 464107, max: 464126, gets: 1209, frees: 806
IDL> a[r1,r2,0] = 1 ; long delay
IDL> Help, /Memory
heap memory used: 472380, max: 400504268, gets: 1902, frees: 1478
IDL> Print,(400504268-464126)/4
100010035
Oh my, nearly half a gigabyte of memory was allocated for the temporary index array just to assign a value to a single element (over and over again). What if I pre-build my index vector:
IDL> a = RandomU(seed,1,1,1)
IDL> r1 = RandomU(seed,10000L) & r2 = RandomU(seed,10000L)
IDL> r = r1 + 1*r2
IDL> Help, /Memory
heap memory used: 504193, max: 504212, gets: 1211, frees: 806
IDL> a[r] = 1 ; no delay
IDL> Help, /Memory
heap memory used: 504221, max: 544282, gets: 1215, frees: 808
What a difference this makes.
Bottom line? Keep in mind this duality in how IDL treats arrays as subscripts, and be very careful when mixing array subscripts with other types. If you mean for, for example,
[ [1,2], [3,4], 0 ] ==> [1,3,0], [2,4,0]
instead of this:
[ [1,2], [3,4], 0 ] ==> [1,3,0], [1,4,0], [2,3,0], [2,4,0]
Then you should use:
[ [1,2], [3,4], [0,0] ]
or just pre-build your indices as a single index vector beforehand.
Last Updated 8 January 2006