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

using System.Drawing;

namespace LeastSquares
{
    public class CurveFunctions
    {
        // The function.
        public static double F(List<double> coeffs, double x)
        {
            double total = 0;
            double xFactor = 1;
            for (int i = 0; i < coeffs.Count; i++)
            {
                total += xFactor * coeffs[i];
                xFactor *= x;
            }
            return total;
        }

        // Return the error squared.
        public static double ErrorSquared(List<PointF> points, List<double> coeffs)
        {
            double total = 0;
            foreach (PointF pt in points)
            {
                double dy = pt.Y - F(coeffs, pt.X);
                total += dy * dy;
            }
            return total;
        }

        // Find the least squares linear fit.
        public static List<double> FindPolynomialLeastSquaresFit(List<PointF> points, int degree)
        {
            // Allocate space for (degree + 1) equations with 
            // (degree + 2) terms each (including the constant term).
            double[,] coeffs = new double[degree + 1, degree + 2];

            // Calculate the coefficients for the equations.
            for (int j = 0; j <= degree; j++)
            {
                // Calculate the coefficients for the jth equation.

                // Calculate the constant term for this equation.
                coeffs[j, degree + 1] = 0;
                foreach (PointF pt in points)
                {
                    coeffs[j, degree + 1] -= Math.Pow(pt.X, j) * pt.Y;
                }

                // Calculate the other coefficients.
                for (int aSub = 0; aSub <= degree; aSub++)
                {
                    // Calculate the dth coefficient.
                    coeffs[j, aSub] = 0;
                    foreach (PointF pt in points)
                    {
                        coeffs[j, aSub] -= Math.Pow(pt.X, aSub + j);
                    }
                }
            }

            // Solve the equations.
            double[] answer = GaussianElimination(coeffs);

            // Return the result converted into a List<double>.
            return answer.ToList<double>();
        }

        // Perform Gaussian elimination on these coefficients.
        // Return the array of values that gives the solution.
        private static double[] GaussianElimination(double[,] coeffs)
        {
            int maxEquation = coeffs.GetUpperBound(0);
            int maxCoeff = coeffs.GetUpperBound(1);
            for (int i = 0; i <= maxEquation; i++)
            {
                // Use coeffs[i, i] to eliminate the ith
                // coefficient in all of the other equations.

                // Find a row with non-zero ith coefficient.
                if (coeffs[i, i] == 0)
                {
                    for (int j = i + 1; j <= maxEquation; j++)
                    {
                        // See if this one works.
                        if (coeffs[j, i] != 0)
                        {
                            // This one works. Swap equations i and j.
                            // This starts at k = i because all
                            // coefficients to the left are 0.
                            for (int k = i; k <= maxCoeff; k++)
                            {
                                double temp = coeffs[i, k];
                                coeffs[i, k] = coeffs[j, k];
                                coeffs[j, k] = temp;
                            }
                            break;
                        }
                    }
                }

                // Make sure we found an equation with
                // a non-zero ith coefficient.
                double coeff_i_i = coeffs[i, i];
                if (coeff_i_i == 0)
                {
                    throw new ArithmeticException(String.Format(
                        "There is no unique solution for these points.",
                        coeffs.GetUpperBound(0) - 1));
                }

                // Normalize the ith equation.
                for (int j = i; j <= maxCoeff; j++)
                {
                    coeffs[i, j] /= coeff_i_i;
                }

                // Use this equation value to zero out
                // the other equations' ith coefficients.
                for (int j = 0; j <= maxEquation; j++)
                {
                    // Skip the ith equation.
                    if (j != i)
                    {
                        // Zero the jth equation's ith coefficient.
                        double coef_j_i = coeffs[j, i];
                        for (int d = 0; d <= maxCoeff; d++)
                        {
                            coeffs[j, d] -= coeffs[i, d] * coef_j_i;
                        }
                    }
                }
            }

            // At this point, the ith equation contains
            // 2 non-zero entries:
            //      The ith entry which is 1
            //      The last entry coeffs[maxCoeff]
            // This means Ai = coeffs[maxCoeff].
            double[] solution = new double[maxEquation + 1];
            for (int i = 0; i <= maxEquation; i++)
            {
                solution[i] = coeffs[i, maxCoeff];
            }

            // Return the solution values.
            return solution;
        }
    }
}
