Thread Subject:

Generating a huge array.

Hi,

I'm trying to generate an array with 4 columns (w1,w2,w3,w4) which is populated with every possible combination from 0 to 1, for each column, such that the row sum's to 1.

For example:
[0,0,0,0;
 0,0,0,0.001;
 0,0,0,0.002;
...
 0,0,0.001,0.001;
 0,0,0.002,0.001;
 0,0,0.003,0.001;
...
 0,0.001,0.001,0.001;
 0,0.002,0.001,0.001;
 0,0.003,0.001,0.001;
... and so on.]

Then perhaps loop through and remove all instances where w1,w2,w3,w4 don't sum to 1.

In Mathematica I would achieve this by doing this...
DeleteCases[
  Flatten[Table[{a, b, c, 1 - a - b - c}, {a, 0, 1, .001}, {b, 0,
     1, .001}, {c, 0, 1, .001}], 2],
  {_, _, _, x_ /; x < 0}];

As you can see, column 4 is generated by 1 minus the sum of the others, ensuring every column either sum's to 1, or is negative. Then the DeleteCases function, loops through this and removes any instances where the last column is negative, resulting in only those array's which sum to 1. Flatten is necessary to remove the multidimensional array.

I was considering finding a way to export the array from Mathematica and import it into MatLab, but there must be some way of doing this in MatLab, besides writing heaps of For statements and similar.

Just so you know, I'm generating every possible combination of weights for a portfolio of 4 assets. I know I can generate the efficient frontier, and have looked through the portfolio library and can't seem to find anything which will help.

Thanks.

"Uriah Stephenson-Ward" <uriah+spam@urijah.org> wrote in message <hpp13j$o1i$1@fred.mathworks.com>...
>
> I'm trying to generate an array with 4 columns (w1,w2,w3,w4) which is populated with every possible combination from 0 to 1, for each column, such that the row sum's to 1.
>
> For example:
> [0,0,0,0;
> 0,0,0,0.001;
> 0,0,0,0.002;
> ...
> 0,0,0.001,0.001;
> 0,0,0.002,0.001;
> 0,0,0.003,0.001;
> ...
> 0,0.001,0.001,0.001;
> 0,0.002,0.001,0.001;
> 0,0.003,0.001,0.001;
> ... and so on.]

So, you want to generate an array with 4 columns that has each combination of numbers from 0.000, 0.001, 0.002, ..., 0.999, 1.000, and then throw out those rows that don't sum to 1? As I look at your example and read your words it looks like you want to generate a 1000 x 1000 x 1000 array and then work with it to delete the rows you don't want. That's an 8 Gigabyte array you want to have in memory. Is that what you are asking or am I misunderstanding you? That's quite a large array. Do I understand your problem statement? You want to end up with an array like this:

