Partitions of multisets (oeis.org/A249620)

<?php

// en.wikiversity.org/wiki/Partitions_of_multisets
// oeis.org/A249620

// Rather quick and dirty script to create the partitions of multisets with up to 8 elements.
// The essential steps are the creation of the partitions of {1..n} from the values in oeis.org/A231428,
// the replacement of 1..n with the elements of the multiset using preg_replace
// and the removal of redundant rows in the new array.
// There are certainly better ways to do this.
// Tilman Piesk, November 2014


// the numbers n and k determine the output
// n is the number of elements in the multiset, and k its signature
$n = 5;  // values between 1 and 8
$k = 3;  // values between 0 and $partnum[$n]-1

// the numbers of partitions of n (oeis.org/A000041)
$partnum = array(1, 1, 2, 3, 5, 7, 11, 15, 22);

// the partitions of an 8-element set (oeis.org/A231428)
// $setpartsDec = array(0,1,2,4,7,8,12,16,18,25,32,33,42,52,63,64,68,80,96,116,128,130,136,160...268435455);

// non-one addends of the partitions of 8 in the order defined by oeis.org/A194602
$intparts = array(
    array(),
    array(2),
    array(3),
    array(2, 2),
    array(4),
    array(3, 2),
    array(5),
    array(2, 2, 2),
    array(4, 2),
    array(3, 3),
    array(6),
    array(3, 2, 2),
    array(5, 2),
    array(4, 3),
    array(7),
    array(2, 2, 2, 2),
    array(4, 2, 2),
    array(3, 3, 2),
    array(6, 2),
    array(5, 3),
    array(4, 4),
    array(8)
);


// $setparts will contain the set partitions in this form: https://oeis.org/A231428/a231428.txt
$setparts = array();

// the binary digits go from $setpartBin to $setpartMat using these coordinates
$coord = array(
    array(1,0),
    array(2,0),array(2,1),
    array(3,0),array(3,1),array(3,2),
    array(4,0),array(4,1),array(4,2),array(4,3),
    array(5,0),array(5,1),array(5,2),array(5,3),array(5,4),
    array(6,0),array(6,1),array(6,2),array(6,3),array(6,4),array(6,5),
    array(7,0),array(7,1),array(7,2),array(7,3),array(7,4),array(7,5),array(7,6)
);

for ($a=0; $a<4140; $a++) {

    // decimal from A231428 to binary string
    $setpartBin = strrev(decbin($setpartsDec[$a]));

    // this is going to be the equivalence relation matrix of the set partition
    $setpartMat = array(
        array(1, 0, 0, 0, 0, 0, 0, 0),
        array(0, 1, 0, 0, 0, 0, 0, 0),
        array(0, 0, 1, 0, 0, 0, 0, 0),
        array(0, 0, 0, 1, 0, 0, 0, 0),
        array(0, 0, 0, 0, 1, 0, 0, 0),
        array(0, 0, 0, 0, 0, 1, 0, 0),
        array(0, 0, 0, 0, 0, 0, 1, 0),
        array(0, 0, 0, 0, 0, 0, 0, 1)
    );

    // get the 1s from $setpartBin and put them in $setpartMat
    for ($b = 0; $b < strlen($setpartBin); $b++) {
        if ($setpartBin[$b] == '1') {
            $i = $coord[$b][0];
            $j = $coord[$b][1];
            $setpartMat[$i][$j] = 1;
            $setpartMat[$j][$i] = 1;
        }
    }

    // remove redundant rows from $setpartMat
    $setpartMat = array_unique($setpartMat, SORT_REGULAR);

    $blocks = array();

    foreach ($setpartMat as $matrow) {
        $block = array();
        for ($b = 0; $b < 8; $b++) {
            if ($matrow[$b] == 1) {
                array_push($block, $b + 1);  // the elements are numbered from 1
            }
        }
        if (count($block) > 1) {  // only non-singleton blocks shall be shown
            array_push($blocks, $block);
        }
    }

    array_push($setparts, $blocks);

}
// $setparts is ready


////////////////////////////////////////////////////////////////////////////////////////////////


// the numbers of partitions of n-sets (oeis.org/A000110)
$bells = array(1, 1, 2, 5, 15, 52, 203, 877, 4140);
$bell = $bells[$n];

// $setpartsBell is an array of the first $bell set partitions
$setpartsBell = array_slice($setparts, 0, $bell);

foreach ($setpartsBell as $b => $setpart) {
    $merge = array();
    foreach ($setpart as $block) {
        $merge = array_merge($merge, $block);
    }
    $singletons= array_diff( range(1, $n), $merge);
    foreach ($singletons as $singleton) {
        array_push($setpart, $singleton);
    }
    sort($setpart);
    $setpartsBell[$b] = $setpart;
}
$setpartsBellStr = var_export($setpartsBell, true);
$setpartsBellStr = preg_replace('/\s/', '', $setpartsBellStr);  // remove white space
$setpartsBellStr = preg_replace('/\d+=>/', '', $setpartsBellStr);  // remove indices


// the signature of the multiset is the $k-th element of $intparts
// it is in the form $intpartAddends = array(3,2) and must be transformed in the form $intpartMultiset = 1,1,1,2,2,3,4,5 (lengh $n)
$intpartAddends = $intparts[$k];
$intpartMultiset = array();
$count = 1;
foreach($intpartAddends as $addend) {
    for ($i=1; $i<=$addend; $i++) {
        array_push($intpartMultiset, $count);
    }
    $count++;
}
// elements for the non-one addends must be added to get length $n
while (count($intpartMultiset)<$n) {
    array_push($intpartMultiset, $count++);
}


// $setpartsBellStr is $setpartsBell as a string, i.e. it is an array that contains partitions of the set {1..n}
// now the numbers 1..n shall be replaced by the corresponding elements of $intpartMultiset to get the partitions of the multiset
$count = 1;
foreach ($intpartMultiset as $val) {
    $setpartsBellStr = preg_replace("/$count/", "$val", $setpartsBellStr);
    $count++;
}

// create the array of partitions of the multiset from the manipulated string
eval( '$multisetparts =' . $setpartsBellStr . ';' );

// remove redundant partitions
foreach ($multisetparts as $key => $val) {
    $derp = $val;
    sort($derp);
    $multisetparts[$key] = $derp;
}

sort($multisetparts);
$multisetparts = array_map("unserialize", array_unique(array_map("serialize", $multisetparts)));


////////////////////////////////////////////////////////////////////////////////////////////////


// print the content of $multisetparts on the screen
// for n=5 and k=3 this looks like http://pastebin.com/g5AJBRgT
$outputStr = "elements in multiset: $n         integer partition: $k (" . implode('+',$intpartAddends);
if (array_sum($intpartAddends)<$n) {
    if ($k==0) {
        $outputStr .= '1' . str_repeat('+1', $n - array_sum($intpartAddends) - 1);
    } else {
        $outputStr .= str_repeat('+1', $n - array_sum($intpartAddends));
    }
}
$outputStr .= ")\n\n";

$i = 1;
$mspLen = count($multisetparts);
$strPadLen = strlen((string)($mspLen));
foreach ($multisetparts as $part) {
    $outputStr .= str_pad($i, $strPadLen, ' ', STR_PAD_LEFT) . ':    ';
    $partLen = count($part);
    foreach ($part as $block) {
        if (!is_array($block)) {
            $block = array($block);
        }
        $outputStr .= '{ ' . implode(' ', $block) . ' }';
    }
    if ($i < $mspLen) {
        $outputStr .= "\n";
    }
    $i++;
}

echo $outputStr;