/*
 * Normaldistribution.cpp
 *
 *  Created on: 20.05.2011
 *      Author: sven
 */

#include "Normaldistribution.h"
#include <cmath>
#include <cstdlib>
#include <iostream>

Normaldistribution::Normaldistribution() {
}

Normaldistribution::~Normaldistribution() {
}

void Normaldistribution::calcRandomVector(const CholeskyMatrix& D, const std::vector<double>& mu , std::vector<double>& y)
{
	std::vector<double> u(nenEndOfCompanies);
	calcStandardNormalDistributedRandomVector(nenEndOfCompanies,u);
	for (int i=0;i<nenEndOfCompanies;i++)
	{
		double sum = 0.0;
		for (int j=0; j< nenEndOfCompanies ; j++)
		{
			sum += (D.at(i).at(j) * u.at(j));
		}
		y.at(i) = sum + mu.at(i);
//		std::cout << y.at(i) << std::endl;
	}
//	std::cout << std::endl;
}

void Normaldistribution::calcStandardNormalDistributedRandomVector(const int size, std::vector<double>& u)
{
	double x1,x2;
	double v1,v2;
	double u1,u2;
	double s;
	int counter = 0;
	while (counter < size)
	{
		do
		{
			x1 = ((double) rand()) / (RAND_MAX + 1.0);
			x2 = ((double) rand()) / (RAND_MAX + 1.0);
			v1 = (2.0 * x1) - 1.0;
			v2 = (2.0 * x2) - 1.0;
			s = (v1*v1) + (v2*v2);
		} while (s >= 1.0);
		u1 = v1 * std::sqrt( - (2.0 / s) * std::log(s) );
		u2 = v2 * std::sqrt( - (2.0 / s) * std::log(s) );
		if (counter < size) {
			u.at(counter) = u1;
			counter++;
		}
		if (counter < size) {
			u.at(counter) = u2;
			counter++;
		}
	}
}