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);
        }

        /// <summary>
        /// Learn the current input.
        /// </summary>
        /// <param name="input"></param>
        /// <exception cref="HtmRuleException"></exception>
        public override void Learn(SparseMatrix input)
        {
            // Limitation due to HTM v1.x design (cannot learn after inference without modifying ouput vector structure)
            if (!IsLearning)
                // TODOlater allow learning after training when using FixedMaxSize nodes
                throw new HtmRuleException("Cannot learn after any other mode than learning", this);

            // Ignore blank input
            //TODOlater? treat any input with identical values for *all* components as blank?
            //TODOlater use DetectBlanks/DetectBlanksMode properties
            if (input.NonZerosCount == 0)
            { return; }

            // Check matrix size
            if (CoincidencesFrequencies.Count > 0 &&
                (CoincidencesFrequencies.Keys.First().RowCount != input.RowCount || CoincidencesFrequencies.Keys.First().ColumnCount != input.ColumnCount))
                throw new HtmRuleException("Cannot learn varying sized inputs", this);

            SparseMatrix existingCoincidence = FindClosestCoincidence(input);

            if (existingCoincidence != null)
                CoincidencesFrequencies[existingCoincidence] += 1;
            else
                if (CoincidencesFrequencies.Count < MaxOutputSize)
                    CoincidencesFrequencies[input] = 1;
        }

        public void ErrorWhenLearningAfterInferenceMode()
        {
            var node = new SpatialNode2DGaussian();
            var mat = new SparseMatrix(4, 4, 4.0);
            var learnAfterInferFails = false;
            var learnAfterTBInferFails = false;
            node.Learn(mat);

            node.Infer(mat);
            try
            {
                node.Learn(mat);
            }
            catch (HtmRuleException e)
            {
                learnAfterInferFails = true;
                Debug.WriteLine(e.Message);
            }

            node.TimeInfer(mat);
            try
            {
                node.Learn(mat);
            }
            catch (HtmRuleException e)
            {
                learnAfterTBInferFails = true;
                Debug.WriteLine(e.Message);
            }

            Assert.IsTrue(learnAfterInferFails);
            Assert.IsTrue(learnAfterTBInferFails);
        }

        public void SparseMatrixSubMatrixMethodWorks()
        {
            var matrix = new SparseMatrix(10, 10, 1.0);
            var sub = matrix.SubMatrix(8, 2, 0, 2);

            Assert.AreEqual(2, sub.RowCount);
            Assert.AreEqual(2, sub.ColumnCount);
        }

        public void SolveLongMatrixThrowsArgumentException()
        {
            var matrix = new SparseMatrix(3, 2);
            var input = new DenseVector(3);

            var solver = new GpBiCg();
            Assert.Throws<ArgumentException>(() => matrix.SolveIterative(input, solver));
        }

        public void SolveLongMatrixThrowsArgumentException()
        {
            var matrix = new SparseMatrix(3, 2);
            var input = new DenseVector(3);

            var solver = new MlkBiCgStab();
            Assert.That(() => matrix.SolveIterative(input, solver), Throws.ArgumentException);
        }

        public void SolveWideMatrixThrowsArgumentException()
        {
            var matrix = new SparseMatrix(2, 3);
            var input = new DenseVector(2);

            var solver = new TFQMR();
            Assert.That(() => matrix.SolveIterative(input, solver), Throws.ArgumentException);
        }

        public void SolveWideMatrixThrowsArgumentException()
        {
            var matrix = new SparseMatrix(2, 3);
            var input = new DenseVector(2);

            var solver = new MlkBiCgStab();
            Assert.Throws<ArgumentException>(() => matrix.SolveIterative(input, solver));
        }

 /// <summary>
 /// Create unit matrix.
 /// </summary>
 /// <param name="size">Matrix size.</param>
 /// <returns>New unit matrix.</returns>
 internal SparseMatrix CreateUnitMatrix(int size)
 {
     var matrix = new SparseMatrix(size);
     for (var i = 0; i < size; i++)
     {
         matrix[i, i] = 2;
     }
     return matrix;
 }

        /// <summary>
        /// Check the result.
        /// </summary>
        /// <param name="preconditioner">Specific preconditioner.</param>
        /// <param name="matrix">Source matrix.</param>
        /// <param name="vector">Initial vector.</param>
        /// <param name="result">Result vector.</param>
        protected override void CheckResult(IPreconditioner<double> preconditioner, SparseMatrix matrix, Vector<double> vector, Vector<double> result)
        {
            Assert.AreEqual(typeof (UnitPreconditioner<double>), preconditioner.GetType(), "#01");

            // Unit preconditioner is doing nothing. Vector and result should be equal
            for (var i = 0; i < vector.Count; i++)
            {
                Assert.IsTrue(vector[i] == result[i], "#02-" + i);
            }
        }

        public void SparseMatrixAddsCorrectlyToZeroComponents()
        {
            var m1 = new SparseMatrix(1, 3);
            var m2 = new SparseMatrix(new double[,] { { 0, 1, 1 } });

            var sum1 = m1 + m2;
            var sum2 = m2 + m1;

            Assert.AreEqual(m2, sum2);
            Assert.AreEqual(m2, sum1);
        }

        public void SparseMatrixSubtractsCorrectlyFromZeroComponents()
        {
            var m1 = new SparseMatrix(new double[,] { { 1, 0, 0, 1 }, { 1, 1, 0, 1 } });
            var m2 = new SparseMatrix(new double[,] { { 1, 1, 1, 0 }, { 1, 0, 0, 1 } });

            var diff1 = new SparseMatrix(new double[,] { { 0, -1, -1, 1 }, { 0, 1, 0, 0 } });
            var diff2 = new SparseMatrix(new double[,] { { 0, 1, 1, -1 }, { 0, -1, 0, 0 } });

            Assert.AreEqual(diff2, m2 - m1);
            Assert.AreEqual(diff1, m1 - m2);
        }

        /// <summary>
        /// Check the result.
        /// </summary>
        /// <param name="preconditioner">Specific preconditioner.</param>
        /// <param name="matrix">Source matrix.</param>
        /// <param name="vector">Initial vector.</param>
        /// <param name="result">Result vector.</param>
        protected override void CheckResult(IPreconditioner<double> preconditioner, SparseMatrix matrix, Vector<double> vector, Vector<double> result)
        {
            Assert.AreEqual(typeof (DiagonalPreconditioner), preconditioner.GetType(), "#01");

            // Compute M * result = product
            // compare vector and product. Should be equal
            var product = new DenseVector(result.Count);
            matrix.Multiply(result, product);
            for (var i = 0; i < product.Count; i++)
            {
                Assert.IsTrue(vector[i].AlmostEqualNumbersBetween(product[i], -Epsilon.Magnitude()), "#02-" + i);
            }
        }

        public static SparseMatrix Extract(this SparseMatrix matrix, int[] rows, int[] columns)
        {
            var extractedMatrix = new SparseMatrix(rows.Length, columns.Length);

            for (int c = 0; c < columns.Length; c++)
            {
                for (int r = 0; r < rows.Length; r++)
                {
                    extractedMatrix[r, c] = matrix[rows[r], columns[c]];
                }
            }

            return extractedMatrix;
        }

        /// <summary>
        /// Do Flash inference for the given input
        /// </summary>
        /// <param name="input"></param>
        /// <returns></returns>
        public override Vector Infer(SparseMatrix input)
        {
            State = NodeState.FlashInference;

            if (LearnedCoincidences.Length == 0)
            {
                return new SparseVector(1);
            }

            var output = new DenseVector(LearnedCoincidences.Length);
            for (int i = 0; i < LearnedCoincidences.Length; ++i)
            {
                var diff = LearnedCoincidences[i] - input;
                var norm = diff.FrobeniusNorm();

                output[i] = Math.Exp(-(norm * norm)  / (2 * SquaredSigma));
            }

            return output;
        }

        /// <summary>
        /// List indexing is 2000x faster than Matrix.Row() or enumeration.
        /// </summary>
        public static void SpeedOfGetRow()
        {
            SparseMatrix myMatrix = new SparseMatrix(1000, 1000);
            SparseVector myVector = SparseVector.OfVector(Vector.Build.Random(1000));
            myVector.CoerceZero(1.8);
            for (int i = 0; i < 1000; i++)
            {
                myMatrix.SetRow(i, myVector);
            }
            List<Vector<double>> myList = new List<Vector<double>>(myMatrix.EnumerateRows());

            Utils.StartTimer();
            for (int repeat = 0; repeat < 10; repeat++)
            {
                for (int i = 0; i < 1000; i++)
                {
                    double foo = myMatrix.Row(i)[0];
                }
            }
            Utils.StopTimer();

            Utils.StartTimer();
            for (int repeat = 0; repeat < 10; repeat++)
            {
                foreach(var row in myMatrix.EnumerateRowsIndexed())
                {
                    double foo = row.Item2[0];
                }
            }
            Utils.StopTimer();

            Utils.StartTimer();
            for (int repeat = 0; repeat < 10; repeat++)
            {
                for (int i = 0; i < 1000; i++)
                {
                    double foo = myList[i][0];
                }
            }
            Utils.StopTimer();
        }

        public SparseMatrix GetMeanMatrix()
        {
            SparseMatrix meanMatix = new SparseMatrix(ratingMatrix.RowCount, ratingMatrix.ColumnCount);

            double totalMean = ratingMatrix.SumMatrix() * (1 / (double)ratingMatrix.NonZerosCount);

            SparseVector movieMeans = ratingMatrix.MeanRowVector();
            SparseVector userMeans = ratingMatrix.MeanColumnVector();

            for (int i = 0; i < meanMatix.RowCount; i++)
            {
                for (int j = 0; j < meanMatix.ColumnCount; j++)
                {
                    double cell = ratingMatrix[i, j];
                    if (cell > 0 || cell < 0)
                        meanMatix[i, j] = ratingMatrix[i, j] - movieMeans[i] - userMeans[j] + totalMean;
                }
            }

            return meanMatix;
        }

        public void CanLearnAndInfer1OnRecognizedSubInputs_5By4Layer0_5Coverage()
        {
            uint maxoutputsize = 5;
            var coverage = 0.5;
            var clone = false;
            uint layerheight = 5;
            uint layerwidth = 4;
            int inputtolayerratio = 2;

            var layer = new Spatial2DLayer(SpatialLayerType.Gaussian, layerheight, layerwidth, coverage, clone, maxoutputsize);
            var input1 = new SparseMatrix((int)layerheight * inputtolayerratio, (int)layerwidth * inputtolayerratio, 1.0);
            var input2 = 2 * input1;

            //
            layer.Learn(input1);
            layer.Learn(input2);

            var infer2 = layer.Infer(input2);
            var infer1 = layer.Infer(input1);
            var input3 = (SparseMatrix)new SparseMatrix((int)layerwidth * inputtolayerratio, (int)layerwidth * inputtolayerratio, 1.0)
                .Stack(new SparseMatrix((int)(layerheight - layerwidth) * inputtolayerratio, (int)layerwidth * inputtolayerratio, 3.0));
            var infer3 = layer.Infer(input3);

            //
            foreach (var infer in new SparseMatrix[] { infer1, infer2, infer3 })
            {
                for (int row = 0; row < layerheight; row++)
                {
                    for (int col = 0; col < layerwidth * layer.MaxNodeOutputSize; col += layer.MaxNodeOutputSize)
                    {
                        var subInput = infer.SubMatrix(row, 1, col, layer.MaxNodeOutputSize);
                        if (infer.Equals(infer3) && row >= layerheight - 2 && row < layerheight)
                            for (int i = 0; i < subInput.ColumnCount; i++)
                                Assert.IsTrue(subInput[0, i] < 0.1);
                        else
                            Assert.IsTrue(subInput[0, 0] == 1.0 || subInput[0, 1] == 1.0);
                    }
                }
            }
        }

        public void CanAddSparseMatricesBothWays()
        {
            var m1 = new SparseMatrix(1, 3);
            var m2 = Matrix<double>.Build.SparseOfArray(new double[,] { { 0, 1, 1 } });
            var sum1 = m1 + m2;
            var sum2 = m2 + m1;
            Assert.IsTrue(sum1.Equals(m2));
            Assert.IsTrue(sum1.Equals(sum2));

            Matrix<double> sparseResult = new SparseMatrix(1, 3);
            sparseResult.Add(m2, sparseResult);
            Assert.IsTrue(sparseResult.Equals(sum1));

            sparseResult = Matrix<double>.Build.SparseOfArray(new double[,] { { 0, 1, 1 } });
            sparseResult.Add(m1, sparseResult);
            Assert.IsTrue(sparseResult.Equals(sum1));

            sparseResult = Matrix<double>.Build.SparseOfArray(new double[,] { { 0, 1, 1 } });
            m1.Add(sparseResult, sparseResult);
            Assert.IsTrue(sparseResult.Equals(sum1));

            sparseResult = Matrix<double>.Build.SparseOfArray(new double[,] { { 0, 1, 1 } });
            sparseResult.Add(sparseResult, sparseResult);
            Assert.IsTrue(sparseResult.Equals(2*sum1));

            Matrix<double> denseResult = new DenseMatrix(1, 3);
            denseResult.Add(m2, denseResult);
            Assert.IsTrue(denseResult.Equals(sum1));

            denseResult = Matrix<double>.Build.DenseOfArray(new double[,] { { 0, 1, 1 } });
            denseResult.Add(m1, denseResult);
            Assert.IsTrue(denseResult.Equals(sum1));

            var m3 = Matrix<double>.Build.DenseOfArray(new double[,] { { 0, 1, 1 } });
            var sum3 = m1 + m3;
            var sum4 = m3 + m1;
            Assert.IsTrue(sum3.Equals(m3));
            Assert.IsTrue(sum3.Equals(sum4));
        }

        public void SparseMatrixIndexedEnumeratorWorks()
        {
            var cols = 7;
            var rows = 7;
            var matrix = new SparseMatrix(rows, cols, new double[] { 2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06, /**/2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06, 34.5104995998366, 3.30253638788279, 2.7461515689052E-06, 2.7461515689052E-06, /**/2.7461515689052E-06, 2.7461515689052E-06, 3.30253638788279, 255, 3.30253638788279, 2.7461515689052E-06, 2.37450036851675E-06, /**/2.7461515689052E-06, 2.7461515689052E-06, 24.4026090935093, 255, 24.4026090935093, 2.7461515689052E-06, 2.7461515689052E-06, /**/2.7461515689052E-06, 2.7461515689052E-06, 3.30253638788279, 34.5104995998366, 3.30253638788279, 2.7461515689052E-06, 2.37450036851675E-06, /**/2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06, /**/2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06, 6.62979636198309E-06, 2.7461515689052E-06, 2.7461515689052E-06, 2.7461515689052E-06 });
            var normalized = matrix.Normalize(1.0, 2.7461515689052E-06);

            //
            var count = 0;
            foreach (var el in matrix.IndexedEnumerator())
            {
                count++;
            }
            var count2 = 0;
            foreach (var el in normalized.IndexedEnumerator())
            {
                count2++;
            }

            //
            Assert.AreEqual(cols * rows, count);
            Assert.AreEqual(cols * rows, count2);
        }