tinguely 0.1.1

use mathru::algebra::linear::Matrix;
use mathru::algebra::linear::Vector;
use mathru::optimization::Gradient;
use crate::model::SupervisedLearn;
use mathru::algebra::abstr::Real;
use rand::Rng;
use mathru::optimization::{Optim};
use mathru::algebra::linear::matrix::Transpose;

/// Logistic regression
///
/// Fore more information: <br>
/// <a href="https://en.wikipedia.org/wiki/Logistic_regression">https://en.wikipedia.org/wiki/Logistic_regression</a>
///
/// # Examples
///
///
pub struct LogisticRegression<T>
{
	beta: Option<Vector<T>>,

	optim: Gradient<T>,
}

impl<T> LogisticRegression<T>
	where T: Real
{
	/// Creates an instance
	///
	/// # Arguments
	/// * 'optim': Optimization algorithm, which is used to optimize the model
	///
	/// # Returns
	/// LogisticRegression instance
	pub fn new(optim: Gradient<T>) -> LogisticRegression<T>
	{
		LogisticRegression
		{
			beta: None,
			optim: optim
		}
	}
	pub fn parameters(self: &Self) -> Option<Vector<T>>
	{
		self.beta.clone()
	}

}

impl<T> SupervisedLearn<Matrix<T>, Vector<T>> for LogisticRegression<T>
	where T: Real
{
    fn predict(self: &Self, x: &Matrix<T>) -> Result<Vector<T>, ()>
	{

		let (m, n) : (usize, usize) = x.dim();
    	let (_input_m, _input_n): (usize, usize) = x.dim();

      	let ones: Matrix<T> = Matrix::<T>::ones(m, n  + 1);

      	let full_input: Matrix<T> = ones.set_slice( x, 0, 1);

		return match self.beta.clone()
		{
			Some(beta) => Ok(LogisticRegressionBase::h_x(&beta, &full_input)),
			None => Err(())
		}
	}

	///
	///
	fn train(self: &mut Self, x: &Matrix<T>, y: &Vector<T>)
	{
		let (x_m, x_n): (usize, usize) = x.dim();
		let (y_m, _y_n): (usize, usize) = x.dim();

		if x_m != y_m
		{
			panic!("Dimension mismatch")
		}

		let ones: Matrix<T> = Matrix::<T>::ones(x_m, x_n + 1);
		let full_input: Matrix<T> = ones.set_slice( x, 0, 1);
		let mut rng = rand::thread_rng();
		let beta_0: Vector<T> = Vector::new_column(x_n + 1, vec![T::from_f64(rng.gen_range(-1.0, 1.0)); x_n + 1]);

		let base: LogisticRegressionBase<T> = LogisticRegressionBase::new(y, &full_input);

		let beta: Vector<T> = self.optim.minimize(&base, &beta_0).arg();
		self.beta = Some(beta);

	}
}


pub struct LogisticRegressionBase<'a, T>
{
	pub x: &'a Matrix<T>,
	pub y: &'a Vector<T>,

}

impl<'a, T> LogisticRegressionBase<'a, T>
	where T: Real
{
	fn new(y: &'a Vector<T>, x: &'a Matrix<T>) -> LogisticRegressionBase<'a, T>
	{
		LogisticRegressionBase
		{
			y: y,
			x: x
		}
	}
}

impl<'a, T> Optim<T> for LogisticRegressionBase<'a, T>
	where T: Real
{
    fn eval(&self, beta: &Vector<T>) -> Vector<T>
	{
		let (m, _n) : (usize, usize) = self.x.dim();

		let mut cost: T = T::zero();

		for i in 0..m
		{
			let x_i: Vector<T> = self.x.get_row(i);
			let y_hat_i: T = self.h_xi(beta, &x_i);

			let y_i: T = *self.y.get(i);
			let cost_i: T = y_i * (y_hat_i.ln()) + (T::one() - y_i) * (T::one() - y_hat_i).ln();
			cost += cost_i;
		}

		return Vector::new_column(1, vec![-cost]);
	}

	fn jacobian(&self, beta: &Vector<T>) -> Matrix<T>
	{
		let y_hat: Vector<T> = LogisticRegressionBase::h_x(beta, self.x);

		let diff = &y_hat - self.y;

		let grad: Vector<T> = &(self.x.clone().transpose()) * &diff;

		return Matrix::from( grad.transpose());
	}

	/// Computes the Hessian at the given value x
	fn hessian(&self, _x: &Vector<T>) -> Matrix<T>
	{
		unimplemented!();
	}

}

impl<'a, T> LogisticRegressionBase<'a, T>
	where T: Real
{
	fn h_xi(self: &'a Self, beta: &'a Vector<T>, x: &'a Vector<T>) -> T
	{
		let (_m, n): (usize, usize) = x.dim();

		let mut y_hat: T = T::zero();

		for k in 0..n
		{
			let x_k: T = x.get(k).clone();
			let beta_k: T = beta.get(k).clone();
			y_hat += beta_k * x_k;
		}

		let v: T = LogisticRegressionBase::sigmoid(&y_hat);

		return v;
	}

	fn sigmoid(z: &T) -> T
	{
		T::one()/(T::one() +  T::exp(-*z))
	}

	fn h_x(beta: &Vector<T>, x: &Matrix<T>) -> Vector<T>
	{
		(x * beta).apply(&LogisticRegressionBase::sigmoid)
	}
}