C# IMatrix.Any方法代码示例

		/// <summary>
		/// Execution of the fast nonnegative least squares algorithm. The algorithm finds a vector x with all elements xi&gt;=0 which minimizes |X*x-y|.
		/// </summary>
		/// <param name="XtX">X transposed multiplied by X, thus a square matrix.</param>
		/// <param name="Xty">X transposed multiplied by Y, thus a matrix with one column and same number of rows as X.</param>
		/// <param name="isRestrictedToPositiveValues">Function that takes the parameter index as argument and returns true if the parameter at this index is restricted to positive values; otherwise the return value must be false.</param>
		/// <param name="tolerance">Used to decide if a solution element is less than or equal to zero. If this is null, a default tolerance of tolerance = MAX(SIZE(XtX)) * NORM(XtX,1) * EPS is used.</param>
		/// <param name="x">Output: solution vector (matrix with one column and number of rows according to dimension of X.</param>
		/// <param name="w">Output: Lagrange vector. Elements which take place in the fit are set to 0. Elements fixed to zero contain a negative number.</param>
		/// <remarks>
		/// <para>
		/// Literature: Rasmus Bro and Sijmen De Jong, 'A fast non-negativity-constrained least squares algorithm', Journal of Chemometrics, Vol. 11, 393-401 (1997)
		/// </para>
		/// <para>
		/// Algorithm modified by Dirk Lellinger 2015 to allow a mixture of restricted and unrestricted parameters.
		/// </para>
		/// </remarks>
		public static void Execution(IROMatrix XtX, IROMatrix Xty, Func<int, bool> isRestrictedToPositiveValues, double? tolerance, out IMatrix x, out IMatrix w)
		{
			if (null == XtX)
				throw new ArgumentNullException(nameof(XtX));
			if (null == Xty)
				throw new ArgumentNullException(nameof(Xty));
			if (null == isRestrictedToPositiveValues)
				throw new ArgumentNullException(nameof(isRestrictedToPositiveValues));

			if (XtX.Rows != XtX.Columns)
				throw new ArgumentException("Matrix should be a square matrix", nameof(XtX));
			if (Xty.Columns != 1)
				throw new ArgumentException(nameof(Xty) + " should be a column vector (number of columns should be equal to 1)", nameof(Xty));
			if (Xty.Rows != XtX.Columns)
				throw new ArgumentException("Number of rows in " + nameof(Xty) + " should match number of columns in " + nameof(XtX), nameof(Xty));

			var matrixGenerator = new Func<int, int, DoubleMatrix>((rows, cols) => new DoubleMatrix(rows, cols));

			// if nargin < 3
			//   tol = 10 * eps * norm(XtX, 1) * length(XtX);
			// end
			double tol = tolerance.HasValue ? tolerance.Value : 10 * DoubleConstants.DBL_EPSILON * MatrixMath.Norm(XtX, MatrixNorm.M1Norm) * Math.Max(XtX.Rows, XtX.Columns);

			//	[m, n] = size(XtX);
			int n = XtX.Columns;

			// P = zeros(1, n);
			// Z = 1:n;
			var P = new bool[n]; // POSITIVE SET: all indices which are currently not fixed are marked with TRUE (Negative set is simply this, but inverted)
			bool initializationOfSolutionRequired = false;
			for (int i = 0; i < n; ++i)
			{
				bool isNotRestricted = !isRestrictedToPositiveValues(i);
				P[i] = isNotRestricted;
				initializationOfSolutionRequired |= isNotRestricted;
			}

			// x = P';
			x = matrixGenerator(n, 1);

			// w = Xty-XtX*x;
			w = matrixGenerator(n, 1);
			MatrixMath.Copy(Xty, w);
			var helper_n_1 = matrixGenerator(n, 1);
			MatrixMath.Multiply(XtX, x, helper_n_1);
			MatrixMath.Subtract(w, helper_n_1, w);

			// set up iteration criterion
			int iter = 0;
			int itmax = 30 * n;

			// outer loop to put variables into set to hold positive coefficients
			// while any(Z) & any(w(ZZ) > tol)
			while (initializationOfSolutionRequired || (P.Any(ele => false == ele) && w.Any((r, c, ele) => false == P[r] && ele > tol)))
			{
				if (initializationOfSolutionRequired)
				{
					initializationOfSolutionRequired = false;
				}
				else
				{
					// [wt, t] = max(w(ZZ));
					// t = ZZ(t);
					int t = -1; // INDEX
					double wt = double.NegativeInfinity;
					for (int i = 0; i < n; ++i)
					{
						if (!P[i])
						{
							if (w[i, 0] > wt)
							{
								wt = w[i, 0];
								t = i;
							}
						}
					}

					// P(1, t) = t;
					// Z(t) = 0;
					P[t] = true;
				}

				// z(PP')=(Xty(PP)'/XtX(PP,PP)');
//.........这里部分代码省略.........