using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using System.Threading;
using Azimuth.DenseLinearAlgebra;

namespace Azimuth.PerformanceTests
{
    public partial class MainForm : Form
    {
        bool QRDecompostion { get; set; }
        int progress;
        Random rand = new Random();
        
        public MainForm()
        {
            InitializeComponent();
        }

        private delegate void ProgressIncrement();

        private void StartTest_Click(object sender, EventArgs e)
        {
            Thread testThread = null;
            int workToBeDone = 0;
            progress = 0;

            if (QRDecompositionRadio.Checked)
            {
                testThread = RunQRTest(out workToBeDone);
            }
            else if (LinearLeastSquaresRadio.Checked)
            {
                testThread = RunLLSTest(out workToBeDone);
            }
            else if (BisectionRadio.Checked)
            {
                testThread = RunBisectionTest(out workToBeDone);
            }
            else if (radioButton3.Checked)
            {
                testThread = RunLLSInequalitiesTest(out workToBeDone);
            }
            else if (matMul.Checked)
            {
                testThread = RunMatMul(out workToBeDone);
            }
            else if (sparseMMCheckbox.Checked)
            {
                testThread = RunSparseMultiply(out workToBeDone);
            }
            else if (choleskyCheck.Checked)
            {
                testThread = RunCholesky(out workToBeDone);
            }
            else if (DantzigButton.Checked)
            {
                testThread = DantzigTest(out workToBeDone);
            }
            else if (radixsortButton.Checked)
            {
                testThread = RunRadixsort(out workToBeDone);
            }
            else if (quicksortButton.Checked)
            {
                testThread = RunQuicksort(out workToBeDone);
            }
            else
            {
                throw new Exception("Not yet implemented.");
            }

            progressBar.Minimum = 0;
            progressBar.Maximum = workToBeDone;

            var uiThread = new Thread(delegate()
            {
                var startTime = System.Diagnostics.Stopwatch.StartNew();
                
                while (testThread.IsAlive)
                {
                    progressBar.Invoke((MethodInvoker) delegate()
                    {
                        progressBar.Value = Math.Min(progressBar.Maximum, Math.Max(progressBar.Minimum, Thread.VolatileRead(ref progress)));
                        progressLabel.Text = Thread.VolatileRead(ref progress).ToString() + "/" + workToBeDone.ToString();
                    });
                    Thread.Sleep(100);
                }

                progressBar.Invoke((MethodInvoker) delegate()
                {
                    progressBar.Value = progressBar.Maximum;
                    progressLabel.Text = Thread.VolatileRead(ref progress).ToString() + "/" + workToBeDone.ToString();
                });

                float durationInSeconds = startTime.ElapsedMilliseconds / 1000.0f;
                string operationCount = (!Scalar.IsCountingFlops) ? "" : "\nOperation count: " + Scalar.FLOPs;
                System.Windows.Forms.MessageBox.Show("Done.  Took " + durationInSeconds + " seconds." + operationCount);
            });

            testThread.IsBackground = true;
            testThread.Start();

            uiThread.IsBackground = true;
            uiThread.Start();
        }

        private Azimuth.DenseLinearAlgebra.RectangularMatrix RandomMatrix(int size)
        {
            return RandomMatrix(size, new Random());
        }

        private Azimuth.DenseLinearAlgebra.RectangularMatrix RandomMatrix(int size, Random rand)
        {
            var matrix = new Azimuth.DenseLinearAlgebra.RectangularMatrix(size, size);
            for (int i = 0; i < matrix.RowCount; ++i)
            {
                for (int j = 0; j < matrix.ColumnCount; ++j)
                {
                    matrix[i, j] = (double)(rand.Next(-10, 10));
                }
            }

            return matrix;
        }

        private Azimuth.SparseLinearAlgebra.SparseRectangularMatrix RandomSparseMatrix(
            int rows, int cols, int nonZeroEntries)
        {
            Random rand = new Random();

            var matrix = new Azimuth.SparseLinearAlgebra.SparseRectangularMatrix(rows, cols);

            for(int i = 0; i < nonZeroEntries; ++i)
            {
                int col = rand.Next(cols);
                int row = rand.Next(rows);

                matrix[row, col] = rand.NextDouble() * 100000;
            }

            return matrix;
        }

