// 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

using System;
using System.Collections.Generic;
using System.Text;

namespace RandomNess
{
    /// <summary>
    /// Random number utility class.
    /// </summary>
    class RandomNess
    {
        public Random random;

        public RandomNess()
        {
            random = new Random();
        }

        /// <summary>
        /// Generates a random number between min and max.
        /// </summary>
        public int RandomNumber(int min, int max)
        {
            return random.Next(min, max);
        }

        /// <summary>
        /// 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 C# 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
        /// </summary>
        /// <param name="size">Number of elements in the set.</param>
        /// <param name="sum">Total value that the sum of all the elements must
        /// equal.</param>
        /// <param name="min">Minimum value of a given element.</param>
        /// <param name="max">Maximum value of a given element.</param>
        public List<double> RandFixedSum(double size, double sum, double min, double max)
        {
            // Check the arguments
            if (size < 1)
            {
                throw new System.ArgumentException("RandomNumber.RandFixedSum was passed arguments that violate the rule: size < 1");
            }
            if (sum < size * min)
            {
                throw new System.ArgumentException("RandomNumber.RandFixedSum was passed arguments that violate the rule: sum < size * min");
            }
            if (sum > size * max)
            {
                throw new System.ArgumentException("RandomNumber.RandFixedSum was passed arguments that violate the rule: sum > size * max");
            }
            if (min >= max)
            {
                throw new System.ArgumentException("RandomNumber.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
            double k = Math.Max(Math.Min(Math.Floor(sum), (size - 1f)), 0f);
            // Must have k <= sum <= k + 1
            sum = Math.Max(Math.Min(sum, k + 1f), k);

            // s1 & s2 will never be negative
            List<double> s1 = new List<double>(), s2 = new List<double>();
            foreach (double s1Element in new Range(k, k - size + 1f, -1f))
            {
                s1.Add(sum - s1Element);
            }
            foreach (double s2Element in new Range(k + size, k + 1f, -1f))
            {
                s2.Add(s2Element - sum);
            }

            // Scale for full 'double' range
            List<List<double>> w = new List<List<double>>();
            for (int i = 0; i < size; i++)
            {
                w.Add(ListFiller.GetList(size + 1, 0f));
            }
            w[0][1] = double.MaxValue;

            List<List<double>> t = new List<List<double>>();
            for (int i = 0; i < size - 1; i++)
            {
                t.Add(ListFiller.GetList(size, 0f));
            }

            // The smallest positive matlab 'double' no.
            double tiny = Math.Pow(2f, -1074f);

            for (int i = 1; i < size; i++)
            {
                List<double> tmp1 = new List<double>(i + 1), tmp2 = new List<double>(i + 1), tmp3 = new List<double>(i + 1), tmp4 = new List<double>(i + 1);

                int sizeIndex = (int)size - (i + 1);
                for (int j = 0; j <= i; j++)
                {
                    tmp1.Add(w[i - 1][j + 1] * (s1[j] / ((double)i + 1f)));
                    tmp2.Add(w[i - 1][j] * (s2[sizeIndex] / ((double)i + 1f)));
                    sizeIndex++;
                }

                for (int j = 0; j <= i; j++)
                {
                    w[i][j + 1] = tmp1[j] + tmp2[j];
                }

                sizeIndex = (int)size - (i + 1);
                for (int j = 0; j <= i; j++)
                {
                    // In case tmp1 & tmp2 are both 0,
                    tmp3.Add(w[i][j + 1] + tiny);
                    // then t is 0 on left & 1 on right
                    tmp4.Add(s2[sizeIndex] > s1[j] ? 1f : 0f);
                    t[i - 1][j] = (tmp2[j] / tmp3[j]) * tmp4[j] + (1f - tmp1[j] / tmp3[j]) * (tmp4[j] == 0f ? 1f : 0f);
                    sizeIndex++;
                }
            }

            // Derive the polytope volume v from the appropriate
            // element in the bottom row of w.
            double v = Math.Pow(size, 3f / 2f) * (w[(int)size - 1][(int)k + 1] / Double.MaxValue) * Math.Pow(max - min, size - 1);

            // Now compute the matrix x.
            List<double> x = ListFiller.GetList(size, 0f);
            List<double> rt = new List<double>((int)size - 1);
            List<double> rs = new List<double>((int)size - 1);

            for (int i = 0; i < size - 1; i++)
            {
                // For random selection of simplex type
                rt.Add(random.NextDouble());
                // For random location within a simplex
                rs.Add(random.NextDouble());
            }

            // Start with sum zero & product 1
            double sm = 0f, pr = 1f;
            int jj = 0;

            // Work backwards in the t table
            for (int i = (int)size - 1; i > 0; i--)
            {
                double e, sx;
                // Use rt to choose a transition
                e = rt[jj] <= t[i - 1][(int)k] ? 1f : 0f;
                // Use rs to compute next simplex coord.
                sx = Math.Pow(rs[jj], 1f / (double)i);
                // Update sum
                sm += (1f - 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.
            List<int> p = new List<int>((int)size);
            for (int i = 0; i < size; i++)
            {
                p.Add(i);
            }

            // Use p to carry out a matrix 'randperm'
            try
            {
                p.Sort(delegate(int p1, int p2) { return ((p1 == p2) ? 0 : RandomNumber(-1, 1)); });
            }
            catch (ArgumentException e)
            {
                // Sort sometimes throws an IComparable exception.
                // This doesn't seem to affect the result of the algorithm,
                // so catch and ignore the error.
            }

            List<double> ret = new List<double>((int)size);
            for (int i = 0; i < size; i++)
            {
                // Permute & rescale x
                ret.Add((double)(max - min) * (double)x[p[i]] + (double)min);
            }

            return ret;
        }

        /// <summary>
        /// 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.
        /// </summary>
        /// <param name="size">Number of elements in the set.</param>
        /// <param name="sum">Total value that the sum of all the elements must
        /// equal.</param>
        /// <param name="min">Minimum value of a given element.</param>
        /// <param name="max">Maximum value of a given element.</param>
        public List<int> RandFixedSum(int size, int sum, int min, int max)
        {
            List<double> listDbl = RandFixedSum((double)size, (double)sum, (double)min, (double)max);
            List<int> listInt = new List<int>(size);
            int sumOfRounded = 0;

            foreach (double randDbl in listDbl)
            {
                int randInt = (int)Math.Round(randDbl);
                listInt.Add(randInt);
                sumOfRounded += randInt;
            }

            // Randomly permute the order in the columns of listInt and rescale.
            List<int> p = new List<int>((int)size);
            for (int i = 0; i < size; i++)
            {
                p.Add(i);
            }

            // Use p to carry out a matrix 'randperm'
            try
            {
                p.Sort(delegate(int p1, int p2) { return ((p1 == p2) ? 0 : RandomNumber(-1, 1)); });
            }
            catch (ArgumentException e)
            {
                // Sort sometimes throws an IComparable exception.
                // This doesn't seem to affect the result of the algorithm,
                // so catch and ignore the error.
            }

            // The total of the rounded integers is lower than the total
            // that we want.
            if (sumOfRounded < sum)
            {
                // Start by searching for an element with value 'min', and
                // increase if unsuccessful.
                int seekAmt = min;
                while (seekAmt < max && sumOfRounded < sum)
                {
                    // Try and raise the rounded total to equal the required
                    // total, by adding onto a single element in the set.
                    int addAmt = sum - sumOfRounded;
                    bool isAdded = false;

                    // No single element in the set can exceed max. In cases
                    // where adding (sum - sumOfRounded) would violate this,
                    // add as much as is permitted, and leave the rest for
                    // another element.
                    if (seekAmt + addAmt > max)
                    {
                        addAmt = max - seekAmt;
                    }

                    for (int i = 0; i < size; i++)
                    {
                        if (listInt[p[i]] == seekAmt)
                        {
                            listInt[p[i]] += addAmt;
                            sumOfRounded += addAmt;
                            isAdded = true;
                            break;
                        }
                    }

                    if (!isAdded)
                    {
                        seekAmt++;
                    }
                }
            }
            // Or the total of the rounded integers is higher than the total
            // that we want.
            else if (sumOfRounded > sum)
            {
                // Start by searching for an element with value 'max', and
                // decrease if unsuccessful.
                int seekAmt = max;
                while (seekAmt > min && sumOfRounded > sum)
                {
                    // Try and lower the rounded total to equal the required
                    // total, by subtracting from a single element in the set.
                    int subtractAmt = sumOfRounded - sum;
                    bool isSubtracted = false;

                    // No single element in the set can be lower than min.
                    // In cases where subtracting (sumOfRounded - sum)
                    // would violate this, subtract as much as is permitted,
                    // and leave the rest for another element.
                    if (seekAmt - subtractAmt < min)
                    {
                        subtractAmt = seekAmt - min;
                    }

                    for (int i = 0; i < size; i++)
                    {
                        if (listInt[p[i]] == seekAmt)
                        {
                            listInt[p[i]] -= subtractAmt;
                            sumOfRounded -= subtractAmt;
                            isSubtracted = true;
                            break;
                        }
                    }

                    if (!isSubtracted)
                    {
                        seekAmt--;
                    }
                }
            }

            return listInt;
        }
    }

    /// <summary>
    /// Creates a List of doubles, with "size" copies of "value".
    /// Similar to PHP's array_fill() function.
    /// </summary>
    class ListFiller
    {
        public static List<double> GetList(double size, double value)
        {
            List<double> filledList = new List<double>((int)size);

            for (int i = 0; i < size; i++)
            {
                filledList.Add(value);
            }

            return filledList;
        }
    }

    /// <summary>
    /// Utility class for an iterable sequence of values.
    /// </summary>
    class Range : System.Collections.IEnumerable
    {
        private readonly double start;
        private double stop;
        private double step;

        public Range(double start, double stop, double step)
        {
            this.start = start;
            this.stop = stop;
            this.step = step;
        }

        public void Do(Action<double> f)
        {
            foreach (double i in this)
            {
                f(i);
            }
        }

        public System.Collections.IEnumerator GetEnumerator()
        {
            double max = Math.Max(start, stop);
            double min = Math.Min(start, stop);
            if (step == default(double))
            {
                step = (start == min) ? 1 : -1;
            }
            double current = start;
            while (current >= min && current <= max)
            {
                yield return current;
                current += step;
            }
        }

        System.Collections.IEnumerator System.Collections.IEnumerable.GetEnumerator()
        {
            return GetEnumerator();
        }
    }
}