/// <summary>
 /// Computes the new state x1 at time t1.
 /// </summary>
 /// <param name="x0">The state x0 at time t0.</param>
 /// <param name="t0">The time t0.</param>
 /// <param name="t1">The target time t1 for which the new state x1 is computed.</param>
 /// <returns>The new state x1 at time t1.</returns>
 public override VectorD Integrate(VectorD x0, double t0, double t1)
 {
     double dt = (t1 - t0);
       VectorD d = FirstOrderDerivative(x0, t0);
       VectorD result = x0 + dt * d;
       return result;
 }

        //--------------------------------------------------------------
        /// <summary>
        /// Creates the principal component analysis for the given list of points.
        /// </summary>
        /// <param name="points">
        /// The list of data points. All points must have the same 
        /// <see cref="VectorD.NumberOfElements"/>.
        /// </param>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="points"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// <paramref name="points"/> is empty.
        /// </exception>
        public PrincipalComponentAnalysisD(IList<VectorD> points)
        {
            if (points == null)
            throw new ArgumentNullException("points");
              if (points.Count == 0)
            throw new ArgumentException("The list of points is empty.");

              // Compute covariance matrix.
              MatrixD covarianceMatrix = StatisticsHelper.ComputeCovarianceMatrix(points);

              // Perform Eigenvalue decomposition.
              EigenvalueDecompositionD evd = new EigenvalueDecompositionD(covarianceMatrix);

              int numberOfElements = evd.RealEigenvalues.NumberOfElements;
              Variances = new VectorD(numberOfElements);
              V = new MatrixD(numberOfElements, numberOfElements);

              // Sort eigenvalues by decreasing value.
              // Since covarianceMatrix is symmetric, we have no imaginary eigenvalues.
              for (int i = 0; i < Variances.NumberOfElements; i++)
              {
            int index = evd.RealEigenvalues.IndexOfLargestElement;

            Variances[i] = evd.RealEigenvalues[index];
            V.SetColumn(i, evd.V.GetColumn(index));

            evd.RealEigenvalues[index] = double.NegativeInfinity;
              }
        }

        public void TestCorrelation2()
        {
            var data = new VectorD[] {
                new VectorD(6, 10),
                new VectorD(10, 13),
                new VectorD(6, 8),
                new VectorD(10, 15),
                new VectorD(5, 8),
                new VectorD(3, 6),
                new VectorD(5, 9),
                new VectorD(9, 10),
                new VectorD(3, 7),
                new VectorD(3, 3),
                new VectorD(11, 18),
                new VectorD(6, 14),
                new VectorD(11, 18),
                new VectorD(9, 11),
                new VectorD(7, 12),
                new VectorD(5, 5),
                new VectorD(8, 7),
                new VectorD(7, 12),
                new VectorD(7, 7),
                new VectorD(9, 7)
            };
            var cor = BaseStatistics.Correlation(data, 0, 1);

            Assert.AreEqual<double>(0.749659, Math.Round(cor, 6));
        }

        public bool FindIntersection(Line2D other, out VectorD intersection)
        {
            /* http://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect
             * Now there are five cases:
             * If r × s = 0 and (q − p) × r = 0, then the two lines are collinear. If in addition, either 0 ≤ (q − p) · r ≤ r · r or 0 ≤ (p − q) · s ≤ s · s, then the two lines are overlapping.
             * If r × s = 0 and (q − p) × r = 0, but neither 0 ≤ (q − p) · r ≤ r · r nor 0 ≤ (p − q) · s ≤ s · s, then the two lines are collinear but disjoint.
             * If r × s = 0 and (q − p) × r ≠ 0, then the two lines are parallel and non-intersecting.
             * If r × s ≠ 0 and 0 ≤ t ≤ 1 and 0 ≤ u ≤ 1, the two line segments meet at the point p + t r = q + u s.
             * Otherwise, the two line segments are not parallel but do not intersect.
             */

            // line1.Start = p
            // line2.Start = q

            var r = Stop - Start;
            var s = other.Stop - other.Start;
            var startPointVector = (other.Start - Start);

            var denom = r.Cross(s); // denom = r × s

            var firstScalar = startPointVector.Cross(s) / denom; // firstScalar = t
            var secondScalar = startPointVector.Cross(r) / denom; // secondScalar = u
            intersection = Start + (firstScalar * r);
            // ReSharper disable CompareOfFloatsByEqualityOperator
            return denom != 0d && firstScalar >= 0d && firstScalar <= 1d && secondScalar >= 0d && secondScalar <= 1d;
            // ReSharper restore CompareOfFloatsByEqualityOperator
        }

 public void Add()
 {
     VectorD v1 = new VectorD(new double[] { 1, 2, 3, 4, 5 });
       VectorD v2 = new VectorD(new double[] { 6, 7, 8, 9, 10 });
       VectorD v3 = VectorD.Add(v1, v2);
       Assert.AreEqual(new VectorD(new double[] { 7, 9, 11, 13, 15 }), v3);
 }

		public void TestVectorDProd()
		{
			var vector = new VectorD(2.5, 3.0, 4.1);
			var prod = Math.Round(vector.Prod(), 2);

			Assert.AreEqual<double>(30.75, prod);
		}

        public void Test1()
        {
            VectorD state = new VectorD(new double[]{ 0, 1 });
              VectorD result = new RungeKutta4IntegratorD(GetFirstOrderDerivatives).Integrate(state, 2, 2.5);

              Assert.AreEqual(0, result[0]);
              Assert.AreEqual(3.06103515625, result[1]);
        }

 public VectorD GetFirstOrderDerivatives(VectorD x, double t)
 {
     // A dummy function: f(x[index], t) = index * t;
       VectorD result = new VectorD(x.NumberOfElements);
       for (int i = 0; i < result.NumberOfElements; i++)
     result[i] = i*t;
       return result;
 }

		public void TestUnion()
		{
			var x = new VectorD(2, 4, 6, 8);
			var y = new VectorD(4, 8, 12, 14, 16);
			var union = x.Union(y);

			Assert.AreEqual<VectorD>(
				new VectorD(2, 4, 6, 8, 12, 14, 16),
				union);
		}

        public void TestIntersect()
        {
            var x = new VectorD(2, 4, 6, 8);
            var y = new VectorD(4, 8, 12, 14, 16);
            var intersect = x.Intersect(y);

            Assert.AreEqual<VectorD>(
                new VectorD(4, 8),
                intersect);
        }

        //--------------------------------------------------------------
        /// <summary>
        /// Creates the eigenvalue decomposition of the given matrix.
        /// </summary>
        /// <param name="matrixA">The square matrix A.</param>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="matrixA"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// <paramref name="matrixA"/> is non-square (rectangular).
        /// </exception>
        public EigenvalueDecompositionD(MatrixD matrixA)
        {
            if (matrixA == null)
            throw new ArgumentNullException("matrixA");
              if (matrixA.IsSquare == false)
            throw new ArgumentException("The matrix A must be square.", "matrixA");

              _n = matrixA.NumberOfColumns;
              _d = new VectorD(_n);
              _e = new VectorD(_n);

              _isSymmetric = matrixA.IsSymmetric;

              if (_isSymmetric)
              {
            _v = matrixA.Clone();

            // Tridiagonalize.
            ReduceToTridiagonal();

            // Diagonalize.
            TridiagonalToQL();
              }
              else
              {
            _v = new MatrixD(_n, _n);

            // Abort if A contains NaN values.
            // If we continue with NaN values, we run into an infinite loop.
            for (int i = 0; i < _n; i++)
            {
              for (int j = 0; j < _n; j++)
              {
            if (Numeric.IsNaN(matrixA[i, j]))
            {
              _e.Set(double.NaN);
              _v.Set(double.NaN);
              _d.Set(double.NaN);
              return;
            }
              }
            }

            // Storage of nonsymmetric Hessenberg form.
            MatrixD matrixH = matrixA.Clone();
            // Working storage for nonsymmetric algorithm.
            double[] ort = new double[_n];

            // Reduce to Hessenberg form.
            ReduceToHessenberg(matrixH, ort);

            // Reduce Hessenberg to real Schur form.
            HessenbergToRealSchur(matrixH);
              }
        }

        public void SolveWithDefaultInitialGuess()
        {
            MatrixD A = new MatrixD(new double[,] { { 4 } });
              VectorD b = new VectorD(new double[] { 20 });

              JacobiMethodD solver = new JacobiMethodD();
              VectorD x = solver.Solve(A, b);

              Assert.IsTrue(VectorD.AreNumericallyEqual(new VectorD(1, 5), x));
              Assert.AreEqual(2, solver.NumberOfIterations);
        }

        public void Test1()
        {
            MatrixD A = new MatrixD(new double[,] { { 4 } });
              VectorD b = new VectorD(new double[] { 20 });

              SorMethodD solver = new SorMethodD();
              VectorD x = solver.Solve(A, null, b);

              Assert.IsTrue(VectorD.AreNumericallyEqual(new VectorD(1, 5), x));
              Assert.AreEqual(2, solver.NumberOfIterations);
        }

        public void Test2()
        {
            MatrixD A = new MatrixD(new double[,] { { 1, 0 },
                                              { 0, 1 }});
              VectorD b = new VectorD(new double[] { 20, 28 });

              JacobiMethodD solver = new JacobiMethodD();
              VectorD x = solver.Solve(A, null, b);

              Assert.IsTrue(VectorD.AreNumericallyEqual(b, x));
              Assert.AreEqual(2, solver.NumberOfIterations);
        }

        public void Test1()
        {
            VectorD state = new VectorD (new double[] { 0, 1, 2, 3, 4, 5 });
              VectorD result = new MidpointIntegratorD(GetFirstOrderDerivatives).Integrate(state, 2, 2.5);

              Assert.AreEqual(0, result[0]);
              Assert.AreEqual(2.125, result[1]);
              Assert.AreEqual(4.25, result[2]);
              Assert.AreEqual(6.375, result[3]);
              Assert.AreEqual(8.5, result[4]);
              Assert.AreEqual(10.625, result[5]);
        }

        public void Test3()
        {
            MatrixD A = new MatrixD(new double[,] { { 2, 0 },
                                              { 0, 2 }});
              VectorD b = new VectorD(new double[] { 20, 28 });

              GaussSeidelMethodD solver = new GaussSeidelMethodD();
              VectorD x = solver.Solve(A, null, b);

              Assert.IsTrue(VectorD.AreNumericallyEqual(b / 2, x));
              Assert.AreEqual(2, solver.NumberOfIterations);
        }

        public void Test1()
        {
            VectorD state = new VectorD (new double[]{ 1, 2, 3, 4, 5, 6 });
              VectorD result = new ExplicitEulerIntegratorD(GetFirstOrderDerivatives).Integrate(state, 2, 2.5);

              Assert.AreEqual(1, result[0]);
              Assert.AreEqual(3, result[1]);
              Assert.AreEqual(5, result[2]);
              Assert.AreEqual(7, result[3]);
              Assert.AreEqual(9, result[4]);
              Assert.AreEqual(11, result[5]);
        }

		public void TestSolveGE()
		{
			var a = new MatrixD(new double[] { 1, 4, 4, 2, 2, 5, 3, 1, 1, -2, -3, 1, 1, 4, 1, 3 }, 4, 4);
			var b = new VectorD(new double[] { -1, -7, -12, 2 });

			var c = SLESolver.Solve(a, b, SLEAlgorithm.GE);
			Assert.AreEqual<int>(4, c.Length);
			Assert.AreEqual<double>(-2.0000, Math.Round(c[0], 4));
			Assert.AreEqual<double>(-1.0000, Math.Round(c[1], 4));
			Assert.AreEqual<double>(1.0000, Math.Round(c[2], 4));
			Assert.AreEqual<double>(2.0000, Math.Round(c[3], 4));
		}

        public void Test4()
        {
            MatrixD A = new MatrixD(new double[,] { { -12, 2 },
                                              { 2, 3 }});
              VectorD b = new VectorD(new double[] { 20, 28 });

              SorMethodD solver = new SorMethodD();
              VectorD x = solver.Solve(A, null, b);

              VectorD solution = MatrixD.SolveLinearEquations(A, b);
              Assert.IsTrue(VectorD.AreNumericallyEqual(solution, x));
        }

        public void Test5()
        {
            MatrixD A = new MatrixD(new double[,] { { -21, 2, -4, 0 },
                                              { 2, 3, 0.1, -1 },
                                              { 2, 10, 111.1, -11 },
                                              { 23, 112, 111.1, -143 }});
              VectorD b = new VectorD(new double[] { 20, 28, -12, 0.1 });

              GaussSeidelMethodD solver = new GaussSeidelMethodD();
              VectorD x = solver.Solve(A, null, b);

              VectorD solution = MatrixD.SolveLinearEquations(A, b);
              Assert.IsTrue(VectorD.AreNumericallyEqual(solution, x));
        }