GithubNeuralNetwork Beta

using System;
using System.Linq;

public class GithubNeuralNetwork
{
    private int[] _layers;
    private Matrix[] _weights;
    private Matrix[] _biases;
    private Func<Matrix, Matrix>[] _activationFunctions;
    private double _learningRate;
    private double _epsilon;
    private Matrix[] _gamma;
    private Matrix[] _beta;
    private double _initialLearningRate;
    private double _decayRate;
    private string _optimizer;
    private Matrix[] _movingMeans;
    private Matrix[] _movingVariances;
    private Matrix[] _mWeights;
    private Matrix[] _vWeights;
    private Matrix[] _mBiases;
    private Matrix[] _vBiases;
    private Matrix[] _mGamma;
    private Matrix[] _vGamma;
    private Matrix[] _mBeta;
    private Matrix[] _vBeta;
    private int _t;

    public GithubNeuralNetwork(double learningRate, double epsilon, string optimizer, double decayRate, params int[] layers)
    {
        _layers = layers;
        _weights = new Matrix[layers.Length - 1];
        _biases = new Matrix[layers.Length - 1];
        _activationFunctions = new Func<Matrix, Matrix>[layers.Length - 1];
        _learningRate = learningRate;
        _epsilon = epsilon;
        _gamma = new Matrix[layers.Length - 1];
        _beta = new Matrix[layers.Length - 1];
        _initialLearningRate = learningRate;
        _decayRate = decayRate;
        _optimizer = optimizer;
        _movingMeans = new Matrix[layers.Length - 1];
        _movingVariances = new Matrix[layers.Length - 1];
        _mWeights = new Matrix[layers.Length - 1];
        _vWeights = new Matrix[layers.Length - 1];
        _mBiases = new Matrix[layers.Length - 1];
        _vBiases = new Matrix[layers.Length - 1];
        _mGamma = new Matrix[layers.Length - 1];
        _vGamma = new Matrix[layers.Length - 1];
        _mBeta = new Matrix[layers.Length - 1];
        _vBeta = new Matrix[layers.Length - 1];
        _t = 1;

        InitializeWeightsAndBiases();
        SetActivationFunctions();
    }

    private void InitializeWeightsAndBiases()
    {
        Random rand = new Random();
        for (int i = 0; i < _weights.Length; i++)
        {
            _weights[i] = XavierInitialization(_layers[i + 1], _layers[i], rand);
            _biases[i] = Matrix.Zeros(_layers[i + 1], 1);
            _gamma[i] = Matrix.Ones(_layers[i + 1], 1);
            _beta[i] = Matrix.Zeros(_layers[i + 1], 1);

            _movingMeans[i] = Matrix.Zeros(_layers[i + 1], 1);
            _movingVariances[i] = Matrix.Ones(_layers[i + 1], 1);

            _mWeights[i] = Matrix.Zeros(_weights[i].RowCount, _weights[i].ColumnCount);
            _vWeights[i] = Matrix.Zeros(_weights[i].RowCount, _weights[i].ColumnCount);
            _mBiases[i] = Matrix.Zeros(_biases[i].RowCount, _biases[i].ColumnCount);
            _vBiases[i] = Matrix.Zeros(_biases[i].RowCount, _biases[i].ColumnCount);
            _mGamma[i] = Matrix.Zeros(_gamma[i].RowCount, _gamma[i].ColumnCount);
            _vGamma[i] = Matrix.Zeros(_gamma[i].RowCount, _gamma[i].ColumnCount);
            _mBeta[i] = Matrix.Zeros(_beta[i].RowCount, _beta[i].ColumnCount);
            _vBeta[i] = Matrix.Zeros(_beta[i].RowCount, _beta[i].ColumnCount);
        }
    }

    private void SetActivationFunctions()
    {
        for (int i = 0; i < _activationFunctions.Length; i++)
        {
            _activationFunctions[i] = MatrixFunctions.ReLU; // Using Rectified Linear Unit (ReLU) as default activation function
        }
    }

