public void TestVelocityOnPlane()
        {
            Frisbee.SimulationState st = new Frisbee.SimulationState
                                             {
                                                 VX = 1,
                                                 VY = 1,
                                                 Theta = Math.PI / 4
                                             };

            Matrix<double> transformation =
                new SparseMatrix(new [,]
                                     {
                                         {st.CosTheta, st.SinTheta*st.SinPhi, -st.SinTheta*st.CosPhi},
                                         {0, st.CosPhi, st.SinPhi},
                                         {st.SinTheta, -st.CosTheta*st.SinPhi, st.CosTheta*st.CosPhi}
                                     });

            SparseVector c3 = new SparseVector(transformation.Row(2));

            SparseVector velocity = new SparseVector(new[] { st.VX, st.VY, st.VZ });
            double velocityMagnitude = velocity.Norm(2);

            double velocityC3 = velocity.DotProduct(c3);

            Vector<double> vp = velocity.Subtract(c3.Multiply(velocityC3));
            double vpMagnitude = vp.Norm(2);
        }

        public static SparseVector Extract(this SparseVector vector, int[] rows)
        {
            var extractedVector = new SparseVector(rows.Length);

            for (int i = 0; i < rows.Length; i++)
            {
                extractedVector[i] = vector[rows[i]];
            }

            return extractedVector;
        }

 public static SparseVector MeanRowVector(this Matrix<double> matrix)
 {
     SparseVector meanVector = new SparseVector(matrix.RowCount);
     for (int i = 0; i < matrix.RowCount; i++)
     {
         SparseVector row = new SparseVector(SparseVectorStorage<double>.OfVector(matrix.Row(i).Storage));
         double elements = row.NonZerosCount;
         if (elements > 0)
         {
             meanVector[i] = row.Sum() * (1 / elements);
         }
         else
         {
             meanVector[i] = 0;
         }
     }
     return meanVector;
 }

        public static SparseVector MeanColumnVector(this Matrix<double> matrix)
        {
            SparseVector meanVector = new SparseVector(matrix.ColumnCount);
            for (int i = 0; i < matrix.ColumnCount; i++)
            {
                SparseVector column = new SparseVector(SparseVectorStorage<double>.OfVector(matrix.Column(i).Storage));
                double elements = column.NonZerosCount;
                if (elements > 0)
                {
                    double sum = column.Sum();
                    meanVector[i] = sum*(1/elements);
                }
                else
                {
                    meanVector[i] = 0;
                }

            }
            return meanVector;
        }

        /// <summary>
        /// Converts the string representation of a real sparse vector to double-precision sparse vector equivalent.
        /// A return value indicates whether the conversion succeeded or failed.
        /// </summary>
        /// <param name="value">
        /// A string containing a real vector to convert.
        /// </param>
        /// <param name="formatProvider">
        /// An <see cref="IFormatProvider"/> that supplies culture-specific formatting information about value.
        /// </param>
        /// <param name="result">
        /// The parsed value.
        /// </param>
        /// <returns>
        /// If the conversion succeeds, the result will contain a complex number equivalent to value.
        /// Otherwise the result will be <c>null</c>.
        /// </returns>
        public static bool TryParse(string value, IFormatProvider formatProvider, out SparseVector result)
        {
            bool ret;
            try
            {
                result = Parse(value, formatProvider);
                ret = true;
            }
            catch (ArgumentNullException)
            {
                result = null;
                ret = false;
            }
            catch (FormatException)
            {
                result = null;
                ret = false;
            }

            return ret;
        }

        /// <summary>
        /// Outer product of two vectors
        /// </summary>
        /// <param name="u">First vector</param>
        /// <param name="v">Second vector</param>
        /// <returns>Matrix M[i,j] = u[i]*v[j] </returns>
        /// <exception cref="ArgumentNullException">If the u vector is <see langword="null" />.</exception>
        /// <exception cref="ArgumentNullException">If the v vector is <see langword="null" />.</exception>
        public static Matrix<double> OuterProduct(SparseVector u, SparseVector v)
        {
            if (u == null)
            {
                throw new ArgumentNullException("u");
            }

            if (v == null)
            {
                throw new ArgumentNullException("v");
            }

            var matrix = new SparseMatrix(u.Count, v.Count);
            for (var i = 0; i < u._storage.ValueCount; i++)
            {
                for (var j = 0; j < v._storage.ValueCount; j++)
                {
                    matrix.At(i, j, u._storage.Values[i] * v._storage.Values[j]);
                }
            }

            return matrix;
        }

        /// <summary>
        /// Returns a negated vector.
        /// </summary>
        /// <returns>The negated vector.</returns>
        /// <remarks>Added as an alternative to the unary negation operator.</remarks>
        public override Vector Negate()
        {
            var result = new SparseVector(this.Count)
                         {
                             NonZeroValues = new double[this.NonZerosCount],
                             NonZeroIndices = new int[this.NonZerosCount],
                             NonZerosCount = this.NonZerosCount
                         };

            Buffer.BlockCopy(this.NonZeroIndices, 0, result.NonZeroIndices, 0, this.NonZerosCount * Constants.SizeOfInt);

            CommonParallel.For(
                0,
                this.NonZerosCount,
                index => result.NonZeroValues[index] = -this.NonZeroValues[index]);

            return result;
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="SparseVector"/> class by
        /// copying the values from another.
        /// </summary>
        /// <param name="other">
        /// The vector to create the new vector from.
        /// </param>
        public SparseVector(SparseVector other)
            : this(other.Count)
        {
            // Lets copy only needed data. Portion of needed data is determined by NonZerosCount value
            NonZeroValues = new double[other.NonZerosCount];
            NonZeroIndices = new int[other.NonZerosCount];
            NonZerosCount = other.NonZerosCount;

            Buffer.BlockCopy(other.NonZeroValues, 0, this.NonZeroValues, 0, other.NonZerosCount * Constants.SizeOfDouble);
            Buffer.BlockCopy(other.NonZeroIndices, 0, this.NonZeroIndices, 0, other.NonZerosCount * Constants.SizeOfInt);
        }

 /// <summary>
 /// Tensor Product (Outer) of this and another vector.
 /// </summary>
 /// <param name="v">The vector to operate on.</param>
 /// <returns>
 /// Matrix M[i,j] = this[i] * v[j].
 /// </returns>
 /// <seealso cref="OuterProduct"/>
 public Matrix TensorMultiply(SparseVector v)
 {
     return OuterProduct(this, v);
 }

 private static double CosineR(SparseVector Vector_a, SparseVector Vector_b)
 {
     return Vector_a.DotProduct(Vector_b) / (Vector_a.L2Norm() * Vector_b.L2Norm());
     //return Distance.Cosine(R.Row(a).ToArray(), R.Row(b).ToArray());
 }

        public void CheckSparseMechanismByZeroMultiply()
        {
            var vector = new SparseVector(10000);

            // Add non-zero elements
            vector[200] = 1.5;
            vector[500] = 3.5;
            vector[800] = 5.5;
            vector[0] = 7.5;

            // Multiply by 0
            vector *= 0;

            var storage = (SparseVectorStorage<double>) vector.Storage;
            Assert.AreEqual(0, vector[200]);
            Assert.AreEqual(0, vector[500]);
            Assert.AreEqual(0, vector[800]);
            Assert.AreEqual(0, vector[0]);
            Assert.AreEqual(0, storage.ValueCount);
        }

        public void CheckSparseMechanismBySettingValues()
        {
            var vector = new SparseVector(10000);
            var storage = (SparseVectorStorage<double>) vector.Storage;

            // Add non-zero elements
            vector[200] = 1.5;
            Assert.AreEqual(1.5, vector[200]);
            Assert.AreEqual(1, storage.ValueCount);

            vector[500] = 3.5;
            Assert.AreEqual(3.5, vector[500]);
            Assert.AreEqual(2, storage.ValueCount);

            vector[800] = 5.5;
            Assert.AreEqual(5.5, vector[800]);
            Assert.AreEqual(3, storage.ValueCount);

            vector[0] = 7.5;
            Assert.AreEqual(7.5, vector[0]);
            Assert.AreEqual(4, storage.ValueCount);

            // Remove non-zero elements
            vector[200] = 0;
            Assert.AreEqual(0, vector[200]);
            Assert.AreEqual(3, storage.ValueCount);

            vector[500] = 0;
            Assert.AreEqual(0, vector[500]);
            Assert.AreEqual(2, storage.ValueCount);

            vector[800] = 0;
            Assert.AreEqual(0, vector[800]);
            Assert.AreEqual(1, storage.ValueCount);

            vector[0] = 0;
            Assert.AreEqual(0, vector[0]);
            Assert.AreEqual(0, storage.ValueCount);
        }

        public void CanScaleAVectorWhenSettingPreviousNonzeroElementsToZero()
        {
            var vector = new SparseVector(20);
            vector[10] = 1.0;
            vector[11] = 2.0;
            vector[11] = 0.0;

            var scaled = new SparseVector(20);
            vector.Multiply(3.0, scaled);

            Assert.AreEqual(3.0, scaled[10]);
            Assert.AreEqual(0.0, scaled[11]);
        }

        public void CanPointwiseMultiplySparseVector()
        {
            var zeroArray = new[] {0.0, 1.0, 0.0, 1.0, 0.0};
            var vector1 = SparseVector.OfEnumerable(Data);
            var vector2 = SparseVector.OfEnumerable(zeroArray);
            var result = new SparseVector(vector1.Count);

            vector1.PointwiseMultiply(vector2, result);

            for (var i = 0; i < vector1.Count; i++)
            {
                Assert.AreEqual(Data[i]*zeroArray[i], result[i]);
            }

            var resultStorage = (SparseVectorStorage<double>) result.Storage;
            Assert.AreEqual(2, resultStorage.ValueCount);
        }

        public void CanDotProductOfTwoSparseVectors()
        {
            var vectorA = new SparseVector(10000);
            vectorA[200] = 1;
            vectorA[500] = 3;
            vectorA[800] = 5;
            vectorA[100] = 7;
            vectorA[900] = 9;

            var vectorB = new SparseVector(10000);
            vectorB[300] = 3;
            vectorB[500] = 5;
            vectorB[800] = 7;

            Assert.AreEqual(50.0, vectorA.DotProduct(vectorB));
        }

        /// <summary>
        /// Creates a vector containing specified elements.
        /// </summary>
        /// <param name="index">The first element to begin copying from.</param>
        /// <param name="length">The number of elements to copy.</param>
        /// <returns>A vector containing a copy of the specified elements.</returns>
        /// <exception cref="ArgumentOutOfRangeException"><list><item>If <paramref name="index"/> is not positive or
        /// greater than or equal to the size of the vector.</item>
        /// <item>If <paramref name="index"/> + <paramref name="length"/> is greater than or equal to the size of the vector.</item>
        /// </list></exception>
        /// <exception cref="ArgumentException">If <paramref name="length"/> is not positive.</exception>
        public override Vector<double> SubVector(int index, int length)
        {
            if (index < 0 || index >= Count)
            {
                throw new ArgumentOutOfRangeException("index");
            }

            if (length <= 0)
            {
                throw new ArgumentOutOfRangeException("length");
            }

            if (index + length > Count)
            {
                throw new ArgumentOutOfRangeException("length");
            }

            var result = new SparseVector(length);
            for (var i = index; i < index + length; i++)
            {
                result.At(i - index, At(i));
            }

            return result;
        }

 private static double PearsonR(SparseVector Vector_a, SparseVector Vector_b)
 {
     double correlation = Correlation.Pearson(Vector_a,Vector_b);
     if (double.IsNaN(correlation))
     {
         // This means one of the row has 0 standard divation,
         // it does not correlate to anyone
         // so I assign the correlatino to be 0. however, strictly speaking, it should be left NaN
         correlation = 0;
     }
     return correlation;
 }

        /// <summary>
        /// Creates a new instance of the Vector class.
        /// </summary>
        /// <param name="data">The array to create this vector from.</param>
        /// <returns>The new <c>Vector</c>.</returns>
        protected override Vector<double> CreateVector(IList<double> data)
        {
            var vector = new SparseVector(data.Count);
            for (var index = 0; index < data.Count; index++)
            {
                vector[index] = data[index];
            }

            return vector;
        }

        /// <summary>
        /// Creates a vector containing specified elements.
        /// </summary>
        /// <param name="index">The first element to begin copying from.</param>
        /// <param name="length">The number of elements to copy.</param>
        /// <returns>A vector containing a copy of the specified elements.</returns>
        /// <exception cref="ArgumentOutOfRangeException"><list><item>If <paramref name="index"/> is not positive or
        /// greater than or equal to the size of the vector.</item>
        /// <item>If <paramref name="index"/> + <paramref name="length"/> is greater than or equal to the size of the vector.</item>
        /// </list></exception>
        /// <exception cref="ArgumentException">If <paramref name="length"/> is not positive.</exception>
        public override Vector SubVector(int index, int length)
        {
            if (index < 0 || index >= this.Count)
            {
                throw new ArgumentOutOfRangeException("index");
            }

            if (length <= 0)
            {
                throw new ArgumentOutOfRangeException("length");
            }

            if (index + length > this.Count)
            {
                throw new ArgumentOutOfRangeException("length");
            }

            var result = new SparseVector(length);
            for (int i = index; i < index + length; i++)
                result[i - index] = this[i];

            return result;
        }

 public void CanCreateSparseMatrix()
 {
     var vector = new SparseVector(3);
     var matrix = Matrix<double>.Build.SameAs(vector, 2, 3);
     Assert.IsInstanceOf<SparseMatrix>(matrix);
     Assert.AreEqual(2, matrix.RowCount);
     Assert.AreEqual(3, matrix.ColumnCount);
 }