<?php
/**
 * @file
 *   Random set of numbers with a fixed sum. Much of this code is a PHP port of
 *   the original Matlab code, as noted.
 *
 *   Author: Jeremy Epstein
 *   Date: 24 Sep 2008
 *   License: GNU General Public License v2
 *   http://www.gnu.org/licenses/gpl-2.0.txt
 */

print_random_set(20, 60, 1, 8);

/**
 * Generates a set of random numbers, such that all numbers are within
 * a given range, and such that the sum of all numbers is equal to
 * a user-defined value.
 *
 * This is a PHP port of the randfixedsum function, written in Matlab
 * by Roger Stafford, and with minimal modifications from the original
 * code. Randfixedsum can be found here:
 * http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=9700&objectType=file
 *
 * @param $size
 *   Number of elements in the set.
 * @param $sum
 *   Total value that the sum of all the elements must equal.
 * @param $min
 *   Minimum value of a given element.
 * @param $max
 *   Maximum value of a given element.
 *
 * @return
 *   Array of floating-point random numbers, or an error message on failure.
 */
function randfixedsum($size, $sum, $min, $max) {
  // Check the arguments
  if ($size < 1) {
    return "randfixedsum was passed arguments that violate the rule: size < 1";
  }
  if ($sum < $size * $min) {
    return "randfixedsum was passed arguments that violate the rule: sum < size * min";
  }
  if ($sum > $size * $max) {
    return "randfixedsum was passed arguments that violate the rule: sum > size * max";
  }
  if ($min >= $max) {
    return "randfixedsum was passed arguments that violate the rule: min >= max";
  }

  // Rescale to a unit cube: 0 <= x(i) <= 1
  $sum = ($sum - $size * $min) / ($max - $min);

  // Construct the transition probability table, t.
  // t(i,j) will be utilized only in the region where j <= i + 1.

  // Must have 0 <= k <= size - 1
  $k = max(min(floor($sum), ($size - 1.0)), 0.0);
  // Must have k <= sum <= k + 1
  $sum = max(min($sum, $k + 1.0), $k);

  // s1 & s2 will never be negative
  $s1 = array();
  $s2 = array();
  foreach (range($k, $k - $size + 1.0, -1.0) as $s1_element) {
    $s1[] = $sum - $s1_element;
  }
  foreach (range($k + $size, $k + 1.0, -1.0) as $s2_element) {
    $s2[] = $s2_element - $sum;
  }

  // Scale for full 'double' range
  $w = array();
  for ($i = 0; $i < $size; $i++) {
    $w[] = array_fill(0, $size + 1, 0.0);
  }
  $massive = 1.7976931348623157E+308;
  $w[0][1] = $massive;

  $t = array();
  for ($i = 0; $i < $size - 1; $i++) {
    $t[] = array_fill(0, $size, 0.0);
  }

  // The smallest positive matlab 'double' no.
  $tiny = pow(2.0, -1074.0);

  for ($i = 1; $i < $size; $i++) {
    $tmp1 = array();
    $tmp2 = array();
    $tmp3 = array();
    $tmp4 = array();

    $sizeIndex = (int)($size - ($i + 1));
    for ($j = 0; $j <= $i; $j++) {
      $tmp1[] = $w[$i - 1][$j + 1] * ($s1[$j] / ((double)$i + 1.0));
      $tmp2[] = $w[$i - 1][$j] * ($s2[$sizeIndex] / ((double)$i + 1.0));
      $sizeIndex++;
    }

    for ($j = 0; $j <= $i; $j++) {
      $w[$i][$j + 1] = $tmp1[$j] + $tmp2[$j];
    }

    $sizeIndex = (int)$size - ($i + 1);
    for ($j = 0; $j <= $i; $j++) {
      // In case tmp1 & tmp2 are both 0,
      $tmp3[] = $w[$i][$j + 1] + $tiny;
      // then t is 0 on left & 1 on right
      $tmp4[] = $s2[$sizeIndex] > $s1[$j] ? 1.0 : 0.0;
      $t[$i - 1][$j] = ($tmp2[$j] / $tmp3[$j]) * $tmp4[$j] + (1.0 - $tmp1[$j] / $tmp3[$j]) * ($tmp4[$j] == 0.0 ? 1.0 : 0.0);
      $sizeIndex++;
    }
  }

  // Derive the polytope volume v from the appropriate
  // element in the bottom row of w.
  $v = pow($size, 3.0 / 2.0) * ($w[(int)$size - 1][(int)$k + 1] / $massive) * pow($max - $min, $size - 1);

  // Now compute the matrix x.
  $x = array_fill(0, $size, 0.0);
  $rt = array();
  $rs = array();

  for ($i = 0; $i < $size - 1; $i++) {
    // For random selection of simplex type
    $rt[] = (float)rand() / (float)getrandmax();
    // For random location within a simplex
    $rs[] = (float)rand() / (float)getrandmax();
  }

  // Start with sum zero & product 1
  $sm = 0.0;
  $pr = 1.0;
  $jj = 0;

  // Work backwards in the t table
  for ($i = (int)$size - 1; $i > 0; $i--) {
    $e = NULL;
    $sx = NULL;
    // Use rt to choose a transition
    $e = $rt[$jj] <= $t[$i - 1][(int)$k] ? 1.0 : 0.0;
    // Use rs to compute next simplex coord.
    $sx = pow($rs[$jj], 1.0 / (double)$i);
    // Update sum
    $sm += (1.0 - $sx) * $pr * $sum / (double)($i + 1);
    // Update product
    $pr *= $sx;
    // Calculate x using simplex coords.
    $x[$jj] = $sm + $pr * $e;
    // Transition adjustment
    $sum -= $e;
    $k -= $e;

    $jj++;
  }

  // Compute the last x
  $x[(int)$size - 1] = $sm + $pr * $sum;

  // Randomly permute the order in the columns of x and rescale.
  $p = array();
  for ($i = 0; $i < $size; $i++) {
    $p[] = $i;
  }

  // Use p to carry out a matrix 'randperm'
  shuffle($p);

  $ret = array();
  for ($i = 0; $i < $size; $i++) {
    // Permute & rescale x
    $ret[] = (double)($max - $min) * (double)$x[$p[$i]] + (double)$min;
  }

  return $ret;
}