    public void SetActivationFunction(int layerIndex, Func<Matrix, Matrix> activationFunction)
    {
        if (layerIndex < 0 || layerIndex >= _activationFunctions.Length)
        {
            throw new ArgumentOutOfRangeException(nameof(layerIndex));
        }
        _activationFunctions[layerIndex] = activationFunction;
    }

    private Matrix XavierInitialization(int rows, int cols, Random rand)
    {
        double scale = Math.Sqrt(2.0 / (rows + cols));
        return Matrix.RandomMatrix(rows, cols, rand) * scale;
    }

    private Matrix LayerNormalization(Matrix x, Matrix gamma, Matrix beta, int layerIndex)
    {
        Matrix mean = MatrixFunctions.Mean(x, axis: 1);
        Matrix variance = MatrixFunctions.Variance(x, axis: 1);

        _movingMeans[layerIndex] = (_movingMeans[layerIndex] * 0.9) + (mean * 0.1);
        _movingVariances[layerIndex] = (_movingVariances[layerIndex] * 0.9) + (variance * 0.1);

        Matrix normalized = (x - mean) / MatrixFunctions.Sqrt(variance + _epsilon);
        return (gamma * normalized) + beta;
    }

    private void Optimizer(Matrix[] gradientsWeights, Matrix[] gradientsBiases, Matrix[] gradientsGamma, Matrix[] gradientsBeta)
    {
        double beta1 = 0.9;
        double beta2 = 0.999;
        double epsilon = 1e-8;

        for (int i = 0; i < _weights.Length; i++)
        {
            if (_optimizer == "Adam")
            {
                _mWeights[i] = (beta1 * _mWeights[i]) + ((1 - beta1) * gradientsWeights[i]);
                _vWeights[i] = (beta2 * _vWeights[i]) + ((1 - beta2) * (gradientsWeights[i] * gradientsWeights[i]));

                _mBiases[i] = (beta1 * _mBiases[i]) + ((1 - beta1) * gradientsBiases[i]);
                _vBiases[i] = (beta2 * _vBiases[i]) + ((1 - beta2) * (gradientsBiases[i] * gradientsBiases[i]));

                _mGamma[i] = (beta1 * _mGamma[i]) + ((1 - beta1) * gradientsGamma[i]);
                _vGamma[i] = (beta2 * _vGamma[i]) + ((1 - beta2) * (gradientsGamma[i] * gradientsGamma[i]));

                _mBeta[i] = (beta1 * _mBeta[i]) + ((1 - beta1) * gradientsBeta[i]);
                _vBeta[i] = (beta2 * _vBeta[i]) + ((1 - beta2) * (gradientsBeta[i] * gradientsBeta[i]));

                Matrix mHatWeights = _mWeights[i] / (1 - Math.Pow(beta1, _t));
                Matrix vHatWeights = _vWeights[i] / (1 - Math.Pow(beta2, _t));

                Matrix mHatBiases = _mBiases[i] / (1 - Math.Pow(beta1, _t));
                Matrix vHatBiases = _vBiases[i] / (1 - Math.Pow(beta2, _t));

                Matrix mHatGamma = _mGamma[i] / (1 - Math.Pow(beta1, _t));
                Matrix vHatGamma = _vGamma[i] / (1 - Math.Pow(beta2, _t));

                Matrix mHatBeta = _mBeta[i] / (1 - Math.Pow(beta1, _t));
                Matrix vHatBeta = _vBeta[i] / (1 - Math.Pow(beta2, _t));

                _weights[i] -= (_learningRate * mHatWeights) / (MatrixFunctions.Sqrt(vHatWeights) + epsilon);
                _biases[i] -= (_learningRate * mHatBiases) / (MatrixFunctions.Sqrt(vHatBiases) + epsilon);
                _gamma[i] -= (_learningRate * mHatGamma) / (MatrixFunctions.Sqrt(vHatGamma) + epsilon);
                _beta[i] -= (_learningRate * mHatBeta) / (MatrixFunctions.Sqrt(vHatBeta) + epsilon);

                _t++;
            }
            else if (_optimizer == "AdaGrad")
            {
                _vWeights[i] += gradientsWeights[i] * gradientsWeights[i];
                _vBiases[i] += gradientsBiases[i] * gradientsBiases[i];
                _vGamma[i] += gradientsGamma[i] * gradientsGamma[i];
                _vBeta[i] += gradientsBeta[i] * gradientsBeta[i];

                _weights[i] -= (_learningRate / (MatrixFunctions.Sqrt(_vWeights[i]) + epsilon)) * gradientsWeights[i];
                _biases[i] -= (_learningRate / (MatrixFunctions.Sqrt(_vBiases[i]) + epsilon)) * gradientsBiases[i];
                _gamma[i] -= (_learningRate / (MatrixFunctions.Sqrt(_vGamma[i]) + epsilon)) * gradientsGamma[i];
                _beta[i] -= (_learningRate / (MatrixFunctions.Sqrt(_vBeta[i]) + epsilon)) * gradientsBeta[i];
            }
            else
            {
                _weights[i] -= _learningRate * gradientsWeights[i];
                _biases[i] -= _learningRate * gradientsBiases[i];
                _gamma[i] -= _learningRate * gradientsGamma[i];
                _beta[i] -= _learningRate * gradientsBeta[i];
            }
        }
    }

