C# GeneralMatrix.Transpose方法代码示例

        /**
        * Computes the Moore–Penrose pseudoinverse using the SVD method.
        *
        * Modified version of the original implementation by Kim van der Linde.
        */
        public static GeneralMatrix pinv(GeneralMatrix x)
        {
            if (x.Rank() < 1)
                return null;

            if (x.ColumnDimension > x.RowDimension)
                return pinv(x.Transpose()).Transpose();

            SingularValueDecomposition svdX = new SingularValueDecomposition(x);
            double[] singularValues = svdX.SingularValues;
            double tol = Math.Max(x.ColumnDimension, x.RowDimension)
                    * singularValues[0] * 2E-16;

            double[] singularValueReciprocals = new double[singularValues.Count()];
            for (int i = 0; i < singularValues.Count(); i++)
                singularValueReciprocals[i] = Math.Abs(singularValues[i]) < tol ? 0
                        : (1.0 / singularValues[i]);

            double[][] u = svdX.GetU().Array;
            double[][] v = svdX.GetV().Array;

            int min = Math.Min(x.ColumnDimension, u[0].Count());

            double[][] inverse = new double[x.ColumnDimension][];

            for (int i = 0; i < x.ColumnDimension; i++) {
                inverse[i] = new double[x.RowDimension];

                for (int j = 0; j < u.Count(); j++)
                    for (int k = 0; k < min; k++)
                        inverse[i][j] += v[i][k] * singularValueReciprocals[k] * u[j][k];
            }
            return new GeneralMatrix(inverse);
        }

        protected override GeneralMatrix CalculateNextHessianApproximation(GeneralMatrix previousH, 
			double[]prevX, double[]curX, double[]prevGrad, double[]curGrad)
        {
            GeneralMatrix currentH = new GeneralMatrix(_nDim,_nDim);
            GeneralMatrix cX = new GeneralMatrix(curX,_nDim);
            GeneralMatrix pX = new GeneralMatrix(prevX,_nDim);
            GeneralMatrix cG = new GeneralMatrix(curGrad,_nDim);
            GeneralMatrix pG = new GeneralMatrix(prevGrad,_nDim);

            GeneralMatrix dX = cX.Subtract(pX);
            GeneralMatrix dG = cG.Subtract(pG);

            double aK1 = 1/(dX.Transpose().Multiply(dG).GetElement(0,0));
            GeneralMatrix aK2 = dX.Multiply(dX.Transpose());

            GeneralMatrix aK = aK2.Multiply(aK1);

            double bK1 = -1/(dG.Transpose().Multiply(previousH).Multiply(dG).GetElement(0,0));
            GeneralMatrix bK2 = previousH.Multiply(dG).Multiply(dG.Transpose()).Multiply(previousH.Transpose());

            GeneralMatrix bK =bK2.Multiply(bK1);

            currentH = previousH.Add(aK).Add(bK);

            return currentH;
        }

 /// <summary>Solve X*A = B, which is also A'*X' = B'</summary>
 /// <param name="B">   right hand side
 /// </param>
 /// <returns>     solution if A is square, least squares solution otherwise.
 /// </returns>
 public virtual GeneralMatrix SolveTranspose(GeneralMatrix B)
 {
     return Transpose().Solve(B.Transpose());
 }