/**
 * Converts the floating-point numbers returned by the randfixedsum
 * algorithm into rounded integers. Adds or subtracts any disrepancy
 * from the total to random numbers in the set, as needed.
 *
 * @param $size
 *   Number of elements in the set.
 * @param $sum
 *   Total value that the sum of all the elements must equal.
 * @param $min
 *   Minimum value of a given element.
 *@param $max
 *   Maximum value of a given element.
 */
function randfixedsum_rounded($size, $sum, $min, $max) {
  $list_dbl = randfixedsum((double)$size, (double)$sum, (double)$min, (double)$max);
  $list_int = array();
  $sum_of_rounded = 0;

  if (!is_array($list_dbl)) {
    return $list_dbl;
  }

  foreach ($list_dbl as $rand_dbl) {
    $rand_int = (int)round($rand_dbl);
    $list_int[] = $rand_int;
    $sum_of_rounded += $rand_int;
  }

  // Randomly permute the order in the columns of list_int and rescale.
  $p = array();
  for ($i = 0; $i < $size; $i++) {
    $p[] = $i;
  }

  // Use p to carry out a matrix 'randperm'
  shuffle($p);

  // The total of the rounded integers is lower than the total
  // that we want.
  if ($sum_of_rounded < $sum) {
    // Start by searching for an element with value 'min', and
    // increase if unsuccessful.
    $seek_amt = $min;
    while ($seek_amt < $max && $sum_of_rounded < $sum) {
      // Try and raise the rounded total to equal the required
      // total, by adding onto a single element in the set.
      $add_amt = $sum - $sum_of_rounded;
      $is_added = FALSE;

      // No single element in the set can exceed max. In cases
      // where adding (sum - $sum_of_rounded) would violate this,
      // add as much as is permitted, and leave the rest for
      // another element.
      if ($seek_amt + $add_amt > $max) {
        $add_amt = $max - $seek_amt;
      }

      for ($i = 0; $i < $size; $i++) {
        if ($list_int[$p[$i]] == $seek_amt) {
          $list_int[$p[$i]] += $add_amt;
          $sum_of_rounded += $add_amt;
          $is_added = TRUE;
          break;
        }
      }

      if (!$is_added) {
        $seek_amt++;
      }
    }
  }
  // Or the total of the rounded integers is higher than the total
  // that we want.
  elseif ($sum_of_rounded > $sum) {
    // Start by searching for an element with value 'max', and
    // decrease if unsuccessful.
    $seek_amt = $max;
    while ($seek_amt > $min && $sum_of_rounded > $sum) {
      // Try and lower the rounded total to equal the required
      // total, by subtracting from a single element in the set.
      $subtract_amt = $sum_of_rounded - $sum;
      $isSubtracted = FALSE;

      // No single element in the set can be lower than min.
      // In cases where subtracting ($sum_of_rounded - sum)
      // would violate this, subtract as much as is permitted,
      // and leave the rest for another element.
      if ($seek_amt - $subtract_amt < $min) {
        $subtract_amt = $seek_amt - $min;
      }

      for ($i = 0; $i < $size; $i++) {
        if ($list_int[$p[$i]] == $seek_amt) {
          $list_int[$p[$i]] -= $subtract_amt;
          $sum_of_rounded -= $subtract_amt;
          $isSubtracted = TRUE;
          break;
        }
      }

      if (!$isSubtracted) {
        $seek_amt--;
      }
    }
  }

  return $list_int;
}

/**
 * Prints a set of random numbers formatted in a table.
 */
function print_random_set($size, $sum, $min, $max) {
?>
<html>
<head>
  <title>My numbers</title>
</head>

<body>

<?php
  $random_set = randfixedsum_rounded($size, $sum, $min, $max);
  if (is_array($random_set)) {
?>
<h1>My numbers</h1>

<table border="1">
  <tr>
    <th rowspan="<?php print $size; ?>" valign="top">Set</th>
<?php
    $total = 0.0;
    $first = TRUE;

    foreach ($random_set as $random_num) {
      if ($first) {
        print "    <td>". $random_num ."</td>\n  </tr>\n";
        $first = FALSE;
      }
      else {
        print "  <tr>\n    <td>". $random_num ."</td>\n  </tr>\n";
      }
      $total += $random_num;
    }

    print "  <tr>\n    <th>Total</td>\n    <th>". $total ."</th>\n  </tr>\n";
?>
</table>
<?php
  }
  else {
    print "<p>Error: ". $random_set ."</p>\n";
  }
?>

</body>
</html>
<?php
}