// Copyright (C) 2014-2015 Guibing Guo
//
// This file is part of LibRec.
//
// LibRec is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// LibRec is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with LibRec. If not, see <http://www.gnu.org/licenses/>.
//

package librec.util;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.HashSet;
import java.util.List;
import java.util.Random;
import java.util.Set;

/**
 * 
 * @author Guo Guibing
 *
 */
public class Randoms {
	private static Random r = new Random(System.currentTimeMillis());

	private static List<Object> _tempList = new ArrayList<>();

	/**
	 * Random generate an integer in [0, range)
	 * 
	 * @param range
	 * @return
	 */
	public static int uniform(int range) {
		return uniform(0, range);
	}

	public static void seed(long seed) {
		r = new Random(seed);
	}

	/**
	 * Random generate an integer in [min, max)
	 * 
	 * @param min
	 * @param max
	 * @return
	 */
	public static int uniform(int min, int max) {
		return min + r.nextInt(max - min);
	}

	/**
	 * return real number uniformly in [0, 1)
	 */
	public static double random() {
		return uniform();
	}

	/**
	 * @return a random number from a given list of numbers
	 */
	public static <T> T random(List<T> data) {
		int idx = uniform(data.size());

		return data.get(idx);
	}

	/**
	 * a random double array with values in [0, 1)
	 * 
	 * @param size
	 *            the size of random array
	 */
	public static double[] doubles(int size) {
		double[] array = new double[size];
		for (int i = 0; i < size; i++)
			array[i] = random();

		return array;
	}

	/**
	 * a random double array with values in [min, max)
	 * 
	 * @param size
	 *            the size of random array
	 */
	public static double[] doubles(double min, double max, int size) {
		double[] array = new double[size];
		for (int i = 0; i < size; i++)
			array[i] = uniform(min, max);

		return array;
	}

	/**
	 * random (uniformly distributed) double in [0, 1)
	 */
	public static double uniform() {
		return uniform(0.0, 1.0);
	}

	/**
	 * random (uniformly distributed) double in [min, max)
	 */
	public static double uniform(double min, double max) {
		return min + (max - min) * r.nextDouble();
	}

	/**
	 * Return a boolean, which is true with probability p, and false otherwise.
	 */
	public static boolean bernoulli(double p) {
		return uniform() < p;
	}

	/**
	 * Return a boolean, which is true with probability .5, and false otherwise.
	 */
	public static boolean bernoulli() {
		return bernoulli(0.5);
	}

	/**
	 * return a real number from a Gaussian distribution with given mean and stddev
	 * 
	 * @param mu
	 *            mean
	 * @param sigma
	 *            stddev
	 * 
	 */
	public static double gaussian(double mu, double sigma) {
		return mu + sigma * r.nextGaussian();
	}

