13,861,366 members

#### Stats

83.8K views
92 bookmarked
Posted 19 Sep 2012
Licenced CPOL

# Yet Another Math Parser (YAMP)

, 30 Sep 2012
Constructing a fast math parser using Reflection to do numerics like Matlab.
 ```using System; using YAMP; namespace YAMP.Numerics { /// /// QR Decomposition. /// For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n /// orthogonal matrix Q and an n-by-n upper triangular matrix R so that /// A = Q*R. /// The QR decompostion always exists, even if the matrix does not have /// full rank, so the constructor will never fail. The primary use of the /// QR decomposition is in the least squares solution of nonsquare systems /// of simultaneous linear equations. This will fail if IsFullRank() /// returns false. /// public class QRDecomposition : DirectSolver { #region Members /// /// Array for internal storage of decomposition. /// double[][] QR; /// /// Row and column dimensions. /// int m, n; /// /// Array for internal storage of diagonal of R. /// double[] Rdiag; #endregion // Class variables #region Constructor /// /// QR Decomposition, computed by Householder reflections. /// /// Rectangular matrix /// Structure to access R and the Householder vectors and compute Q. public QRDecomposition(MatrixValue A) { // Initialize. QR = A.GetRealArray(); m = A.DimensionY; n = A.DimensionX; Rdiag = new double[n]; // Main loop. for (int k = 0; k < n; k++) { // Compute 2-norm of k-th column without under/overflow. var nrm = 0.0; for (int i = k; i < m; i++) nrm = NumericHelpers.Hypot(nrm, QR[i][k]); if (nrm != 0.0) { // Form k-th Householder vector. if (QR[k][k] < 0) nrm = -nrm; for (int i = k; i < m; i++) QR[i][k] /= nrm; QR[k][k] += 1.0; // Apply transformation to remaining columns. for (int j = k + 1; j < n; j++) { var s = 0.0; for (int i = k; i < m; i++) s += QR[i][k] * QR[i][j]; s = (-s) / QR[k][k]; for (int i = k; i < m; i++) QR[i][j] += s * QR[i][k]; } } Rdiag[k] = -nrm; } } #endregion // Constructor #region Public Properties /// /// Is the matrix full rank? /// /// true if R, and hence A, has full rank. virtual public bool FullRank { get { for (int j = 0; j < n; j++) { if (Rdiag[j] == 0) return false; } return true; } } /// /// Return the Householder vectors /// /// Lower trapezoidal matrix whose columns define the reflections. virtual public MatrixValue H { get { var X = new MatrixValue(m, n); for (int i = 1; i <= m; i++) { for (int j = 1; j <= n; j++) { if (i >= j) X[i, j].Value = QR[i - 1][j - 1]; else X[i, j].Value = 0.0; } } return X; } } /// /// Return the upper triangular factor /// /// R virtual public MatrixValue R { get { var X = new MatrixValue(n, n); for (int i = 1; i <= n; i++) { for (int j = 1; j <= n; j++) { if (i < j) X[i, j].Value = QR[i - 1][j - 1]; else if (i == j) X[i, j].Value = Rdiag[i - 1]; else X[i, j].Value = 0.0; } } return X; } } /// /// Generate and return the (economy-sized) orthogonal factor /// /// Q virtual public MatrixValue Q { get { var X = new MatrixValue(m, n); for (int k = n; k > 0; k--) { for (int i = 1; i <= m; i++) Q[i, k].Value = 0.0; Q[k, k].Value = 1.0; for (int j = k; j <= n; j++) { var l = k - 1; if (QR[l][l] != 0) { var s = 0.0; for (int i = k; i <= m; i++) s += QR[i - 1][l] * Q[i, j].Value; s = (-s) / QR[l][l]; for (int i = k; i <= m; i++) Q[i, j].Value += s * QR[i - 1][l]; } } } return X; } } #endregion // Public Properties #region Public Methods /// /// Least squares solution of A*X = B /// /// A Matrix with as many rows as A and any number of columns. /// X that minimizes the two norm of Q*R*X-B. /// Matrix row dimensions must agree. /// Matrix is rank deficient. public override MatrixValue Solve(MatrixValue B) { if (B.DimensionY != m) throw new DimensionException(B.DimensionY, m); if (!this.FullRank) throw new MatrixFormatException("full rank"); // Copy right hand side var nx = B.DimensionX; var X = B.GetRealArray(); // Compute Y = transpose(Q)*B for (int k = 0; k < n; k++) { for (int j = 0; j < nx; j++) { var s = 0.0; for (int i = k; i < m; i++) s += QR[i][k] * X[i][j]; s = (-s) / QR[k][k]; for (int i = k; i < m; i++) X[i][j] += s * QR[i][k]; } } // Solve R*X = Y; for (int k = n - 1; k >= 0; k--) { for (int j = 0; j < nx; j++) X[k][j] /= Rdiag[k]; for (int i = 0; i < k; i++) { for (int j = 0; j < nx; j++) X[i][j] -= X[k][j] * QR[i][k]; } } return new MatrixValue(X, n, nx).SubMatrix(0, n, 0, nx); } #endregion // Public Methods } }```

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

## Share

 Architect Germany
Florian lives in Munich, Germany. He started his programming career with Perl. After programming C/C++ for some years he discovered his favorite programming language C#. He did work at Siemens as a programmer until he decided to study Physics.

During his studies he worked as an IT consultant for various companies. After graduating with a PhD in theoretical particle Physics he is working as a senior technical consultant in the field of home automation and IoT.

Florian has been giving lectures in C#, HTML5 with CSS3 and JavaScript, software design, and other topics. He is regularly giving talks at user groups, conferences, and companies. He is actively contributing to open-source projects. Florian is the maintainer of AngleSharp, a completely managed browser engine.