Computing all partitions of a vector

[John Mullahy]

Dear Colleagues: Is anyone aware of a Mata strategy for computing all possible partitions of a vector? A template of what I have in mind is

v=4,5,6
z=allpartitions(v)

where the z that is returned is the 5x3 matrix[

       1   2   3
    +-------------+
  1 |  1   1   1  |
  2 |  1   1   2  |
  3 |  1   2   1  |
  4 |  1   2   2  |
  5 |  1   2   3  |
    +-------------+

That is, the elements in each row of z indicate the partition membership index. In this example there are five partitions; in set notation

Partition 1: {4,5,6}
Partition 2: {4,5},{6}
Partition 3: {4,6},{5}
Partition 4: {4},{5,6}
Partition 5: {4},{5},{6}

The number of partitions of a given vector (=set) is given by Bell's Number https://oeis.org/A000110

It would be easy enough to brute force code this by hand for small sets, e.g. #S=3 or #S=4. But beyond this it would be nice to have a general function that did the heavy lifting.

Any suggestions or references would be appreciated. If I need to write something from scratch that's understandable, but if something already exists then all the better.

Thanks in advance.

[Mike Lacy]

I wonder if some of the stuff in Ben Jann's -moremata- might work here. There is a function -mmsubset()- within it that says it will "Obtain subsets, compositions, and partitions." I can't quickly tell if those things cover similar ground but describe it differently, but I think they are close enough to be useful. If you have that package installed, check out -help mmsubset()- .

 net describe moremata, from("http://fmwww.bc.edu/RePEc/bocode/m")

Hope this is not a wild goose chase.

[John Mullahy]

Thanks very much, Mike. I shall explore... it sounds promising.

[John Mullahy]

Unfortunately I didn't find exactly what I needed in –moremata– so decided to brute force a strategy. It's quite possible (indeed likely) that I'm the only one who'll ever find this useful, but I figured I'd close the loop by sharing what I did. (Did I mention brute force??)

Here's what I came up with, which involves defining four functions only one or two of which—pataa() and perhaps paacont(); patmat() and tga() operate behind the scenes—the user need explicitly engage.

cap mata mata drop patmat() tga() pataa() paacont()

set matastrict off

mata

/* patmat() computes all the unique partitions and returns what           */
/*  are essentially indexes that will be used to reference elements of v  */

function patmat(v) {
real ctch,bin,pmat,bind
 
ctch=J(1,cols(v),.)
pmat=J(rows(v),cols(v),.)

for (n=1;n<=rows(v);n++) {

bin=v[n,1]
ctch=1

for (j=2;j<=cols(v);j++) {
 bind=select(1..cols(bin),bin:==v[n,j])
 if (sum(bin:==v[n,j])>0) ctch=ctch,ctch[min(bind)]
 else ctch=ctch,(max(ctch)+1)
 bin=bin,v[n,j]
}

pmat[n,.]=ctch

}

return(uniqrows(pmat))

}



/* tga() returns a cols(v)^cols(v) matrix that is used as in input in pataa()  */


function tga(g) {
real matrix tv,tm,v

v=1..g
tv=tm=v

if (g==1) return(tv)

else if (g>1) {
for (j=1;j<=(g-1);j++) {
 tm=J(1,g,tm)
 tj=J(1,0,.)
 for (k=1;k<=g;k++) {
  tj=tj,J(1,g^j,v[k])
 }
 tm=tm\tj
}

return(sort(tm',1..cols(tm')))

}

}


/* pataa() is the function that computes the partitions given the output    */
/*  from patmat() and tga() and stores them in an associative array that is */
/*  returned by the function. Presumably the elements of this associative   */
/*  array will be retrieved to be used in subsequent compuatations of       */
/*  interest to the user.                                                   */
/*                                                                          */
/* The number of partitions of an m-vector v is given by B(m) where B(.) is */
/*  the Bell number https://oeis.org/A000110                                */


function pataa(v) {

real t,p,aa

aa = asarray_create("real", 2)

t=tga(cols(v))
p=patmat(t)

for (j=1;j<=rows(p);j++) {
 for (k=1;k<=max(p[j,.]);k++) {
   asarray(aa,(j,k),(select(v,p[j,.]:==k)))
 }
}

return(aa)

}



/* paacont() is a utility function (no pun intended) that displays the partitions  */

function paacont(aa) {
 real aak,saa
 aak=sort(asarray_keys(aa),(1,2))
 
 " "
 "(Note: For each partition each row represents one element of that partition)"
 
 for (j=1;j<=max(aak[.,1]);j++) {
  " "
  "Partition "+strofreal(j)
  saa=select(aak,aak[.,1]:==j)
  for (k=1;k<=max(saa[.,2]);k++) {
   st_matrix("z",asarray(aa,(j,k)))
   stata("matrix list z, noblank noheader nonames")
  }
 }
 
}


end

For example, to obtain all partitions of the row vector

v = (2,5,8,11)

and have these stored in an associative array z one summons pataa():

z = pataa(v)

.
Then to display the partitions:

paacont(z)

which displays

  
  (Note: For each partition each row represents one element of that partition)
   
  Partition 1
   2   5   8  11
   
  Partition 2
  2  5  8
  11
   
  Partition 3
   2   5  11
  8
   
  Partition 4
  2  5
   8  11
   
  Partition 5
  2  5
  8
  11
   
  Partition 6
   2   8  11
  5
   
  Partition 7
  2  8
   5  11
   
  Partition 8
  2  8
  5
  11
   
  Partition 9
   2  11
  5  8
   
  Partition 10
  2
   5   8  11
   
  Partition 11
  2
  5  8
  11
   
  Partition 12
   2  11
  5
  8
   
  Partition 13
  2
   5  11
  8
   
  Partition 14
  2
  5
   8  11
   
  Partition 15
  2
  5
  8
  11

To be useful in actual computations the user will likely need to retrieve the keys (indexes) from the associative array in order to retrieve the partitions themselves. This can be done by:

aak=asarray_keys(z)
aak=sort(aak,(1,2))
aak

which for this example returns

         1    2
     +-----------+
   1 |   1    1  |
   2 |   2    1  |
   3 |   2    2  |

         . . .

  36 |  15    3  |
  37 |  15    4  |
     +-----------+