	/**
	 * Randomly sample 1 point from Gamma Distribution with the given parameters. The code is from Mahout
	 * (http://mahout.apache.org/), available under Apache 2 license.
	 * 
	 * @param alpha
	 *            alpha parameter for Gamma Distribution.
	 * @param scale
	 *            scale parameter for Gamma Distribution.
	 * @return a sample point randomly drawn from the given distribution.
	 */
	public static double gamma(double alpha, double scale) {
		double rate = 1 / scale;

		if (alpha <= 0.0 || rate <= 0.0) {
			throw new IllegalArgumentException();
		}

		double gds;
		double b = 0.0;

		// CASE A: Acceptance rejection algorithm gs
		if (alpha < 1.0) {
			b = 1.0 + 0.36788794412 * alpha; // Step 1
			while (true) {
				double p = b * r.nextDouble();
				// Step 2. Case gds <= 1
				if (p <= 1.0) {
					gds = Math.exp(Math.log(p) / alpha);
					if (Math.log(r.nextDouble()) <= -gds) {
						return gds / rate;
					}
				}
				// Step 3. Case gds > 1
				else {
					gds = -Math.log((b - p) / alpha);
					if (Math.log(r.nextDouble()) <= ((alpha - 1.0) * Math.log(gds))) {
						return gds / rate;
					}
				}
			}
		}
		// CASE B: Acceptance complement algorithm gd (gaussian distribution,
		// box muller transformation)
		else {
			double ss = 0.0;
			double s = 0.0;
			double d = 0.0;

			// Step 1. Preparations
			if (alpha != -1.0) {
				ss = alpha - 0.5;
				s = Math.sqrt(ss);
				d = 5.656854249 - 12.0 * s;
			}

			// Step 2. Normal deviate
			double v12;
			double v1;

			do {
				v1 = 2.0 * r.nextDouble() - 1.0;
				double v2 = 2.0 * r.nextDouble() - 1.0;
				v12 = v1 * v1 + v2 * v2;
			} while (v12 > 1.0);

			double t = v1 * Math.sqrt(-2.0 * Math.log(v12) / v12);
			double x = s + 0.5 * t;
			gds = x * x;

			if (t >= 0.0) { // Immediate acceptance
				return gds / rate;
			}

			double u = r.nextDouble();
			if (d * u <= t * t * t) { // Squeeze acceptance
				return gds / rate;
			}

			double q0 = 0.0;
			double si = 0.0;
			double c = 0.0;

			// Step 4. Set-up for hat case
			if (alpha != -1.0) {
				double rr = 1.0 / alpha;
				double q9 = 0.0001710320;
				double q8 = -0.0004701849;
				double q7 = 0.0006053049;
				double q6 = 0.0003340332;
				double q5 = -0.0003349403;
				double q4 = 0.0015746717;
				double q3 = 0.0079849875;
				double q2 = 0.0208333723;
				double q1 = 0.0416666664;

				q0 = ((((((((q9 * rr + q8) * rr + q7) * rr + q6) * rr + q5) * rr + q4) * rr + q3) * rr + q2) * rr + q1)
						* rr;

				if (alpha > 3.686) {
					if (alpha > 13.022) {
						b = 1.77;
						si = 0.75;
						c = 0.1515 / s;
					} else {
						b = 1.654 + 0.0076 * ss;
						si = 1.68 / s + 0.275;
						c = 0.062 / s + 0.024;
					}
				} else {
					b = 0.463 + s - 0.178 * ss;
					si = 1.235;
					c = 0.195 / s - 0.079 + 0.016 * s;
				}
			}

			double v, q;
			double a9 = 0.104089866;
			double a8 = -0.112750886;
			double a7 = 0.110368310;
			double a6 = -0.124385581;
			double a5 = 0.142873973;
			double a4 = -0.166677482;
			double a3 = 0.199999867;
			double a2 = -0.249999949;
			double a1 = 0.333333333;

			// Step 5. Calculation of q
			if (x > 0.0) {
				// Step 6.
				v = t / (s + s);
				if (Math.abs(v) > 0.25) {
					q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
				}
				// Step 7. Quotient acceptance
				else {
					q = q0
							+ 0.5
							* t
							* t
							* ((((((((a9 * v + a8) * v + a7) * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1)
							* v;
				}
				if (Math.log(1.0 - u) <= q) {
					return gds / rate;
				}
			}

			double e7 = 0.000247453;
			double e6 = 0.001353826;
			double e5 = 0.008345522;
			double e4 = 0.041664508;
			double e3 = 0.166666848;
			double e2 = 0.499999994;
			double e1 = 1.000000000;

			// Step 8. Double exponential deviate t
			while (true) {
				double sign_u;
				double e;
				do { // Step 9. Rejection of t
					e = -Math.log(r.nextDouble());
					u = r.nextDouble();
					u = u + u - 1.0;
					sign_u = (u > 0) ? 1.0 : -1.0;
					t = b + (e * si) * sign_u;
				} while (t <= -0.71874483771719);

				// Step 10. New q(t)
				v = t / (s + s);

				if (Math.abs(v) > 0.25) {
					q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
				} else {
					q = q0
							+ 0.5
							* t
							* t
							* ((((((((a9 * v + a8) * v + a7) * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1)
							* v;
				}

				// Step 11.
				if (q <= 0.0) {
					continue;
				}

				// Step 12. Hat acceptance
				double w;
				if (q > 0.5) {
					w = Math.exp(q) - 1.0;
				} else {
					w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) * q + e1) * q;
				}

				if (c * u * sign_u <= w * Math.exp(e - 0.5 * t * t)) {
					x = s + 0.5 * t;
					return x * x / rate;
				}
			}
		}
	}

	/**
	 * Return an integer with a Poisson distribution with mean lambda.
	 */
	public static int poisson(double lambda) {
		// using algorithm given by Knuth
		// see http://en.wikipedia.org/wiki/Poisson_distribution
		int k = 0;
		double p = 1.0;
		double L = Math.exp(-lambda);
		do {
			k++;
			p *= uniform();
		} while (p >= L);
		return k - 1;
	}

	/**
	 * Return a real number with a Pareto distribution with parameter alpha.
	 */
	public static double pareto(double alpha) {
		return Math.pow(1 - uniform(), -1.0 / alpha) - 1.0;
	}

	/**
	 * Return a real number with a Cauchy distribution.
	 */
	public static double cauchy() {
		return Math.tan(Math.PI * (uniform() - 0.5));
	}

	/**
	 * Return a number from a discrete distribution: i with probability a[i]. Precondition: array entries are
	 * nonnegative and their sum (very nearly) equals 1.0.
	 */
	public static int discrete(double[] a) {
		double EPSILON = 1E-6;
		double sum = 0.0;
		for (int i = 0; i < a.length; i++) {
			if (a[i] < 0.0)
				throw new IllegalArgumentException("array entry " + i + " is negative: " + a[i]);
			sum = sum + a[i];
		}
		if (sum > 1.0 + EPSILON || sum < 1.0 - EPSILON)
			throw new IllegalArgumentException("sum of array entries not equal to one: " + sum);

		// the for loop may not return a value when both r is (nearly) 1.0 and when the cumulative sum is less than 1.0 (as a result of floating-point roundoff error)
		while (true) {
			double r = uniform();
			sum = 0.0;
			for (int i = 0; i < a.length; i++) {
				sum = sum + a[i];
				if (sum > r)
					return i;
			}
		}
	}

	/**
	 * Return a real number from an exponential distribution with rate lambda.
	 */
	public static double exp(double lambda) {
		return -Math.log(1 - uniform()) / lambda;
	}

	/**
	 * generate next random integer in a range besides exceptions
	 * 
	 * @param range
	 *            the range located of next integer
	 * @param exceptions
	 *            the exception values when generating integers, sorted first
	 * @return next no-repeated random integer
	 */
	public static int nextInt(int range, int... exceptions) {
		return nextInt(0, range, exceptions);
	}

	/**
	 * generate next random integer in a range [min, max) besides exceptions
	 * 
	 * @param min
	 *            the minimum of range
	 * @param max
	 *            the maximum of range
	 * @param exceptions
	 *            the exception values when generating integers, sorted first
	 * @return next no-repeated random integer
	 */

	public static int nextInt(int min, int max, int... exceptions) {
		int next;
		while (true) {
			next = min + r.nextInt(max - min);
			if (exceptions != null && exceptions.length > 0 && Arrays.binarySearch(exceptions, next) >= 0) {
				continue;
			}
			if (_tempList.contains(next))
				continue;
			else {
				_tempList.add(next);
				break;
			}
		}
		return next;
	}

	/**
	 * Generate no repeat {@code size} indexes from {@code min} to {@code max}
	 * 
	 */
	public static int[] indexs(int size, int min, int max) {
		Set<Integer> used = new HashSet<>();
		int[] index = new int[size];

		for (int i = 0; i < index.length; i++) {
			while (true) {
				int ind = uniform(min, max);

				if (!used.contains(ind)) {
					index[i] = ind;
					used.add(ind);
					break;
				}
			}
		}

		return index;
	}

	public static void clearCache() {
		_tempList.clear();
	}

	/**
	 * generate next integers array with no repeated elements
	 * 
	 * @param length
	 *            the length of the array
	 * @param range
	 *            the index range of the array, default [0, range)
	 * @param exceptions
	 *            the exceptions when generating values
	 * @return ascending sorted integer array
	 * @throws Exception
	 *             if the range is less than length, an exception will be thrown
	 */
	public static int[] nextIntArray(int length, int range) throws Exception {
		return nextIntArray(length, 0, range, null);
	}

	public static int[] nextIntArray(int length, int range, int... exceptions) throws Exception {
		return nextIntArray(length, 0, range, exceptions);
	}

	public static int[] nextIntArray(int length, int min, int max) throws Exception {
		return nextIntArray(length, min, max, null);
	}

	public static int[] nextIntArray(int length, int min, int max, int... exceptions) throws Exception {
		int maxLen = max - min; // because max itself is not counted
		if (maxLen < length)
			throw new Exception("The range is less than legth");

		int[] index = new int[length];
		if (maxLen == length) {
			for (int i = 0; i < length; i++)
				index[i] = min + i;
		} else {
			Randoms.clearCache();
			for (int i = 0; i < index.length; i++)
				index[i] = Randoms.nextInt(min, max, exceptions);
			Arrays.sort(index);
		}
		return index;
	}

	/**
	 * Generate a set of random (unique) integers in the range [min, max) with length {@code length}
	 * 
	 * @return a set of unique integers
	 */
	public static List<Integer> randInts(int length, int min, int max) throws Exception {
		int len = max - min;
		if (len < length)
			throw new Exception("The range is less than legth");

		Set<Integer> ints = new HashSet<>();

		while (true) {
			int rand = min + r.nextInt(max - min);
			ints.add(rand);

			if (ints.size() >= length)
				break;
		}

		List<Integer> res = new ArrayList<>(ints);
		Collections.sort(res);

		return res;
	}

	/**
	 * Get a normalize array of probabilities
	 * 
	 * @param size
	 *            array size
	 */
	public static double[] randProbs(int size) {
		if (size < 1)
			throw new IllegalArgumentException("The size param must be greate than zero");

		double[] pros = new double[size];

		int sum = 0;
		for (int i = 0; i < pros.length; i++) {
			//avoid zero
			pros[i] = r.nextInt(size) + 1;
			sum += pros[i];
		}

		//normalize
		for (int i = 0; i < pros.length; i++) {
			pros[i] = pros[i] / sum;
		}

		return pros;
	}

	public static int[] ints(int range, int size) {
		int[] data = new int[size];

		for (int i = 0; i < size; i++)
			data[i] = uniform(range);

		return data;
	}

	public static int[] ints(int min, int max, int size) {
		int[] data = new int[size];

		for (int i = 0; i < size; i++)
			data[i] = uniform(min, max);

		return data;
	}

	public static List<Double> list(int size) {
		return list(size, 0, 1, false);
	}

	public static List<Double> list(int size, int min, int max) {
		return list(size, min, max, false);
	}

	public static List<Double> list(int size, int min, int max, boolean isInteger) {
		List<Double> list = new ArrayList<>(size);
		for (int i = 0; i < size; i++) {
			if (isInteger)
				list.add(uniform(min, max) + 0.0);
			else
				list.add(uniform(min + 0.0, max + 0.0));
		}

		return list;
	}

	/**
	 * Generate a permutation from min to max
	 * 
	 * @param min
	 *            the minimum value
	 * @param max
	 *            the maximum value
	 * @return a permutation
	 */
	public static List<Integer> permute(int min, int max) {
		List<Integer> list = new ArrayList<>();

		int len = max - min + 1;
		for (int i = 0; i < len; i++) {
			while (true) {
				int index = uniform(min, max + 1);

				if (!list.contains(index)) {
					list.add(index);
					break;
				}
			}
		}

		return list;
	}
}
