// SFunk.cs : Defines the entry point for the console application.
/// (C) 2017, Andrew Polar under GPL ver. 3.
//   LICENSE
//
//   This program is free software; you can redistribute it and/or
//   modify it under the terms of the GNU General Public License as
//   published by the Free Software Foundation; either version 3 of
//   the License, or (at your option) any later version.
//
//   This program is distributed in the hope that it will be useful, but
//   WITHOUT ANY WARRANTY; without even the implied warranty of
//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//   General Public License for more details at
//   Visit <http://www.gnu.org/copyleft/gpl.html>.
//
// Last modification 01.25.2017
//
// The concept is taken from http://sifter.org/~simon/journal/20061211.html
//

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace EffectiveComputations
{
    class Column
    {
        public List<double> Elements;
    }

    class Row
    {
        public List<Column> Columns;
    }

    class Matrix
    {
        public List<Row> Rows;
    }

    class SFunk
    {
        private int m_rows;
        private int m_cols;
        private int m_dimension;
        private double[,] model;
        private double[,] u;
        private double[,] m;

        const double m_learningRate = 0.001;  //if you see negative values in cofactors, change it to smaller value
        const double m_regularization = 0.005;

        public void GetDecomposition(Matrix datamatrix, int dimension)
        {
            m_rows = datamatrix.Rows.Count;
            m_cols = datamatrix.Rows[0].Columns.Count;
            model = new double[m_rows, m_cols];
            m_dimension = dimension;
            u = new double[m_rows, m_dimension];
            m = new double[m_dimension, m_cols];
            IntitalizeAverage(datamatrix);
            RunDecomposition(datamatrix);
        }

        private void IntitalizeAverage(Matrix datamatrix)
        {
            double average = 0.0;
            int counter = 0;
            foreach (Row row in datamatrix.Rows)
            {
                foreach (Column column in row.Columns)
                {
                    foreach (double element in column.Elements)
                    {
                        average += element;
                        ++counter;
                    }
                }
            }
            average /= (double)(counter);
            average /= (double)(m_dimension);
            average = Math.Sqrt(average);

            for (int i = 0; i < m_rows; ++i)
            {
                for (int j = 0; j < m_dimension; ++j)
                {
                    u[i, j] = average;
                }
            }

            for (int i = 0; i < m_dimension; ++i)
            {
                for (int j = 0; j < m_cols; ++j)
                {
                    m[i, j] = average;
                }
            }
        }

        private double ComputeRate(int row, int col)
        {
            double rate = 0.0;
            for (int i = 0; i < m_dimension; ++i)
            {
                rate += u[row, i] * m[i, col];
            }
            return rate;
        }

        private double GetRMSE(Matrix datamatrix)
        {
            double RMSE = 0.0;
            int i = 0;
            foreach (Row row in datamatrix.Rows)
            {
                int j = 0;
                foreach (Column column in row.Columns)
                {
                    double value = 0.0;
                    for (int k = 0; k < m_dimension; ++k)
                    {
                        value += u[i, k] * m[k, j];
                    }

                    foreach (double element in column.Elements)
                    {
                        RMSE += (element - value) * (element - value);
                    }
                    ++j;
                }
                ++i;
            }

            RMSE /= m_rows;
            RMSE /= m_cols;
            return Math.Sqrt(RMSE);
        }

        public void ShowLeft()
        {
            Console.WriteLine("Left matrix");
            for (int i = 0; i < m_rows; ++i)
            {
                for (int j = 0; j < m_dimension; ++j)
                {
                    Console.Write("{0:0.000} ", u[i, j]);
                }
                Console.WriteLine();
            }
            Console.WriteLine();
        }

        public void ShowRight()
        {
            Console.WriteLine("Right matrix");
            for (int i = 0; i < m_dimension; ++i)
            {
                for (int j = 0; j < m_cols; ++j)
                {
                    Console.Write("{0:0.000} ", m[i, j]);
                }
                Console.WriteLine();
            }
            Console.WriteLine();
        }

        public void ShowModel()
        {
            Console.WriteLine("Left * Right product");
            for (int i = 0; i < m_rows; ++i)
            {
                for (int j = 0; j < m_cols; ++j)
                {
                    model[i, j] = 0.0;
                    for (int k = 0; k < m_dimension; ++k)
                    {
                        model[i, j] += u[i, k] * m[k, j];
                    }
                }
            }

            for (int i = 0; i < m_rows; ++i)
            {
                for (int j = 0; j < m_cols; ++j)
                {
                    Console.Write("{0:0.000} ", model[i, j]);

                }
                Console.WriteLine();
            }
            Console.WriteLine();
        }

        public void ShowModelRound()
        {
            for (int i = 0; i < m_rows; ++i)
            {
                for (int j = 0; j < m_cols; ++j)
                {
                    Console.Write("{0:0.0} ", Math.Round(model[i, j]));

                }
                Console.WriteLine();
            }
            Console.WriteLine();
        }
        public void ShowOriginalData(Matrix datamatrix)
        {
            Console.WriteLine("Original data represents 3d matrix\n");
            int nLayers = 0;
            foreach (Row row in datamatrix.Rows)
            {
                foreach (Column column in row.Columns)
                {
                    if (nLayers < column.Elements.Count) nLayers = column.Elements.Count;
                }
            }

            for (int layer = 0; layer < nLayers; ++layer)
            {
                Console.WriteLine("Layer {0}", layer);
                foreach (Row row in datamatrix.Rows)
                {
                    foreach (Column column in row.Columns)
                    {
                        bool isavailable = false;
                        int cnt = 0;
                        foreach (double element in column.Elements)
                        {
                            if (cnt == layer)
                            {
                                Console.Write("{0:0.0} ", element);
                                isavailable = true;
                            }
                            ++cnt;
                        }
                        if (false == isavailable)
                        {
                            Console.Write("N/A  ");
                        }
                    }
                    Console.WriteLine();
                }
                Console.WriteLine();
            }
            Console.WriteLine();
        }

        public double[,] GetLeft()
        {
            return u;
        }

        public double[,] GetRight()
        {
            return m;
        }

        private void RunDecomposition(Matrix datamatrix)
        {
            double RMSE = 0.0;
            double prevRMSE = 0.0;
            int epoch = 0;
            while (true)
            {
                for (int dim = 0; dim < m_dimension; ++dim)
                {
                    int i = 0;
                    foreach (Row row in datamatrix.Rows)
                    {
                        int j = 0;
                        foreach (Column column in row.Columns)
                        {
                            foreach (double element in column.Elements)
                            {
                                double err = (element - ComputeRate(i, j));
                                double uOld = u[i, dim];
                                u[i, dim] += m_learningRate * (err * m[dim, j] - m_regularization * u[i, dim]);
                                m[dim, j] += m_learningRate * (err * uOld - m_regularization * m[dim, j]);
                            }
                            ++j;
                        }
                        ++i;
                    }
                }
                if (++epoch >= 8000)
                {
                    break;
                }
                RMSE = GetRMSE(datamatrix);
                if (Math.Abs(prevRMSE - RMSE) / (RMSE) < 1e-7)
                {
                    break;
                }
                else
                {
                    prevRMSE = RMSE;
                }
            }
            Console.WriteLine("Number of epochs {0}, RMSE {1}\n", epoch, RMSE);
        }
    }

    class Program
    {
        static void Main(string[] args)
        {
            //////// Data assignment for self-test ////////////////
            double[,] data = 
            {
                {30.00, 45.00, 40.00, 11.00},
                {30.00, 45.00, 33.00, 11.00},
                {30.00, 45.00, 64.00, 20.00},
                {24.00, 36.00, 71.00, 25.00},
                {36.00, 54.00, 79.00, 32.00}
            };

            int rows = data.GetLength(0);
            int cols = data.GetLength(1);

            Random rnd = new Random();

            Matrix matrix = new Matrix();
            matrix.Rows = new List<Row>();
            for (int i = 0; i < rows; ++i)
            {
                matrix.Rows.Add(new Row());
                matrix.Rows[i].Columns = new List<Column>();
                for (int j = 0; j < cols; ++j)
                {
                    matrix.Rows[i].Columns.Add(new Column());
                    matrix.Rows[i].Columns[j].Elements = new List<double>();
                    //we make 3d matrix
                    matrix.Rows[i].Columns[j].Elements.Add(data[i, j]);
                    matrix.Rows[i].Columns[j].Elements.Add(data[i, j] + rnd.Next() % 4);
                    matrix.Rows[i].Columns[j].Elements.Add(data[i, j] + rnd.Next() % 4);
                }
            }

            //here we remove two element arrays
            matrix.Rows[0].Columns[1].Elements.Clear();
            matrix.Rows[3].Columns[2].Elements.Clear();

            ///////// Start decomposition ////////////////

            SFunk sf = new SFunk();
            sf.GetDecomposition(matrix, 3);
            sf.ShowLeft();
            sf.ShowRight();
            sf.ShowModel();
            sf.ShowOriginalData(matrix);
        }
    }
}