        private Thread RunAtan2(out int doneValue)
        {
            doneValue = 1000;

            var matrix1 = RandomMatrix(doneValue);

            var workerThread = new Thread(delegate()
            {
                for (int i = 0; i < 1000 - 1; ++i)
                {
                    var col = matrix1.ExtractColumn(i);
                    DenseVector.Atan2(col, matrix1.ExtractColumn(i + 1), ref col);
                }
            });

            return workerThread;
        }

        private Thread RunMatMul(out int doneValue)
        {
            /*
             * Here's Matlab/Octave code for comparison:
                rand1 = rand(1000,1000);
                rand2 = rand(1000,1000);
                t = cputime; mat*mat2; e = cputime - t
             * 
             * Octave w/out SSE took about 1.2 secs on a i7 @ 2.67 Ghz
             * (that is, each CPU tick can produce ~1 Flop,
             *  and we need 2 * N^3 Flops total).
             * */
            doneValue = 1000;

            var matrix1 = RandomMatrix(doneValue);
            var matrix2 = RandomMatrix(doneValue);

            var workerThread = new Thread(delegate()
            {
                var mat = matrix1.MatrixMultiply(matrix2, ref progress);
            });

            return workerThread;
        }

        private Thread RunQuicksort(out int doneValue)
        {
            doneValue = 4;
            //doneValue *= doneValue;

            var list = new List<uint>(10000000);

            for (int i = 0; i < 10000000; ++i)
            {
                list.Add((uint)rand.Next());
            }

            var workerThread = new Thread(delegate()
            {
                for (int i = 0; i < 4; ++i)
                {
                    var clone = list.ToList();
                    clone.Sort();
                    System.Threading.Interlocked.Increment(ref progress);
                }                
            });

            return workerThread;
        }

        private Thread RunRadixsort(out int doneValue)
        {
            int work = 4;
            doneValue = work;

            var list = new uint[10000000];

            for (int i = 0; i < 10000000; ++i)
            {
                list[i] = ((uint)rand.Next());
            }

            var workerThread = new Thread(delegate()
            {
                for (int i = 0; i < work; ++i)
                {
                    var newList = Azimuth.DataStructures.RadixSort.Sort(list.ToArray());
                    System.Threading.Interlocked.Increment(ref progress);
                }

                progress = work;
            });

            return workerThread;
        }

        private Thread DantzigTest(out int doneValue)
        {
            int size = 9;
            int iterCount = 10000000;
            doneValue = iterCount;

            int unsolvable = 0;

            var mat1 = RandomPSDSymmetricMatrix(size);

            var rand = new Random();

            var q = new DenseLinearAlgebra.DenseVector(size);
            for (int j = 0; j < size; ++j)
            {
                q[j] = -1;// (-2 * rand.NextDouble()) * 100000;
            }

            var workerThread = new Thread(delegate()
            {
                for (int i = 0; i < iterCount; ++i)
                {
                    if (iterCount > 1)
                    {
                        mat1 = RandomPSDSymmetricMatrix(size, rand);
                    }

                    DenseVector z;

                    bool solved = DantzigsAlgorithm.SolveLCP(out z, mat1, q);

                    if (solved)
                    {
                        var w = z.PreMatrixMultiplyBy(mat1);
                        w.Add(q);

                        for (int j = 0; j < w.Length; ++j)
                        {
                            if (w[j] < 1e-5)
                            {
                                w[j] = 0;
                            }
                        }

                        for (int j = 0; j < w.Length; ++j)
                        {
                            if (w[j].Value < -1e-5)
                            {
                                throw new Exception("Not complementary!");
                            }
                            if (!(z[j] >= 0))
                            {
                                throw new Exception("Not complementary!");
                            }
                            if (Math.Abs(w[j] * z[j]) > 1e-5f)
                            {
                                throw new Exception("Not complementary!");
                            }
                        }
                    }
                    else
                    {
                        EnsureNoAnswer(mat1);
                        unsolvable++;
                    }

                    System.Threading.Interlocked.Increment(ref progress);
                }

                MessageBox.Show(String.Format("{0} unsolvable", (Scalar)unsolvable / (Scalar)iterCount));
            });

            return workerThread;
        }

