#include "Gaussian2d.h"

Gaussian2d::Gaussian2d()
{
	m_mux = 0.0;
	m_muy = 0.0;
	m_Cxx = 1.0;
	m_Cyy = 1.0;
	m_Cxy = 1.0;

	m_factor = 0.0;
	m_beta_Cxy = 0.0;
	m_inv_sqrt_Cxx = 1.0;
}

Gaussian2d::Gaussian2d(double mux, double muy, double Cxx, double Cyy, double Cxy)
{
	m_mux = mux;
	m_muy = muy;
	m_Cxx = Cxx;
	m_Cyy = Cyy;
	m_Cxy = Cxy;

	m_factor = 0.0;
	m_beta_Cxy = 0.0;
	m_inv_sqrt_Cxx = 1.0 / std::sqrt(Cxx);
}

std::unique_ptr<Gaussian2d> Gaussian2d::construct(double mux, double muy, double Cxx, double Cyy, double Cxy)
{
	return std::unique_ptr<Gaussian2d>(new Gaussian2d(mux, muy, Cxx, Cyy, Cxy));
}

void Gaussian2d::fix_beta(double beta)
{
	m_factor = std::exp(beta * (m_muy + 0.5 * beta * m_Cyy));
	m_beta_Cxy = beta * m_Cxy;
}


double Gaussian2d::e_beta_expectation(double beta, double K, double mux, double muy, double Cxx, double Cyy, double Cxy)
{
	return std::exp(beta * (muy + 0.5 * beta * Cyy)) * Algorithm::normal_CDF((mux + beta * Cxy - K) / std::sqrt(Cxx));
}

double Gaussian2d::e_beta_expectation(double beta, double K) const
{
	return std::exp(beta * (m_muy + 0.5 * beta * m_Cyy)) * Algorithm::normal_CDF((m_mux + beta * m_Cxy - K) * m_inv_sqrt_Cxx);
}

double Gaussian2d::e_beta_expectation(double K) const
{
	return m_factor * Algorithm::normal_CDF((m_mux + m_beta_Cxy - K) * m_inv_sqrt_Cxx);
}