C# Matrix.Column方法代码示例

 public Matrix<double> NormalizeColumns(Matrix<double> dataToNormalize, IList<int> columnsToNormalize = null)
 {
     columnsToNormalize = columnsToNormalize ?? Enumerable.Range(0, dataToNormalize.ColumnCount).ToList();
     var normalizedColumns = new ConcurrentBag<Tuple<int, Vector<double>>>();
     var columnStds = new double[dataToNormalize.ColumnCount];
     var columnMeans = new double[dataToNormalize.ColumnCount];
     for (int colIdx = 0; colIdx < dataToNormalize.ColumnCount; colIdx++)
     {
         columnStds[colIdx] = dataToNormalize.Column(colIdx).StandardDeviation();
         columnMeans[colIdx] = dataToNormalize.Column(colIdx).Mean();
     }
     Parallel.For(0, dataToNormalize.ColumnCount, colIdx =>
     {
         Vector<double> vectorToAdd;
         if (!columnsToNormalize.Contains(colIdx))
         {
             vectorToAdd = dataToNormalize.Column(colIdx);
         }
         else
         {
             var columnVector = dataToNormalize.Column(colIdx);
             var columnStd = columnStds[colIdx];
             var columnMean = columnMeans[colIdx];
             vectorToAdd = columnVector.Subtract(columnMean).Divide(columnStd);
         }
         normalizedColumns.Add(new Tuple<int, Vector<double>>(colIdx, vectorToAdd));
     });
     return Matrix<double>.Build.DenseOfColumnVectors(normalizedColumns.OrderBy(tpl => tpl.Item1).Select(tpl => tpl.Item2));
 }

        private static Tuple<DenseMatrix, DenseVector, DenseVector> FeatureNormalize(Matrix<double> x)
        {
            var mu = DenseVector.Create(x.ColumnCount, i => x.Column(i).Mean());
            var sigma = DenseVector.Create(x.ColumnCount, i => x.Column(i).StandardDeviation());

            return Tuple.Create(
                DenseMatrix.Create(x.RowCount, x.ColumnCount, (i, j) => (x[i, j] - mu[j]) / sigma[j]),
                mu,
                sigma);
        }

 static int ChooseBestFeatureToSplit(Matrix<double> dataSet, string[] labels)
 {
     var numFeatures = dataSet.ColumnCount;
     var baseEntropy = CalcShannonEnt(dataSet, labels);
     var bestInfoGain = 0.0;
     var bestFeature = -1;
     for(var i = 0;i<numFeatures;i++)
     {
         var featList = dataSet.Column(i);
         var uniqueVals = featList.Distinct();
         var newEntropy = 0.0;
         foreach(var value in uniqueVals)
         {
             var subDataSet = SplitDataSet(dataSet, i, value);
             var prob = (float)subDataSet.RowCount / (float)dataSet.RowCount;
             newEntropy += prob * CalcShannonEnt(subDataSet, labels);
         }
         var infoGain = baseEntropy - newEntropy;
         if (infoGain > bestInfoGain)
         {
             bestInfoGain = infoGain;
             bestFeature = i;
         }
     }
     return bestFeature;
 }

 public static NormDataSet ToNorm(Matrix<double> matrix)
 {
     var minArray = new double[matrix.ColumnCount];
     var rangeArray = new double[matrix.ColumnCount];
     for(var i = 0;i<matrix.ColumnCount;i++)
     {
         minArray[i] = matrix.Column(i).Minimum();
         rangeArray[i] = matrix.Column(i).Maximum() - minArray[i];
     }
     var minMatrix = Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount, (i,j)=> minArray[j]);
     var rangeMatrix = Matrix.Build.Dense(matrix.RowCount, matrix.ColumnCount, (i,j)=> rangeArray[j]);
     var normMatrix = matrix - minMatrix;
     normMatrix = Matrix.op_DotDivide(normMatrix, rangeMatrix);
     return new NormDataSet
     {
         NormData = normMatrix,
         MinVals = Vector.Build.DenseOfArray(minArray),
         Ranges = Vector.Build.DenseOfArray(rangeArray)
     };
 }

        public Molecule(double[,] inputAtoms, bool getDirections)
        {
            numAtoms = inputAtoms.GetLength(0);
            atoms = DenseMatrix.OfArray(inputAtoms);
            centroid = atoms.ColumnSums() / (double)numAtoms;

            if (getDirections)
            {
                for (int i = 0; i < 3; i++)
                {
                    atoms.SetColumn(i, atoms.Column(i) - centroid[i]);
                }
                direction = atoms.Svd(true).VT.Row(0);
                if ((atoms.Row(atoms.RowCount - 1) - atoms.Row(0)).DotProduct(direction) < 0)
                {
                    direction *= -1;
                }
            }
        }

        /// <summary>
        /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the
        /// solution matrix and X is the unknown matrix.
        /// </summary>
        /// <param name="matrix">The coefficient <see cref="Matrix"/>, <c>A</c>.</param>
        /// <param name="input">The solution <see cref="Matrix"/>, <c>B</c>.</param>
        /// <param name="result">The result <see cref="Matrix"/>, <c>X</c></param>
        public void Solve(Matrix matrix, Matrix input, Matrix result)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

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

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

            if (matrix.RowCount != input.RowCount || input.RowCount != result.RowCount || input.ColumnCount != result.ColumnCount)
            {
                throw Matrix.DimensionsDontMatch<ArgumentException>(matrix, input, result);
            }

            for (var column = 0; column < input.ColumnCount; column++)
            {
                var solution = Solve(matrix, (Vector)input.Column(column));
                foreach (var element in solution.GetIndexedEnumerator())
                {
                    result.At(element.Item1, column, element.Item2);
                }
            }
        }

        /// <summary>
        /// Multiplies this matrix with another matrix and places the results into the result matrix.
        /// </summary>
        /// <param name="other">The matrix to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        protected override void DoMultiply(Matrix<Complex32> other, Matrix<Complex32> result)
        {
            var sparseOther = other as SparseMatrix;
            var sparseResult = result as SparseMatrix;
            if (sparseOther != null && sparseResult != null)
            {
                DoMultiplySparse(sparseOther, sparseResult);
                return;
            }

            var diagonalOther = other.Storage as DiagonalMatrixStorage<Complex32>;
            if (diagonalOther != null && sparseResult != null)
            {
                var diagonal = diagonalOther.Data;
                if (other.ColumnCount == other.RowCount)
                {
                    Storage.MapIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], Zeros.AllowSkip, ExistingData.Clear);
                }
                else
                {
                    result.Storage.Clear();
                    Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], 0, 0, RowCount, 0, 0, ColumnCount, Zeros.AllowSkip, ExistingData.AssumeZeros);
                }
                return;
            }

            result.Clear();
            var columnVector = new DenseVector(other.RowCount);

            var rowPointers = _storage.RowPointers;
            var columnIndices = _storage.ColumnIndices;
            var values = _storage.Values;

            for (var row = 0; row < RowCount; row++)
            {
                var startIndex = rowPointers[row];
                var endIndex = rowPointers[row + 1];

                if (startIndex == endIndex)
                {
                    continue;
                }

                for (var column = 0; column < other.ColumnCount; column++)
                {
                    // Multiply row of matrix A on column of matrix B
                    other.Column(column, columnVector);

                    var sum = Complex32.Zero;
                    for (var index = startIndex; index < endIndex; index++)
                    {
                        sum += values[index] * columnVector[columnIndices[index]];
                    }

                    result.At(row, column, sum);
                }
            }
        }

        /// <summary>
        /// Multiplies this matrix with another matrix and places the results into the result matrix.
        /// </summary>
        /// <param name="other">The matrix to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        protected override void DoMultiply(Matrix<double> other, Matrix<double> result)
        {
            result.Clear();
            var columnVector = new DenseVector(other.RowCount);

            var rowPointers = _storage.RowPointers;
            var columnIndices = _storage.ColumnIndices;
            var values = _storage.Values;
            var valueCount = _storage.ValueCount;

            for (var row = 0; row < RowCount; row++)
            {
                // Get the begin / end index for the current row
                var startIndex = rowPointers[row];
                var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
                if (startIndex == endIndex)
                {
                    continue;
                }

                for (var column = 0; column < other.ColumnCount; column++)
                {
                    // Multiply row of matrix A on column of matrix B
                    other.Column(column, columnVector);

                    var sum = 0d;
                    for (var index = startIndex; index < endIndex; index++)
                    {
                        sum += values[index] * columnVector[columnIndices[index]];
                    }

                    result.At(row, column, sum);
                }
            }
        }

        void UpdateCameraMatrix(Matrix<double> camera, CameraIndex index)
        {
            var RQ = camera.SubMatrix(0, 3, 0, 3).QR();

            var calib = RQ.R;
            if(Math.Abs(calib[2, 2] - 1) > 1e-6)
            {
                double scale = calib[2, 2];
                camera.MultiplyThis(scale);
                // NotifyPropertyChanged("CameraLeft");
                RQ = camera.SubMatrix(0, 3, 0, 3).QR();
            }
            calib = RQ.R;
            var rot = RQ.Q;

            // If fx < 0 then set fx = -fx and [r11,r12,r13] = -[r11,r12,r13]
            // As first row of rotation matrix is multiplied only with fx, then changing sign of both
            // fx and this row won't change matrix M = K*R, and so camera matrix (same with fy, but we need to change skew also)
            if(calib[0, 0] < 0)
            {
                calib[0, 0] = -calib[0, 0];
                rot[0, 0] = -rot[0, 0];
                rot[0, 1] = -rot[0, 1];
                rot[0, 2] = -rot[0, 2];
            }
            if(calib[1, 1] < 0)
            {
                calib[1, 1] = -calib[1, 1];
                calib[0, 1] = -calib[0, 1];
                rot[1, 0] = -rot[1, 0];
                rot[1, 1] = -rot[1, 1];
                rot[1, 2] = -rot[1, 2];
            }

            var trans = -camera.SubMatrix(0, 3, 0, 3).Inverse().Multiply(camera.Column(3));

            SetCalibrationMatrix(index, calib);
            SetRotationMatrix(index, rot);
            SetTranslationVector(index, trans);

            ComputeEssentialFundamental();
        }

        public void Test_EstimateCameraMatrix_Minimised()
        {
            PrepareCameraMatrix();
            var points = GenerateCalibrationPoints_Random(100);
            var noisedPoints = AddNoise(points, _varianceReal, _varianceImage, 200);
            PrepareCalibrator(noisedPoints);

            _calib.HomoPoints();
            _calib.NormalizeImagePoints();
            _calib.NormalizeRealPoints();
            _calib._pointsNormalised = true;
            _calib.CameraMatrix = _calib.FindLinearEstimationOfCameraMatrix();
            //    _calib.FindNormalisedVariances();
            _calib.DenormaliseCameraMatrix();
            _calib._pointsNormalised = false;

            _eCM = _calib.CameraMatrix;
            double totalDiff = 0.0;
            for(int p = 0; p < _pointsCount; ++p)
            {
                var cp = points[p];
                Vector<double> rp = new DenseVector(4);
                rp[0] = cp.RealX;
                rp[1] = cp.RealY;
                rp[2] = cp.RealZ;
                rp[3] = 1.0;
                var imagePoint = _eCM * rp;

                Vector2 ip = new Vector2(imagePoint[0] / imagePoint[2], imagePoint[1] / imagePoint[2]);
                totalDiff += (ip - cp.Img).Length();
                Assert.IsTrue((ip - cp.Img).Length() < 1.0,
                    "Point after linear estimation too far : " + (ip - cp.Img).Length());
            }

            _calib.HomoPoints();
             _calib.NormalizeImagePoints();
             _calib.NormalizeRealPoints();
            _calib.CameraMatrix = _calib.FindLinearEstimationOfCameraMatrix();
            _calib._pointsNormalised = true;
            _calib.FindNormalisedVariances();
            _calib.UseCovarianceMatrix = true;
            var lecm = _eCM.Clone();

            // Disturb camera matrix a little
              //          _calib.CameraMatrix = AddNoise(_calib.CameraMatrix);

            _calib._miniAlg.DoComputeJacobianNumerically = true;
            _calib._miniAlg.NumericalDerivativeStep = 1e-4;
            _calib.MinimizeError();

            // _calib.DenormaliseCameraMatrix();
            _calib.DecomposeCameraMatrix();
            _eCM = _calib.CameraMatrix;

            var errVec = _eCM.PointwiseDivide_NoNaN(lecm);
            double err = errVec.L2Norm();

            double scaleK = 1.0 / _calib.CameraInternalMatrix[2, 2];
            _eCM.MultiplyThis(-scaleK);

            var eK = _calib.CameraInternalMatrix.Multiply(scaleK);
            var eR = -_calib.CameraRotationMatrix;
            var eC = -(_eCM.SubMatrix(0, 3, 0, 3).Inverse() * _eCM.Column(3));

            Matrix<double> eExt = new DenseMatrix(3, 4);
            eExt.SetSubMatrix(0, 0, eR);
            eExt.SetColumn(3, -eR * eC);

            var eCM = eK * eExt;

            // var errVec = _CM.PointwiseDivide_NoNaN(_eCM);
            // double err = errVec.L2Norm();
            // Assert.IsTrue(
            //    Math.Abs(err - Math.Sqrt(12)) < Math.Sqrt(12) / 1000.0 || // max 0.1% diffrence
            //    (_eCM - _CM).FrobeniusNorm() < 1e-3);

            double estDiff = 0;
            for(int p = 0; p < _pointsCount; ++p)
            {
                var cp = points[p];
                Vector<double> rp = new DenseVector(4);
                rp[0] = cp.RealX;
                rp[1] = cp.RealY;
                rp[2] = cp.RealZ;
                rp[3] = 1.0;
                var imagePoint = _eCM * rp;

                Vector2 ip = new Vector2(imagePoint[0] / imagePoint[2], imagePoint[1] / imagePoint[2]);
                estDiff += (ip - cp.Img).Length();
                Assert.IsTrue((ip - cp.Img).Length() < 1.5,
                        "Point after error minimalisation too far : " + (ip - cp.Img).Length());
            }

            var minialg = _calib._miniAlg;
            // Test conovergence :
            // ||mX-rX|| = ||mX-eX|| + ||rX-eX|| (or squared??)
            // rX - real point from 'points'
            // mX - measured point, noised
            // eX = estimated X from result vector for 3d points and ePeX for image point
            double len2_mr = 0;
//.........这里部分代码省略.........

        public void op_Column_intNegative()
        {
            var matrix = new Matrix<string>();

            Assert.Throws<ArgumentOutOfRangeException>(() => matrix.Column(-1).ToList());
        }

        public void op_Column_int()
        {
            var matrix = new Matrix<string>(3, 3);
            matrix[0, 0] = "0,0";
            matrix[0, 1] = "0,1";
            matrix[0, 2] = "0,2";
            matrix[1, 0] = "1,0";
            matrix[1, 1] = "1,1";
            matrix[1, 2] = "1,2";
            matrix[2, 0] = "2,0";
            matrix[2, 1] = "2,1";
            matrix[2, 2] = "2,2";

            for (var x = 0; x < 3; x++)
            {
                var y = 0;
                foreach (var actual in matrix.Column(x))
                {
                    var expected = "{0},{1}".FormatWith(x, y);

                    Assert.Equal(expected, actual);
                    y++;
                }
            }
        }

        public static void PCA(DenseMatrix input, out Vector<double> latent, out Matrix<double> score, out Matrix<double> coeff)
        {
            int n = input.RowCount;
            int p = input.ColumnCount;

            //de-mean input
            var tmpInput = DenseMatrix.OfMatrix(input);
            for (int i = 0; i < tmpInput.ColumnCount; i++)
            {
                double avg = tmpInput.Column(i).Average();
                tmpInput.SetColumn(i, tmpInput.Column(i).Subtract(avg));
            }

            var svd = tmpInput.Svd(true);
            var sigma = svd.S;
            var tmpCoeff = svd.VT.Transpose();

            score = DenseMatrix.Create(n, p, (_, __) => 0);
            var U = svd.U.SubMatrix(0, n, 0, p);
            for (int i = 0; i < U.RowCount; i++)
            {
                score.SetRow(i, U.Row(i).PointwiseMultiply(sigma));
            }

            sigma = sigma.Divide(Math.Sqrt(n - 1));
            latent = sigma.PointwiseMultiply(sigma);

            //give the largest absolute value in each column a positive sign
            var maxIndices = tmpCoeff.EnumerateColumns().Select(x => x.AbsoluteMaximumIndex());
            var colSigns = maxIndices.Select((x, j) => Math.Sign(tmpCoeff[x, j])).ToList();
            for (int j = 0; j < tmpCoeff.ColumnCount; j++)
            {
                tmpCoeff.SetColumn(j, tmpCoeff.Column(j) * colSigns[j]);
                score.SetColumn(j, score.Column(j) * colSigns[j]);
            }

            coeff = tmpCoeff;
        }

        public void Test_EstimateCameraMatrix_NoisedLinear()
        {
            PrepareCameraMatrix();
            var points = GenerateCalibrationPoints_Random();
            PrepareCalibrator(
                AddNoise(points, _varianceReal, _varianceImage));

            _calib.HomoPoints();

            _calib.NormalizeImagePoints();
            _calib.NormalizeRealPoints();

            _calib.CameraMatrix = _calib.FindLinearEstimationOfCameraMatrix();

            _calib.DenormaliseCameraMatrix();
            _calib.DecomposeCameraMatrix();
            _eCM = _calib.CameraMatrix;

            double scaleK = 1.0 / _calib.CameraInternalMatrix[2, 2];
            _eCM.MultiplyThis(-scaleK);

            var eK = _calib.CameraInternalMatrix.Multiply(scaleK);
            var eR = -_calib.CameraRotationMatrix;
            var eC = -(_eCM.SubMatrix(0, 3, 0, 3).Inverse() * _eCM.Column(3));

            Matrix<double> eExt = new DenseMatrix(3, 4);
            eExt.SetSubMatrix(0, 0, eR);
            eExt.SetColumn(3, -eR * eC);

            var eCM = eK * eExt;

            var errVec = _CM.PointwiseDivide_NoNaN(_eCM);
            double err = errVec.L2Norm();
            Assert.IsTrue(
                Math.Abs(err - Math.Sqrt(12)) < Math.Sqrt(12) / 50.0 || // max 2% diffrence
                (_eCM - _CM).FrobeniusNorm() < 1e-3);

            for(int p = 0; p < _pointsCount; ++p)
            {
                var cp = points[p];
                Vector<double> rp = new DenseVector(4);
                rp[0] = cp.RealX;
                rp[1] = cp.RealY;
                rp[2] = cp.RealZ;
                rp[3] = 1.0;
                var imagePoint = _eCM * rp;

                Vector2 ip = new Vector2(imagePoint[0] / imagePoint[2], imagePoint[1] / imagePoint[2]);
                Assert.IsTrue((ip - cp.Img).Length() < cp.Img.Length() / 10.0);
            }
        }

        /// <summary>
        /// Multiplies this matrix with another matrix and places the results into the result matrix.
        /// </summary>
        /// <param name="other">The matrix to multiply with.</param>
        /// <param name="result">The result of the multiplication.</param>
        protected override void DoMultiply(Matrix<float> other, Matrix<float> result)
        {
            result.Clear();
            var columnVector = new DenseVector(other.RowCount);

            var rowPointers = _storage.RowPointers;
            var columnIndices = _storage.ColumnIndices;
            var values = _storage.Values;

            for (var row = 0; row < RowCount; row++)
            {
                var startIndex = rowPointers[row];
                var endIndex = rowPointers[row + 1];

                if (startIndex == endIndex)
                {
                    continue;
                }

                for (var column = 0; column < other.ColumnCount; column++)
                {
                    // Multiply row of matrix A on column of matrix B
                    other.Column(column, columnVector);

                    var sum = 0f;
                    for (var index = startIndex; index < endIndex; index++)
                    {
                        sum += values[index] * columnVector[columnIndices[index]];
                    }

                    result.At(row, column, sum);
                }
            }
        }