#pragma once

#include "core/omp_config.h"

#include "utils/vector.h"
#include "utils/matrix.h"
#include "modules/neural_networks/layers/Layer.h"

#include "numerics/max.h"
#include "numerics/sub.h"
#include "numerics/exp.h"
#include "numerics/div.h"



namespace neural_networks{

	template <typename T>
	struct Activation_Softmax : Layer<T>{

			//utils::Matrix<T> exp_values;
			//utils::Matrix<T> probabilities;
			utils::Matrix<T> outputs;

			utils::Matrix<T> dinputs;

	        void forward(const utils::Matrix<T>& inputs){
	
	        	// Get unnormalized probabilities
			    utils::Matrix<T> exp_values = numerics::exp(numerics::sub_colwise(inputs, numerics::max_rowwise(inputs)));

			    // Normalize them for each sample
			    utils::Matrix<T> probabilities = numerics::div_colwise(exp_values, numerics::sum_rowwise(exp_values));
			    outputs = probabilities;
	        }

	        void backward(const utils::Matrix<T>& dvalues){

	        	const uint64_t rows = dvalues.rows();
	        	const uint64_t cols = dvalues.cols();

	        	if ((dinputs.rows() != rows) || dinputs.cols() != cols){
	        		dinputs.resize(rows, cols);
	        	}

	        	for (uint64_t i = 0; i < rows; ++i){
	        		T dot = T{0};
	        		for (uint64_t j = 0; j < cols; ++j){
	        			dot += outputs(i,j) * dvalues(i,j);
	        		}
	        		for (uint64_t j = 0; j < cols; ++j){
	        			dinputs(i,j) = outputs(i,j) * (dvalues(i,j) - dot);
	        		}
	        	}
	        }
	};

} // end namespace neural_networks