    private Matrix FeedForward(Matrix input, bool training)
    {
        Matrix outputs = input;
        Matrix[] dropoutMasks = new Matrix[_weights.Length];

        for (int i = 0; i < _weights.Length; i++)
        {
            if (training)
            {
                outputs = outputs * _weights[i] + _biases[i];
                outputs = LayerNormalization(outputs, _gamma[i], _beta[i], i);
                outputs = outputs.Map(_activationFunctions[i]);
            }
            else
            {
                outputs = outputs * _weights[i] + _biases[i];
                outputs = outputs.Map(_activationFunctions[i]);
            }
        }
        return outputs;
    }

    private void Backpropagation(Matrix input, Matrix target)
    {
        Matrix[] outputs = new Matrix[_weights.Length + 1];
        outputs[0] = input;

        for (int i = 0; i < _weights.Length; i++)
        {
            outputs[i + 1] = outputs[i] * _weights[i] + _biases[i];
            outputs[i + 1] = outputs[i + 1].Map(_activationFunctions[i]);
        }

        Matrix[] errors = new Matrix[_weights.Length];
        errors[_weights.Length - 1] = outputs[^1] - target;

        for (int i = _weights.Length - 2; i >= 0; i--)
        {
            errors[i] = (_weights[i + 1].Transpose() * errors[i + 1]).MapDerivative(_activationFunctions[i]);
        }

        Matrix[] gradientsWeights = new Matrix[_weights.Length];
        Matrix[] gradientsBiases = new Matrix[_weights.Length];
        Matrix[] gradientsGamma = new Matrix[_weights.Length];
        Matrix[] gradientsBeta = new Matrix[_weights.Length];

        for (int i = 0; i < _weights.Length; i++)
        {
            gradientsWeights[i] = errors[i] * outputs[i].Transpose();
            gradientsBiases[i] = errors[i];
            gradientsGamma[i] = errors[i] * _movingMeans[i];
            gradientsBeta[i] = errors[i] * _movingVariances[i];
        }

        Optimizer(gradientsWeights, gradientsBiases, gradientsGamma, gradientsBeta);
    }

    private void LearningRateScheduler(int epoch)
    {
        _learningRate = _initialLearningRate / (1 + _decayRate * epoch);
    }

    public void Train(Matrix[] inputs, Matrix[] targets, int epochs, int batchSize)
    {
        Random rand = new Random();

        for (int epoch = 0; epoch < epochs; epoch++)
        {
            for (int i = 0; i < inputs.Length; i += batchSize)
            {
                Matrix[] batchInputs = inputs.Skip(i).Take(batchSize).ToArray();
                Matrix[] batchTargets = targets.Skip(i).Take(batchSize).ToArray();

                for (int j = 0; j < batchSize; j++)
                {
                    Matrix outputs = FeedForward(batchInputs[j], true);
                    Backpropagation(batchInputs[j], batchTargets[j]);
                }
            }

            LearningRateScheduler(epoch);
        }
    }

    public Matrix Predict(Matrix input)
    {
        return FeedForward(input, false);
    }
}