[ 0.000, 0.000, 0.000, 1.000;
  0.000, 0.000, 0.001, 0.999;
               :
            etc.
  all combinations possible where rows sum to 1

Correct?

Because of the size involved, you may have to piece up the program in chunks.

James Tursa

"James Tursa" <aclassyguy_with_a_k_not_a_c@hotmail.com> wrote in message <hpp2rr$h1q$1@fred.mathworks.com>...
> "Uriah Stephenson-Ward" <uriah+spam@urijah.org> wrote in message <hpp13j$o1i$1@fred.mathworks.com>...
> >
> > I'm trying to generate an array with 4 columns (w1,w2,w3,w4) which is populated with every possible combination from 0 to 1, for each column, such that the row sum's to 1.
> >
> > For example:
> > [0,0,0,0;
> > 0,0,0,0.001;
> > 0,0,0,0.002;
> > ...
> > 0,0,0.001,0.001;
> > 0,0,0.002,0.001;
> > 0,0,0.003,0.001;
> > ...
> > 0,0.001,0.001,0.001;
> > 0,0.002,0.001,0.001;
> > 0,0.003,0.001,0.001;
> > ... and so on.]
>
> So, you want to generate an array with 4 columns that has each combination of numbers from 0.000, 0.001, 0.002, ..., 0.999, 1.000, and then throw out those rows that don't sum to 1? As I look at your example and read your words it looks like you want to generate a 1000 x 1000 x 1000 array and then work with it to delete the rows you don't want. That's an 8 Gigabyte array you want to have in memory. Is that what you are asking or am I misunderstanding you? That's quite a large array. Do I understand your problem statement? You want to end up with an array like this:
>
> [ 0.000, 0.000, 0.000, 1.000;
> 0.000, 0.000, 0.001, 0.999;
> :
> etc.
> all combinations possible where rows sum to 1
>
> Correct?
>
> Because of the size involved, you may have to piece up the program in chunks.
>
> James Tursa

Hi James,

That's correct. Perhaps instead of .001 increments, just .01 increments.

How could this be done in MatLab the easiest?

Thanks.

"Uriah Stephenson-Ward" <uriah+spam@urijah.org> wrote in message <hpp3ih$pp4$1@fred.mathworks.com>...
>
> That's correct. Perhaps instead of .001 increments, just .01 increments.
>
> How could this be done in MatLab the easiest?
>
> Thanks.

I'm not sure about easiest, but I tend to resort to mex routines for stuff like this. Below is the file. Call it weight1.c and generate the function like this:

mex -setup
(then press Enter and select a C compiler, such as lcc)
mex weight1.c

Then to generate your weights array do this (note the transpose at the end):

weights = weight1';

James Tursa

--------------------------------------------------------------------------------------------
weight1.c

#include "mex.h"
#define N 100
#define NDOUBLE 100.0
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    char i, j, k;
    char *sp, *sp0;
    double *pr;
    mwSize x, n;
    
    sp = sp0 = mxMalloc( N*N*N*N*sizeof(*sp) );
    for( i=0; i<N; i++ ) {
        for( j=0; j<N; j++ ) {
            for( k=0; k<N; k++ ) {
                if( i+j+k <= N ) {
                    *sp++ = i;
                    *sp++ = j;
                    *sp++ = k;
                    *sp++ = N - i - j - k;
                }
            }
        }
    }
    n = (sp - sp0) / 4;
    sp = sp0;
    plhs[0] = mxCreateDoubleMatrix(4, n, mxREAL);
    pr = mxGetPr(plhs[0]);
    for( x=0; x<4*n; x++ ) {
        *pr++ = *sp++ / NDOUBLE;
    }
    mxFree(sp0);
}

"Uriah Stephenson-Ward" <uriah+spam@urijah.org> wrote in message <hpp13j$o1i$1@fred.mathworks.com>...
> .........
> I'm trying to generate an array with 4 columns (w1,w2,w3,w4) which is populated with every possible combination from 0 to 1, for each column, such that the row sum's to 1.
> .........

  Doing all that removal seems very wasteful of time and computer memory. How about something like this:

 n = 1000; % Temporarily obtain sums w1+w2+w3+w4 = n
 M = zeros((n+1)*(n+2)*(n+3)/6,4);
 ix = 0;
 for w1 = 0:n
  for w2 = 0:n-w1
   for w3 = 0:n-w1-w2
    ix = ix+1;
    M(ix,:) = [w1,w2,w3,n-w1-w2-w3];
   end
  end
 end
 M = M/n; % Now the rows must sum to 1.

The value (n+1)*(n+2)*(n+3)/6 is the exact number of rows required. None have to be thrown away. I have done this using integers until the last step to avoid encountering round off errors with your decimal fractions.
 
Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hpp78k$ef6$1@fred.mathworks.com>...
> "Uriah Stephenson-Ward" <uriah+spam@urijah.org> wrote in message <hpp13j$o1i$1@fred.mathworks.com>...
> > .........
> > I'm trying to generate an array with 4 columns (w1,w2,w3,w4) which is populated with every possible combination from 0 to 1, for each column, such that the row sum's to 1.
> > .........
>
> Doing all that removal seems very wasteful of time and computer memory. How about something like this:
>
> n = 1000; % Temporarily obtain sums w1+w2+w3+w4 = n
> M = zeros((n+1)*(n+2)*(n+3)/6,4);
> ix = 0;
> for w1 = 0:n
> for w2 = 0:n-w1
> for w3 = 0:n-w1-w2
> ix = ix+1;
> M(ix,:) = [w1,w2,w3,n-w1-w2-w3];
> end
> end
> end
> M = M/n; % Now the rows must sum to 1.
>
> The value (n+1)*(n+2)*(n+3)/6 is the exact number of rows required. None have to be thrown away. I have done this using integers until the last step to avoid encountering round off errors with your decimal fractions.
>
> Roger Stafford

Excellent. Thanks for this guys.

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hpp78k$ef6$1@fred.mathworks.com>...
> "Uriah Stephenson-Ward" <uriah+spam@urijah.org> wrote in message <hpp13j$o1i$1@fred.mathworks.com>...
> > .........
> > I'm trying to generate an array with 4 columns (w1,w2,w3,w4) which is populated with every possible combination from 0 to 1, for each column, such that the row sum's to 1.
> > .........
>
> Doing all that removal seems very wasteful of time and computer memory. How about something like this:
>
> n = 1000; % Temporarily obtain sums w1+w2+w3+w4 = n
> M = zeros((n+1)*(n+2)*(n+3)/6,4);
> ix = 0;
> for w1 = 0:n
> for w2 = 0:n-w1
> for w3 = 0:n-w1-w2
> ix = ix+1;
> M(ix,:) = [w1,w2,w3,n-w1-w2-w3];
> end
> end
> end
> M = M/n; % Now the rows must sum to 1.
>
> The value (n+1)*(n+2)*(n+3)/6 is the exact number of rows required. None have to be thrown away. I have done this using integers until the last step to avoid encountering round off errors with your decimal fractions.
>
> Roger Stafford

Interestingly, the Mathematica script creates an array <174504x4 double> and this script creates an array <176851x4 double>, I'm wondering if Mathematica was dropping some possible values, or vice a versa.

Thanks.

After looking at Roger's post I realized I had a typo in my code. The for loop tests were ending prematurely and need <= tests instead of < tests.

#include "mex.h"
#define N 100
#define NDOUBLE 100.0
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    char i, j, k;
    char *sp, *sp0;
    double *pr;
    mwSize x, n;
    
    sp = sp0 = mxMalloc( N*N*N*N*sizeof(*sp) );
    for( i=0; i<=N; i++ ) {
        for( j=0; j<=N; j++ ) {
            for( k=0; k<=N; k++ ) {
                if( i+j+k <= N ) {
                    *sp++ = i;
                    *sp++ = j;
                    *sp++ = k;
                    *sp++ = N - i - j - k;
                }
            }
        }
    }
    n = (sp - sp0) / 4;
    sp = sp0;
    plhs[0] = mxCreateDoubleMatrix(4, n, mxREAL);
    pr = mxGetPr(plhs[0]);
    for( x=0; x<4*n; x++ ) {
        *pr++ = *sp++ / NDOUBLE;
    }
    mxFree(sp0);
}

