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

        //--------------------------------------------------------------
        /// <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="VectorF.NumberOfElements"/>.
        /// </param>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="points"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// <paramref name="points"/> is empty.
        /// </exception>
        public PrincipalComponentAnalysisF(IList<VectorF> points)
        {
            if (points == null)
            throw new ArgumentNullException("points");
              if (points.Count == 0)
            throw new ArgumentException("The list of points is empty.");

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

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

              int numberOfElements = evd.RealEigenvalues.NumberOfElements;
              Variances = new VectorF(numberOfElements);
              V = new MatrixF(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] = float.NegativeInfinity;
              }
        }

 /// <summary>Creates a new instance of NodeDestDistInt.</summary>
 /// <param name="node">The source Node.</param>
 /// <param name="dest">The destination of the new node.</param>
 /// <param name="dist">The distance from the source node to the nearest intersection.</param>
 /// <param name="intersections">The number of intersections.</param>
 public NodeDestDistInt(Node node, VectorF dest, FInt dist, int intersections)
 {
     this.Node = node;
     this.Dest = dest;
     this.Dist = dist;
     this.Intersections = intersections;
 }

 /// <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 VectorF Integrate(VectorF x0, float t0, float t1)
 {
     float dt = (t1 - t0);
       VectorF d = FirstOrderDerivative(x0, t0);
       VectorF result = x0 + dt * d;
       return result;
 }

        public void Test1()
        {
            VectorF state = new VectorF(new float[]{ 0, 1 });
              VectorF result = new RungeKutta4IntegratorF(GetFirstOrderDerivatives).Integrate(state, 2, 2.5f);

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

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

        //--------------------------------------------------------------
        /// <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 EigenvalueDecompositionF(MatrixF 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 VectorF(_n);
              _e = new VectorF(_n);

              _isSymmetric = matrixA.IsSymmetric;

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

            // Tridiagonalize.
            ReduceToTridiagonal();

            // Diagonalize.
            TridiagonalToQL();
              }
              else
              {
            _v = new MatrixF(_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(float.NaN);
              _v.Set(float.NaN);
              _d.Set(float.NaN);
              return;
            }
              }
            }

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

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

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

        public void SolveWithDefaultInitialGuess()
        {
            MatrixF A = new MatrixF(new float[,] { { 4 } });
              VectorF b = new VectorF(new float[] { 20 });

              JacobiMethodF solver = new JacobiMethodF();
              VectorF x = solver.Solve(A, b);

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

        public void Test1()
        {
            MatrixF A = new MatrixF(new float[,] { { 4 } });
              VectorF b = new VectorF(new float[] { 20 });

              SorMethodF solver = new SorMethodF();
              VectorF x = solver.Solve(A, null, b);

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

        public void Test1()
        {
            VectorF state = new VectorF (new float[]{ 1, 2, 3, 4, 5, 6 });
              VectorF result = new ExplicitEulerIntegratorF(GetFirstOrderDerivatives).Integrate(state, 2, 2.5f);

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

        public void Test4()
        {
            MatrixF A = new MatrixF(new float[,] { { -12, 2 },
                                             { 2, 3 }});
              VectorF b = new VectorF(new float[] { 20, 28 });

              SorMethodF solver = new SorMethodF();
              VectorF x = solver.Solve(A, null, b);

              VectorF solution = MatrixF.SolveLinearEquations(A, b);
              Assert.IsTrue(VectorF.AreNumericallyEqual(solution, x));
        }

        public void Test1()
        {
            VectorF state = new VectorF (new float[] { 0, 1, 2, 3, 4, 5 });
              VectorF result = new MidpointIntegratorF(GetFirstOrderDerivatives).Integrate(state, 2, 2.5f);

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

        public void Test2()
        {
            MatrixF A = new MatrixF(new float[,] { { 1, 0 },
                                             { 0, 1 }});
              VectorF b = new VectorF(new float[] { 20, 28 });

              JacobiMethodF solver = new JacobiMethodF();
              VectorF x = solver.Solve(A, null, b);

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

        public void Test3()
        {
            MatrixF A = new MatrixF(new float[,] { { 2, 0 },
                                             { 0, 2 }});
              VectorF b = new VectorF(new float[] { 20, 28 });

              GaussSeidelMethodF solver = new GaussSeidelMethodF();
              VectorF x = solver.Solve(A, null, b);

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

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

              GaussSeidelMethodF solver = new GaussSeidelMethodF();
              VectorF x = solver.Solve(A, null, b);

              VectorF solution = MatrixF.SolveLinearEquations(A, b);
              Assert.IsTrue(VectorF.AreNumericallyEqual(solution, x));
        }

        public void Test()
        {
            // Make a random list.
              RandomHelper.Random = new Random(77);
              List<VectorF> points = new List<VectorF>();
              for (int i = 0; i < 10; i++)
              {
            var vector = new VectorF(4);
            RandomHelper.Random.NextVectorF(vector, -1, 10);
            points.Add(vector);
              }

              PrincipalComponentAnalysisF pca = new PrincipalComponentAnalysisF(points);

              Assert.Greater(pca.Variances[0], pca.Variances[1]);
              Assert.Greater(pca.Variances[1], pca.Variances[2]);
              Assert.Greater(pca.Variances[2], pca.Variances[3]);
              Assert.Greater(pca.Variances[3], 0);

              Assert.IsTrue(pca.V.GetColumn(0).IsNumericallyNormalized);
              Assert.IsTrue(pca.V.GetColumn(1).IsNumericallyNormalized);
              Assert.IsTrue(pca.V.GetColumn(2).IsNumericallyNormalized);
              Assert.IsTrue(pca.V.GetColumn(3).IsNumericallyNormalized);

              // Compute covariance matrix and check if it is diagonal in the transformed space.
              MatrixF cov = StatisticsHelper.ComputeCovarianceMatrix(points);
              MatrixF transformedCov = pca.V.Transposed * cov * pca.V;

              for (int row = 0; row < transformedCov.NumberOfRows; row++)
            for (int column = 0; column < transformedCov.NumberOfColumns; column++)
              if (row != column)
            Assert.IsTrue(Numeric.IsZero(transformedCov[row, column]));

              // The principal components must be Eigenvectors which means that multiplying with the covariance
              // matrix does only change the length!
              VectorF v0 = pca.V.GetColumn(0);
              VectorF v0Result = cov * v0;
              Assert.IsTrue(VectorF.AreNumericallyEqual(v0.Normalized, v0Result.Normalized));
              VectorF v1 = pca.V.GetColumn(1);
              VectorF v1Result = cov * v1;
              Assert.IsTrue(VectorF.AreNumericallyEqual(v1.Normalized, v1Result.Normalized));
              VectorF v2 = pca.V.GetColumn(2);
              VectorF v2Result = cov * v2;
              Assert.IsTrue(VectorF.AreNumericallyEqual(v2.Normalized, v2Result.Normalized));
              VectorF v3 = pca.V.GetColumn(3);
              VectorF v3Result = cov * v3;
              Assert.IsTrue(VectorF.AreNumericallyEqual(v3.Normalized, v3Result.Normalized));
        }

        public void Absolute()
        {
            VectorF v = new VectorF(new float[] { -1, -2, -3, -4 });
              v.Absolute();
              Assert.AreEqual(1, v[0]);
              Assert.AreEqual(2, v[1]);
              Assert.AreEqual(3, v[2]);
              Assert.AreEqual(4, v[3]);

              v = new VectorF(new float[] { 1, 2, 3, 4 });
              v.Absolute();
              Assert.AreEqual(1, v[0]);
              Assert.AreEqual(2, v[1]);
              Assert.AreEqual(3, v[2]);
              Assert.AreEqual(4, v[3]);
        }

        /// <summary>
        /// Solves the specified linear system of equations <i>Ax=b</i>.
        /// </summary>
        /// <param name="matrixA">The matrix A.</param>
        /// <param name="initialX">
        /// The initial guess for x. If this value is <see langword="null"/>, a zero vector will be used
        /// as initial guess.
        /// </param>
        /// <param name="vectorB">The vector b.</param>
        /// <returns>The solution vector x.</returns>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="matrixA"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="vectorB"/> is <see langword="null"/>.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// <paramref name="matrixA"/> is not a square matrix.
        /// </exception>
        /// <exception cref="ArgumentException">
        /// The number of elements of <paramref name="initialX"/> does not match.
        /// </exception>
        public override VectorF Solve(MatrixF matrixA, VectorF initialX, VectorF vectorB)
        {
            // TODO: We can possible improve the method by reordering after each step.
              // This can be done randomly or we sort by the "convergence" of the elements.
              // See book Physics-Based Animation.

              NumberOfIterations = 0;

              if (matrixA == null)
            throw new ArgumentNullException("matrixA");
              if (vectorB == null)
            throw new ArgumentNullException("vectorB");
              if (matrixA.IsSquare == false)
            throw new ArgumentException("Matrix A must be a square matrix.", "matrixA");
              if (matrixA.NumberOfRows != vectorB.NumberOfElements)
            throw new ArgumentException("The number of rows of A and b do not match.");
              if (initialX != null && initialX.NumberOfElements != vectorB.NumberOfElements)
            throw new ArgumentException("The number of elements of the initial guess for x and b do not match.");

              VectorF xOld = initialX ?? new VectorF(vectorB.NumberOfElements);
              VectorF xNew = new VectorF(vectorB.NumberOfElements);
              bool isConverged = false;

              // Make iterations until max iteration count or the result has converged.
              for (int i = 0; i < MaxNumberOfIterations && !isConverged; i++)
              {
            for (int j = 0; j < vectorB.NumberOfElements; j++)
            {
              float delta = 0;
              for (int k = 0; k < j; k++)
            delta += matrixA[j, k] * xOld[k];

              for (int k = j + 1; k < vectorB.NumberOfElements; k++)
            delta += matrixA[j, k] * xOld[k];

              xNew[j] = (vectorB[j] - delta) / matrixA[j, j];
            }

            // Test convergence
            isConverged = VectorF.AreNumericallyEqual(xOld, xNew, Epsilon);

            xOld = xNew.Clone();
            NumberOfIterations = i + 1;
              }

              return xNew;
        }

        /// <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 VectorF Integrate(VectorF x0, float t0, float t1)
        {
            float dt = (t1 - t0);

              VectorF d1 = FirstOrderDerivative(x0, t0);
              VectorF tmp = x0 + dt / 2 * d1;

              VectorF d2 = FirstOrderDerivative(tmp, t0 + dt / 2);
              tmp = x0 + dt / 2 * d2;

              VectorF d3 = FirstOrderDerivative(tmp, t0 + dt / 2);
              tmp = x0 + dt * d3;

              VectorF d4 = FirstOrderDerivative(tmp, t1);
              VectorF result = x0 + dt / 6 * d1 + dt / 3 * d2 + dt / 3 * d3 + dt / 6 * d4;
              return result;
        }

        /// <summary>Creates a new instance of FogOfWar.</summary>
        /// <param name="rows">The number of rows.</param>
        /// <param name="cols">The number of columns.</param>
        /// <param name="worldSize">The size of the world.</param>
        /// <param name="owner">The player that owns this fog of war.</param>
        public FogOfWar(int rows, int cols, VectorF worldSize, Player owner)
        {
            clear();

            Owner = owner;
            Grid = new VisOption[rows, cols];

            // set grid square size
            SqrSize = new VectorF(worldSize.X / (FInt)cols, worldSize.Y / (FInt)rows);
            SqrSizeHalf = SqrSize / FInt.F2;

            // begin unexplored
            resetTo(VisOption.Unexplored);

            // only use drawn if human
            drawn = (owner.Type == Player.PlayerType.Human) ? new bool[rows, cols] : null;
        }