        private void EnsureNoAnswer(SquareMatrix M)
        {
            var q = new DenseVector(M.Size);

            for (int i = 0; i < q.Length; ++i)
            {
                q[i] = -1;
            }

            for (int i = 0; i < (1 << M.Size); ++i)
            {
                var pivots = new HashSet<int>();

                var Mclone = M.Clone();
                var qclone = q.Clone();

                for (int j = 0; j < M.Size; ++j)
                {
                    bool notPivoting = (i & (1 << j)) == 0;

                    if (notPivoting)
                    {
                        continue;
                    }
                    else
                    {
                        if (Mclone[j,j] == 0)
                        {
                            // This pivoting results in a singular matrix
                            break;
                        }
                        else
                        {
                            Azimuth.DenseLinearAlgebra.DantzigsAlgorithm.PrincipalPivot(ref Mclone, ref qclone, j, pivots);
                        }
                    }
                }

                if (q.All(x => x >= 0))
                {
                    // Found an answer
                    throw new Exception("Found an answer");
                }
            }
        }

        private Thread RunLLSInequalitiesTest(out int doneValue)
        {
            const int problemSize = 1024;
            doneValue = problemSize * 3;
            progress = 0;

            Random rand = new Random(8741);

            var matrix = RandomMatrix(problemSize, rand);

            var constraints = RandomMatrix(problemSize, rand);

            var B = RandomMatrix(problemSize, rand);

            var b = new DenseLinearAlgebra.DenseVector(problemSize);
            for(int i = 0; i < problemSize; ++i)
            {
                b[i] = rand.NextDouble() * 100000;
            }

            var h = new DenseLinearAlgebra.DenseVector(problemSize);
            for (int i = 0; i < problemSize; ++i)
            {
                h[i] = rand.NextDouble() * 100000;
            }

            var workerThread = new Thread(delegate()
            {
                var lls = new Azimuth.DenseLinearAlgebra.LinearLeastSquares(matrix, ref progress);

                DenseLinearAlgebra.DenseVector vector;
                lls.SolveSystemWithConstraints(out vector, b, constraints, h, ref progress);
            });

            return workerThread;
        }

        private Thread RunQRTest(out int doneValue)
        {            
            doneValue = 1000;
            progress = 0;

            var matrix = RandomMatrix(doneValue);
            
            var workerThread = new Thread(delegate()
            {
                Azimuth.DenseLinearAlgebra.QRDecomposition qr = 
                    new Azimuth.DenseLinearAlgebra.QRDecomposition(matrix, ref progress);
            });

            return workerThread;
        }

        private Thread RunLLSTest(out int doneValue)
        {
            int size = 1024;
            doneValue = size * 2;
            progress = 0;

            var matrix = RandomMatrix(size);
            
            var workerThread = new Thread(delegate()
            {
                var lls = new Azimuth.DenseLinearAlgebra.LinearLeastSquares(matrix, ref progress);
            });

            return workerThread;
        }

        private Thread RunSparseMultiply(out int doneValue)
        {
            doneValue = 1000;

            var sparseMatrix1 = RandomSparseMatrix(doneValue, doneValue, doneValue);
            var sparseMatrix2 = RandomSparseMatrix(doneValue, doneValue, doneValue);

            var workerThread = new Thread(delegate()
            {
                sparseMatrix1.MatrixMultiplyIntoSparse(sparseMatrix2, ref progress);
            });

            return workerThread;
        }

        private Thread RunCholesky(out int doneValue)
        {
            doneValue = 1000;
            var mat1 = RandomPSDSymmetricMatrix(doneValue);

            var workerThread = new Thread(delegate()
            {
                var chol = new Azimuth.DenseLinearAlgebra.CholeskyDecomposition(
                    mat1, ref progress);
            });

            return workerThread;
        }

        private SquareMatrix RandomPSDSymmetricMatrix(int size)
        {
            return RandomPSDSymmetricMatrix(size, new Random());
        }

        private SquareMatrix RandomPSDSymmetricMatrix(int size, Random rand)
        {
            var mat1 = RandomMatrix(size, rand);

            return new SquareMatrix(mat1.MatrixMultiply(mat1.Transpose()));
        }

        private Thread RunBisectionTest(out int doneValue)
        {
            doneValue = 10000;
            progress = 0;

            var workerThread = new Thread(delegate()
            {
                for (progress = 0; progress < 10000; ++progress)
                {
                    Azimuth.RootFinding.Bisection.FindRoot(x => (x - 4) * (x + 4) * (x - 3) * (x + 2) * (x + 5),
                        new Interval(-10100, 45678));
                }
            });

            return workerThread;
        }
    }
}