"Uriah Stephenson-Ward" <uriah+spam@urijah.org> wrote in message <hpp8ip$2cs$1@fred.mathworks.com>...
>
> Interestingly, the Mathematica script creates an array <174504x4 double> and this script creates an array <176851x4 double>, I'm wondering if Mathematica was dropping some possible values, or vice a versa.

My corrected code and Roger's code (with n = 100) both independently come up with 176851, so I would conclude that the 174504 number from the Mathematica script is wrong.

James Tursa

I have such function on FEX:
http://www.mathworks.com/matlabcentral/fileexchange/17818-all-permutations-of-integers-with-sum-criteria

>> AllVL1(4,100,'',NaN) % Number of solution

ans =

      176851

>> v=AllVL1(4,100)/100;
>> size(v)

ans =

      176851 4

>> v(1:10,:) % Check first 10 solutions

ans =

         0 0 0 1.0000
         0 0 0.0100 0.9900
         0 0 0.0200 0.9800
         0 0 0.0300 0.9700
         0 0 0.0400 0.9600
         0 0 0.0500 0.9500
         0 0 0.0600 0.9400
         0 0 0.0700 0.9300
         0 0 0.0800 0.9200
         0 0 0.0900 0.9100

>> AllVL1(4,1000,'',NaN) % For 0.001 this is the number of solutions

ans =

   167668501

>> Gb=167668501*4*8/1024^3 % The size in Gbytes of such array

Gb =

    4.9969

% Bruno

"Uriah Stephenson-Ward" <uriah+spam@urijah.org> wrote in message <hpp8ip$2cs$1@fred.mathworks.com>...
> ........
> Interestingly, the Mathematica script creates an array <174504x4 double> and this script creates an array <176851x4 double>, I'm wondering if Mathematica was dropping some possible values, or vice a versa.
>
> Thanks.

  Offhand I would guess that Mathematica suffered from some improper round offs with the decimal fractions you gave it. With your inequality x < 0, you may have thrown away some combinations whose x was just a tiny round off amount below zero and that therefore ought to have been saved. Perhaps if you had said x < -0.00000001, it might have worked properly. Remember, on a binary machine (which we all have nowadays) the fractions 0.001 and 0.01 and most other decimal fractions cannot be exactly represented as binary fractions and this can cause mischief. That is no fault of Mathemetica, but rather the algorithm and the machine used. With n equal to 100 you should have gotten precisely 101*102*103/6 = 176851, as James has said. It is something that can be proved with mathematical rigor.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hppeh1$jdo$1@fred.mathworks.com>...