//.........这里部分代码省略.........
                try_success("ArrayMultiplyEquals... ", "");
            }
            catch (System.SystemException e)
            {
                errorCount = try_failure(errorCount, "ArrayMultiplyEquals... ", "(A = R, A = A.*B, but A./B != R)");
                System.Console.Out.WriteLine(e.Message);
            }

            /// <summary>LA methods:
            /// Transpose
            /// Multiply
            /// Condition
            /// Rank
            /// Determinant
            /// trace
            /// Norm1
            /// norm2
            /// normF
            /// normInf
            /// Solve
            /// solveTranspose
            /// Inverse
            /// chol
            /// Eigen
            /// lu
            /// qr
            /// svd
            ///
            /// </summary>

            print("\nTesting linear algebra methods...\n");
            A = new GeneralMatrix(columnwise, 3);
            T = new GeneralMatrix(tvals);
            T = A.Transpose();
            try
            {
                check(A.Transpose(), T);
                try_success("Transpose...", "");
            }
            catch (System.SystemException e)
            {
                errorCount = try_failure(errorCount, "Transpose()...", "Transpose unsuccessful");
                System.Console.Out.WriteLine(e.Message);
            }
            A.Transpose();
            try
            {
                check(A.Norm1(), columnsummax);
                try_success("Norm1...", "");
            }
            catch (System.SystemException e)
            {
                errorCount = try_failure(errorCount, "Norm1()...", "incorrect norm calculation");
                System.Console.Out.WriteLine(e.Message);
            }
            try
            {
                check(A.NormInf(), rowsummax);
                try_success("normInf()...", "");
            }
            catch (System.SystemException e)
            {
                errorCount = try_failure(errorCount, "normInf()...", "incorrect norm calculation");
                System.Console.Out.WriteLine(e.Message);
            }
            try

        public void TestTranspose()
        {
            GeneralMatrix _gm = new GeneralMatrix(2,2);
            _gm.SetElement(0,0,1);
            _gm.SetElement(0,1,2);
            _gm.SetElement(1,0,3);
            _gm.SetElement(1,1,4);
            GeneralMatrix _ngm = _gm.Transpose();

            Assert.AreEqual(1,0,2);
        }

        //encrypt line of input
        static void Encrypt(string input, GeneralMatrix key, Dictionary<char, int> hillDict, StringBuilder encText)
        {
            string ct = "";
            string[] pt = input.Trim().Split(' ');
            for (int i = 0; i < pt.Length; i += 3)
            {
                //encrypt 3 letters at a time
                double[] temp = new double[3];
                temp[0] = Convert.ToInt32(pt[i]);
                temp[1] = Convert.ToInt32(pt[i + 1]);
                temp[2] = Convert.ToInt32(pt[i + 2]);

                //create plain text matrix, transpose and encrypt it
                GeneralMatrix ptMat = new GeneralMatrix(new double[] { temp[0], temp[1], temp[2] }, 3);
                GeneralMatrix trnasPTMat = ptMat.Transpose();
                GeneralMatrix ctMat = key.Multiply(ptMat);
                for (int x = 0; x < ctMat.RowDimension; x++)
                {
                    for (int y = 0; y < ctMat.ColumnDimension; y++)
                    {
                        var tempElement = Convert.ToInt32(ctMat.GetElement(x, y)) % 31;
                        ctMat.SetElement(x, y, tempElement);
                    }
                }

                for (int x = 0; x < ctMat.RowDimension; x++)
                {
                    for (int y = 0; y < ctMat.ColumnDimension; y++)
                    {
                        ct += hillDict.FirstOrDefault(z => z.Value == ctMat.GetElement(x, y)).Key;
                    }
                }
            }
            //append to string builder
            encText.AppendLine(ct);
        }

        /// <summary>
        /// (
        /// </summary>
        /// <param name="matrix"></param>
        /// <returns></returns>
        public static GeneralMatrix ExpandUtility(GeneralMatrix matrix)
        {
            double val=0.0;
            int n = matrix.RowDimension;
            int m = matrix.ColumnDimension;

            if (n!=m)
                throw new ArgumentException("Criteria matrix must be symmetrical");

            GeneralMatrix newMatrix = matrix.Transpose();

            //for all transposed elements calculate their inverse values
            //set diagonal elements to 0
            for (int i=0; i<n; i++)
                for (int j=0; j<=i; j++)
                {
                    val = newMatrix.GetElement(i,j);
                    if (val==0.0)
                        throw new ArgumentException("Criteria comparison values van't be 0");
                        newMatrix.SetElement(i,j,1/val);
                    if (i==j)
                        newMatrix.SetElement(i,j,0);
                }
            //add transposed, inverse matrix to the original one
            //create fully expanded matrix
            return newMatrix.Add(matrix);
        }

		// this is the straight forward implementation of the kabsch algorithm.
		// see http://en.wikipedia.org/wiki/Kabsch_algorithm for a detailed explanation.
		public void Evaluate(int SpreadMax)
		{
			FOutput.SliceCount = 1;
			Matrix4x4 mOut;
			
			if (FInputP.SliceCount > 0 && FInputQ.SliceCount > 0 && FInputEnabled[0]){

				// ======================== STEP 1 ========================
				// translate both sets so that their centroids coincides with the origin
				// of the coordinate system.
				
				double[] meanP = new double[3]{0.0,0.0,0.0}; // mean of first point set
				for(int i=0; i < FInputP.SliceCount;i++){
					meanP[0] += FInputP[i].x;
					meanP[1] += FInputP[i].y;
					meanP[2] += FInputP[i].z;
				}
				meanP[0] /= FInputP.SliceCount;
				meanP[1] /= FInputP.SliceCount;
				meanP[2] /= FInputP.SliceCount;
				double[][]centroidP = new double[3][]{new double[]{meanP[0]},new double[]{meanP[1]},new double[]{meanP[2]}};
				GeneralMatrix mCentroidP = new GeneralMatrix(centroidP);
				
				double[][] arrayP = new double[FInputP.SliceCount][];
				for(int i=0; i < FInputP.SliceCount;i++){
					arrayP[i] = new double[3];
					arrayP[i][0] = FInputP[i].x - meanP[0]; // subtract the mean values from the incoming pointset
					arrayP[i][1] = FInputP[i].y - meanP[1];
					arrayP[i][2] = FInputP[i].z - meanP[2];
				}
				// this is the matrix of the first pointset translated to the origin of the coordinate system
				GeneralMatrix P = new GeneralMatrix(arrayP); 
	
				double[] meanQ = new double[3]{0.0,0.0,0.0}; // mean of second point set
				for(int i=0; i < FInputQ.SliceCount;i++){
					meanQ[0] += FInputQ[i].x;
					meanQ[1] += FInputQ[i].y;
					meanQ[2] += FInputQ[i].z;
				}
				meanQ[0] /= FInputQ.SliceCount;
				meanQ[1] /= FInputQ.SliceCount;
				meanQ[2] /= FInputQ.SliceCount;
				double[][]centroidQ = new double[3][]{new double[]{meanQ[0]},new double[]{meanQ[1]},new double[]{meanQ[2]}};
				GeneralMatrix mCentroidQ = new GeneralMatrix(centroidQ);
				
				double[][] arrayQ = new double[FInputQ.SliceCount][];
				for(int i=0; i < FInputQ.SliceCount;i++){
					arrayQ[i] = new double[3];
					arrayQ[i][0] = FInputQ[i].x - meanQ[0]; // subtract the mean values from the incoming pointset
					arrayQ[i][1] = FInputQ[i].y - meanQ[1];
					arrayQ[i][2] = FInputQ[i].z - meanQ[2];
				}
				// this is the matrix of the second pointset translated to the origin of the coordinate system
				GeneralMatrix Q = new GeneralMatrix(arrayQ);
				
				
				// ======================== STEP2 ========================
				// calculate a covariance matrix A and compute the optimal rotation matrix
				GeneralMatrix A = P.Transpose() * Q;
				
				SingularValueDecomposition svd = A.SVD();
				GeneralMatrix U = svd.GetU();
				GeneralMatrix V = svd.GetV();
				
				// calculate determinant for a special reflexion case.
				double det = (V * U.Transpose()).Determinant();
				double[][] arrayD = new double[3][]{ new double[]{1,0,0},
													new double[] {0,1,0},
													new double[] {0,0,1}
				};		
				arrayD[2][2] = det < 0 ? -1 : 1; // multiply 3rd column with -1 if determinant is < 0
				GeneralMatrix D = new GeneralMatrix(arrayD);
				
				// now we can compute the rotation matrix:
				GeneralMatrix R = V * D * U.Transpose();
				
				// ======================== STEP3 ========================
				// calculate the translation:
				GeneralMatrix T = mCentroidP - R.Inverse() * mCentroidQ;
				
				// ================== OUTPUT TRANSFORM ===================
				
				mOut.m11 = (R.Array)[0][0];
				mOut.m12 = (R.Array)[0][1];
				mOut.m13 = (R.Array)[0][2];
				mOut.m14 = 0;
				mOut.m21 = (R.Array)[1][0];
				mOut.m22 = (R.Array)[1][1];
				mOut.m23 = (R.Array)[1][2];
				mOut.m24 = 0;			
				mOut.m31 = (R.Array)[2][0];
				mOut.m32 = (R.Array)[2][1];
				mOut.m33 = (R.Array)[2][2];
				mOut.m34 = 0;
				mOut.m41 = (T.Array)[0][0];
				mOut.m42 = (T.Array)[1][0];
				mOut.m43 = (T.Array)[2][0];
				mOut.m44 = 1;
//.........这里部分代码省略.........

        /// <summary>
        /// Solves between two point sets
        /// </summary>
        /// <param name="points">Point set</param>
        /// <returns>Affine matrix for each point set</returns>
        public Matrix Solve(IReadOnlyList<CameraToCameraPoint> points)
        {
            if (points == null)
                throw new ArgumentNullException("points");
            if (points.Count < 1)
                throw new ArgumentException("points", "No points provided");

            double[] meanP = new double[3] { 0.0, 0.0, 0.0 }; // mean of first point set
            for (int i = 0; i < points.Count; i++)
            {
                Vector3 orig = points[i].Origin;
                meanP[0] += orig.X;
                meanP[1] += orig.Y;
                meanP[2] += orig.Z;
            }

            double invCount = 1.0 / (double)points.Count;

            meanP[0] *= invCount;
            meanP[1] *= invCount;
            meanP[2] *= invCount;
            double[][] centroidP = new double[3][] { new double[] { meanP[0] }, new double[] { meanP[1] }, new double[] { meanP[2] } };
            GeneralMatrix mCentroidP = new GeneralMatrix(centroidP);

            double[][] arrayP = new double[points.Count][];
            for (int i = 0; i < points.Count; i++)
            {
                Vector3 orig = points[i].Origin;
                arrayP[i] = new double[3];
                arrayP[i][0] = orig.X - meanP[0]; // subtract the mean values from the incoming pointset
                arrayP[i][1] = orig.Y - meanP[1];
                arrayP[i][2] = orig.Z - meanP[2];
            }
            // this is the matrix of the first pointset translated to the origin of the coordinate system
            GeneralMatrix P = new GeneralMatrix(arrayP);

            double[] meanQ = new double[3] { 0.0, 0.0, 0.0 }; // mean of second point set
            for (int i = 0; i < points.Count; i++)
            {
                Vector3 dest = points[i].Destination;
                meanQ[0] += dest.X;
                meanQ[1] += dest.Y;
                meanQ[2] += dest.Z;
            }
            meanQ[0] *= invCount;
            meanQ[1] *= invCount;
            meanQ[2] *= invCount;
            double[][] centroidQ = new double[3][] { new double[] { meanQ[0] }, new double[] { meanQ[1] }, new double[] { meanQ[2] } };
            GeneralMatrix mCentroidQ = new GeneralMatrix(centroidQ);

            double[][] arrayQ = new double[points.Count][];
            for (int i = 0; i < points.Count; i++)
            {
                Vector3 dest = points[i].Destination;
                arrayQ[i] = new double[3];
                arrayQ[i][0] = dest.X - meanQ[0]; // subtract the mean values from the incoming pointset
                arrayQ[i][1] = dest.Y - meanQ[1];
                arrayQ[i][2] = dest.Z - meanQ[2];
            }
            // this is the matrix of the second pointset translated to the origin of the coordinate system
            GeneralMatrix Q = new GeneralMatrix(arrayQ);

            // ======================== STEP2 ========================
            // calculate a covariance matrix A and compute the optimal rotation matrix
            GeneralMatrix A = P.Transpose() * Q;

            SingularValueDecomposition svd = A.SVD();
            GeneralMatrix U = svd.GetU();
            GeneralMatrix V = svd.GetV();

            // calculate determinant for a special reflexion case.
            double det = (V * U.Transpose()).Determinant();
            double[][] arrayD = new double[3][]{ new double[]{1,0,0},
                                                    new double[] {0,1,0},
                                                    new double[] {0,0,1}
                };
            arrayD[2][2] = det < 0 ? -1 : 1; // multiply 3rd column with -1 if determinant is < 0
            GeneralMatrix D = new GeneralMatrix(arrayD);

            // now we can compute the rotation matrix:
            GeneralMatrix R = V * D * U.Transpose();

            // ======================== STEP3 ========================
            // calculate the translation:
            GeneralMatrix T = mCentroidP - R.Inverse() * mCentroidQ;

            return new Matrix((float)(R.Array)[0][0],
               (float)(R.Array)[0][1],
               (float)(R.Array)[0][2],
                0.0f,
               (float)(R.Array)[1][0],
               (float)(R.Array)[1][1],
               (float)(R.Array)[1][2],
               0.0f,
               (float)(R.Array)[2][0],
//.........这里部分代码省略.........