/*
	This file acompanies the article 
		"Generating Random Numbers with Arbitrary Distributions" 
	that can be found here:
		http://code-spot.co.za/2008/09/21/generating-random-numbers-with-arbitrary-distributions/
	
	@author
		Herman Tulleken (herman.tulleken@gmail.com)

	This program illustrates how to generate random numbers that 
	follow a distribution described by a piece-wise linear curve.
*/

#include "stdio.h"
#include "time.h"
#include "stdlib.h"

/**
	The header files below are part of the Special Numbers Library. 
	
	Download from
		http://code.google.com/p/specialnumbers/
*/
#include "ResponseCurve.h"
#include "XYResponseCurve.h"

using namespace luma::numbers;

int main()
{
	srand((unsigned int) time(0));

	const unsigned int oldSampleCount = 9;
	const unsigned int newSampleCount = 1000;

	// Some arbitrary distribution curve
	float inputSamples[] = {-20, -10, 0, 10, 20, 30, 40, 50, 60};
	float outputSamples[] = {2, 10, 80, 75, 60, 30, 10, 5, 1};

	float accumulativeOutputSamples[oldSampleCount] = {0};

	// Step 1. Calculate accumulative output

	accumulativeOutputSamples[0] = outputSamples[0];

	for(int i = 1; i < oldSampleCount; i++)
	{
		accumulativeOutputSamples[i] = accumulativeOutputSamples[i - 1] + outputSamples[i];
	}

	float newInputMin = 2.0f;
	float newInputMax = 273.0f;

	float newOutputMax = 60.0f;
	float newOutputMin = -20.0f;

	// Step2. Load inverse into XY response curve
	XYResponseCurve<float, oldSampleCount> xyCurve(accumulativeOutputSamples, inputSamples);

	float newOutputSamples[newSampleCount] = {0};

	// Step 3. Gather samples for ordinary reponse curve
	for(unsigned int i = 0; i < newSampleCount; i++)
	{
		float input = (i / ((float) newSampleCount - 1)) * (newInputMax - newInputMin) + newInputMin;
		newOutputSamples[i] = xyCurve(input);

		// Used for debugging.
		//printf("%f %f\n", input, newOutputSamples[i]);
	}

	// Construct ordinary response curve from samples.
	ResponseCurve<float, newSampleCount> curve(newInputMin, newInputMax, newOutputSamples);

	for(unsigned int i = 0; i < newSampleCount; i++)
	{
		float input = (i / ((float) newSampleCount - 1)) * (newInputMax - newInputMin) + newInputMin;
		float output = curve(input);

		// Used for debugging.
		//printf("%f %f\n", input, output);
	}

	//test the distribution

	const int testSampleCount = 100;
	int count[testSampleCount] = {0};
	
	// generate 10 000 random numbers, and check count occurrences in distribution bands.
	for(int i = 0; i < 10000; i++)
	{
		int uniformRandInt = rand();
		float uniformRandVal =  uniformRandInt / (float) RAND_MAX; //random value between 0 and 1 

		//Step 4. Map random value to the appropriate input value for the response curve
		float newInput = uniformRandVal*(newInputMax - newInputMin) + newInputMin;

		//The random value that follows the desired distribution.
		float curvedRandVal = curve(newInput);	

		// Calculate the distribution band of the random value.
		int countIndex = (int)((curvedRandVal - newOutputMin) / (newOutputMax - newOutputMin) * testSampleCount);

		// Used for debugging.
		//printf("%d %f %f %f %d\n", rr, r, newInput, dr, countIndex);

		if (countIndex < testSampleCount)
		{
			count[countIndex]++;
		}
		else // a rare occurence, but possible. Just bundle with the last distribution band.
		{
			count[testSampleCount - 1]++;
		}
	}
	
	// Print out, so that the values can be pasted into Excel or Calc.
	for(int i = 0; i < testSampleCount; i++)
	{
		printf("%d %d\n", i, count[i]);
	}		
	
}