> "Uriah Stephenson-Ward" <uriah+spam@urijah.org> wrote in message <hpp8ip$2cs$1@fred.mathworks.com>...
> > ........
> > Interestingly, the Mathematica script creates an array <174504x4 double> and this script creates an array <176851x4 double>, I'm wondering if Mathematica was dropping some possible values, or vice a versa.
> >
> > Thanks.
>
> Offhand I would guess that Mathematica suffered from some improper round offs with the decimal fractions you gave it. With your inequality x < 0, you may have thrown away some combinations whose x was just a tiny round off amount below zero and that therefore ought to have been saved. Perhaps if you had said x < -0.00000001, it might have worked properly. Remember, on a binary machine (which we all have nowadays) the fractions 0.001 and 0.01 and most other decimal fractions cannot be exactly represented as binary fractions and this can cause mischief. That is no fault of Mathemetica, but rather the algorithm and the machine used. With n equal to 100 you should have gotten precisely 101*102*103/6 = 176851, as James has said. It is something that can be proved with mathematical rigor.
>
> Roger Stafford

This is an interesting discussion, in particular since Wolfram promotes Mathematica on the strength of its symbolic mathematics capabilities (not sniping at Wolfram here, just stating the obvious). If you change the sum to 100 and the step size to 1 does Mathematica still return the incorrect answer?

"Mark Shore" <mshore@magmageosciences.ca> wrote in message <hpqb7s$64m$1@fred.mathworks.com>...
> This is an interesting discussion, in particular since Wolfram promotes Mathematica on the strength of its symbolic mathematics capabilities (not sniping at Wolfram here, just stating the obvious). If you change the sum to 100 and the step size to 1 does Mathematica still return the incorrect answer?
----------------
Hi Mark,

  That was the only way I could account for Mathematica not getting the correct number of surviving entries in its 'Table' - by doing at least some of its calculations in binary numeric form rather than symbolic, which would of necessity be subject to round off error. I don't really know much about Mathematica, so that observation has to be considered speculative. I have to assume that with all integral values in its 'Table' it really ought to get the count right. It's an experiment Uriah could presumably carry out.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hpqd9g$1rf$1@fred.mathworks.com>...
> "Mark Shore" <mshore@magmageosciences.ca> wrote in message <hpqb7s$64m$1@fred.mathworks.com>...
> > This is an interesting discussion, in particular since Wolfram promotes Mathematica on the strength of its symbolic mathematics capabilities (not sniping at Wolfram here, just stating the obvious). If you change the sum to 100 and the step size to 1 does Mathematica still return the incorrect answer?
> ----------------
> Hi Mark,
>
> That was the only way I could account for Mathematica not getting the correct number of surviving entries in its 'Table' - by doing at least some of its calculations in binary numeric form rather than symbolic, which would of necessity be subject to round off error. I don't really know much about Mathematica, so that observation has to be considered speculative. I have to assume that with all integral values in its 'Table' it really ought to get the count right. It's an experiment Uriah could presumably carry out.
>
> Roger Stafford

It looks like a good guess. I have Mathematica on another computer (almost never use it - thought I'd be doing some symbolic math for my own edification, but didn't...), so I fired it up and ran Uriah's script for a range of values that _should_ all produce the same answer:

1000/10 (176851,4)
100/1 (176851,4)
10/0.1 (172153,4)
1/0.01 (173809,4)
0.1/0.001 (176263,4)

Also interesting is the fact that none of the array sizes generated match that from Uriah, suggesting that this could alternatively be a CPU or OS related bug. I'm using v7.0.1 on Windows XP (32 bit) with a T9400 Core 2 Duo mobile processor.

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hppeh1$jdo$1@fred.mathworks.com>...
> "Uriah Stephenson-Ward" <uriah+spam@urijah.org> wrote in message <hpp8ip$2cs$1@fred.mathworks.com>...
> > ........
> > Interestingly, the Mathematica script creates an array <174504x4 double> and this script creates an array <176851x4 double>, I'm wondering if Mathematica was dropping some possible values, or vice a versa.
> >
> > Thanks.
>
> Offhand I would guess that Mathematica suffered from some improper round offs with the decimal fractions you gave it. With your inequality x < 0, you may have thrown away some combinations whose x was just a tiny round off amount below zero and that therefore ought to have been saved. Perhaps if you had said x < -0.00000001, it might have worked properly. Remember, on a binary machine (which we all have nowadays) the fractions 0.001 and 0.01 and most other decimal fractions cannot be exactly represented as binary fractions and this can cause mischief. That is no fault of Mathemetica, but rather the algorithm and the machine used. With n equal to 100 you should have gotten precisely 101*102*103/6 = 176851, as James has said. It is something that can be proved with mathematical rigor.
>
> Roger Stafford

Roger is correct.

Dimensions[
 DeleteCases[
  Flatten[Table[{a, b, c, 1 - a - b - c}, {a, 0, 1, .01}, {b, 0,
     1, .01}, {c, 0, 1, .01}], 2], {_, _, _, x_ /; x < -0.000001}]]

returns (176851,4)

So a basic lesson in round-off errors and relational operators is relearned. It just looked different in Mathematica, I guess.