in

.NET Opensource

Community for opensource projects by Christoph Rüegg

Rewrite of Matrix Class

Last post 05-06-2007 20:35 by Christoph Rüegg. 6 replies.
Page 1 of 1 (7 items)
Sort Posts: Previous Next
  • 05-04-2007 13:21

    • MovGP0
    • Top 10 Contributor
      Male
    • Joined on 11-21-2005
    • Austria
    • Posts 49

    Rewrite of Matrix Class

    I'm started to write a new Matrix Class that is more powerful than the current one.

    This is just a try to see how a matrix class that could get implemented.

    Design Issues:

    • There is a separate partial class that implements general interfaces to make this class as useful as possible.
    • The Class is (currently) not optimized for performance, but for functionality
    • There is a separate partial class that implements calculations for determinants. The best algorithm should get determined by checking the type and form of the matrix. ie. for 2x2 and 3x3 Matrices the rule of sarrus is most effective and used automatically.
    • A separate partial class should check the form of the matrix (triangular, or similar)
    • The class is written Generic, so that it can handle different datatypes. Special datatypes like complex numbers might need subclasses that implement special functions.
    • more isses will follow as the class grows
     Please take a quick look at this classes and say me your meaning, so that I can correct issues in early developoment stage. Thanks!
  • 05-04-2007 13:22 In reply to

    • MovGP0
    • Top 10 Contributor
      Male
    • Joined on 11-21-2005
    • Austria
    • Posts 49

    Matrix Class

    /*
     * Generic Matrix Class for Math.NET
     * written by MovGP0 (mailto:movgp0@gmail.com)
     * licensed under the terms of the CC GNU LGPL 2.1
     * (http://creativecommons.org/licenses/LGPL/2.1/)
     */
    using System;
    using System.Collections.Generic;
    using System.Linq;
    using System.Text;

    namespace MathNet.Numerics.LinearAlgebra
    {
        /// <summary>This is a generic Matrix Class.</summary>
        /// <typeparam name="T">
        /// Type of the Matrix Elwements.
        /// Usually the type should be real, double, or rational.
        /// Note: There should be an extended Matrix-Class for complex values.
        /// </typeparam>
        partial class NMatrix<T>
        {
            #region Array
            private T[,] array;
            private uint m
            {
                get{
                    return this.array.GetUpperBound(0);
                }
            }
            private uint n
            {
                get
                {
                    return this.array.GetUpperBound(1);
                }
            }
            #endregion
            #region Constructors

            /// <summary>
            /// Creates a new square Matrix.
            /// </summary>
            /// <param name="Dimension">The Dimension N of the NxN Matrix.</param>
            public NMatrix(uint Dimension)
            {
                NMatrix(Dimension, Dimension);
            }

            /// <summary>Creates a new Matrix from a 2-dimensional Array</summary>
            /// <param name="Array">
            /// The Array where the Matrix is built from.
            /// The array must be of rectangular shape and has
            /// to be one or two dimensional.
            /// </param>
            public NMatrix(T [,] Array )
            {//TODO: check C# Documentation if Array is forced to be 2D because of the "[,]"

                // TODO: remove if Array is forced to be 2D
                // Check if the Array is 1D or 2D
                if (Array.Rank < 1 || Array.Rank > 2)
                {
                    string message = String.Format(
                        "Array must be of dimension 1 or 2.{0}" +
                        "Dimension of the Array is {1}.",
                        Environment.NewLine, Array.Rank);
                    throw new ArgumentException(message);
                }

                #region copy the elements
                uint x,y = 1;

                // create a new array of fitting size
                uint[] lowerBonds = { 1, 1 }; // index starting at 1,1 instead of 0,0
                uint[] lengths = { // size of the array
                    // size of first dimension (lenght of the vector)
                    Array.GetLength(0),

                    // if Array is 1D then the new matrix is a vector (1 column)
                    //TODO: simplify to "Array.GetLength(1)" if Array is forced to be 2D
                    (Array.Rank == 2) ? Array.GetLength(1) : 1
                };
                array = Array.CreateInstance(typeof(T), lengths, lowerBonds);

                for(int i = Array.GetLowerBound(0); i <= Array.GetUpperBound(0); i++)
                {
                    if (Array.Rank == 2) // array is 2D
                    {
                        for (int j = Array.GetLowerBound(1); j <= Array.GetUpperBound(1); j++)
                        {
                            array[y, x] == Array[i,j];
                            x++; y++;
                        }
                    }
                    //TODO: remove if Array is forced to be 2D
                    else // array is 1D
                    {
                        array[y, 1] == ArrayIdea;
                        y++;
                    }
                }
                #endregion
            }

            /// <summary>
            /// Creates a new Matrix of a given size.
            /// After the size got set it can't get changed.
            /// </summary>
            /// <param name="m">Number of Columns.</param>
            /// <param name="n">Number of Rows.</param>
            public NMatrix(uint m, uint n)
            {
                // check size of the matrix
                if(m=0 || n=0)
                {
                    string message = String.Format(
                        "Size of a Matrix can't be zero.{0}"+
                        "Number of Columns was {1} and the Number of Rows was {2}.",
                        Environment.NewLine, m, n);
                    throw new ArgumentNullException(message);
                }

                // create a new array of the given Type
                uint[] lowerBonds = { 1, 1 }; // index starting at 1,1 instead of 0,0
                uint[] lengths = { n, m }; // size of the array
                array = Array.CreateInstance(typeof(T), lengths, lowerBonds);

                // initialize the array
                for(uint i = 1; i <= m; i++)
                {
                    for (uint j = 1; j <= n; j++)
                    {
                        array[i, j] = (T)0;
                    }
                }

            }

            /// <summary>Creates a new Identity-Matrix.</summary>
            /// <param name="Dimension">The Dimension N the matrix shall have.</param>
            /// <returns>The NxN Identity-Matrix.</returns>
            public static NMatrix<T> IdentityMatrix(uint Dimension)
            {
                // check arguments
                if(Dimension == 0)
                {
                    throw new ArgumentNullException("Dimension of a Identity-Matrix can't be zero.");
                }

                // generate a new empty matrix
                NMatrix matrix = new NMatrix<T>(Dimension, Dimension);

                // fill diagonal with ones
                for (uint i = 1; i <= Dimension; i++)
                {
                    matrix[i, i] = (T)1;
                }

                // return the identity-matrix
                return matrix;
            }

            #endregion
            #region Properties
            
            /// <summary>
            /// Gives access to a specific element in the matrix.
            /// </summary>
            /// <param name="i">The Number of the column. First column is 1.</param>
            /// <param name="j">The Number of the row. First row is 1. </param>
            /// <returns>The Element on Position [i,j].</returns>
            public T this[uint i, uint j]
            {
                get
                {
                    CheckPosition(i, j);
                    array[i, j] = value;
                }
                set
                {
                    CheckPosition(i, j);
                    array[i, j] = value;
                }
            }

            /// <summary>Helper function to determine if a given position is in the range of the Matrix.</summary>
            /// <param name="i">The Row of the matrix to access.</param>
            /// <param name="j">The Column of the matrix to access.</param>
            private void CheckPosition(uint i, uint j)
            {
                if (i > m || j > n)
                {
                    string message = String.Format(
                        "The given Position is out of the size of the matrix.{0}" +
                        "The matrix has size {1},{2} and the requested position was {3},{4}.",
                        Environment.NewLine, this.m, this.n, i, j);
                    throw new ArgumentOutOfRangeException(message);
                }
                if (i == 0 || j == 0)
                {
                    string message = String.Format(
                        "Index of a Matrix starts by definition at 1,1.{0}" +
                        "The requested position was {1},{2}.",
                        Environment.NewLine, i, j);
                    throw new ArgumentNullException(message);  
                }
            }

            /// <summary>Returns the dimension (size) of the matrix. </summary>
            /// <param name="m">The number of Columns the matrix has. </param>
            /// <param name="n">The number of Rows the matrix has. </param>
            public void Dimension(out uint m, out uint n)
            {
                m = this.m;
                n = this.n;
            }

            public uint Dimension
            {
                get{
                    if(!IsSquare())
                    {
                        throw new ArithmeticException("Matrix is not square." +
                                  "Please use \"Dimension(out i, out j)\" instead.");
                    }

                    return this.m;
                }
            }

            /// <summary>Determines the Rank of the Matrix.</summary>
            public uint Rank
            {
                get
                {
                    throw new NotImplementedException();
                }
            }
            #endregion
            #region Operators

            /// <summary>Adds two Matrices of same Dimensions. </summary>
            /// <param name="A">The one matrix to add.</param>
            /// <param name="B">The other matrix to add.</param>
            /// <returns>A new matrix that is the sum of the other two.</returns>
            /// <remarks>
            /// This Function is inefficient because it creates a new Matrix.
            /// Use <code>A.Add(B)</code> or <code>A + B</code> where possible.
            /// </remarks>
            public static NMatrix<T> Add(NMatrix<T> A, NMatrix<T> B)
            {
                throw NotImplementedException();
                NMatrix<T> C = (NMatrix<T>) A.Clone();
                return C.Add(B);
            }

            /// <summary>Calls the Add-Function to add another Matrix of same Dimensions to this one. </summary>
            /// <param name="A">The matrix containing the values to add.</param>
            public void operator +(NMatrix<T> A)
            {
                this.Add(A);
            }

            /// <summary>Adds another Matrix of same Dimensions to this one. </summary>
            /// <param name="A">The matrix containing the values to add.</param>
            public void Add(NMatrix<T> A)
            {
                throw new NotImplementedException();

                uint tm, tn, am, an;

                // check if the matrices have the same dimensions
                this.Dimension(out tm, out tn);
                A.Dimension(out am, an);

                if (tm != am || tn != am)
                {
                    string message = String.Format(
                        "The Matrices can't get added, because they don't have equal dimensions.{0}" +
                        "This Matrix has Dimension {1},{2}, while the other Matrix has Dimension {3},{4}.",
                        Environment.NewLine, tm, tn, am, an);
                    throw new ArithmeticException(message);
                }

                // sum up the elements.
                for (uint i = 1; i <= tm; i++)
                {
                    for (uint j = 1; j <= tn; j++)
                    {
                        array[i, j] += A[i, j];
                    }
                }
            }

            #endregion
        }
    }
  • 05-04-2007 13:25 In reply to

    • MovGP0
    • Top 10 Contributor
      Male
    • Joined on 11-21-2005
    • Austria
    • Posts 49

    Re: Matrix Class

    /*
     * Generic Matrix Class for Math.NET
     * written by MovGP0 (mailto:movgp0@gmail.com)
     * licensed under the terms of the CC GNU LGPL 2.1
     * (http://creativecommons.org/licenses/LGPL/2.1/)
     */
    using System;
    using System.Collections.Generic;
    using System.Linq;
    using System.Text;

    namespace MathNet.Numerics.LinearAlgebra
    {
        // implementation of algorithms for calculating Determinants
        partial class NMatrix<T>
        {
            /// <summary>
            /// Returns the Determinant of the Matrix.
            /// The Matrix needs to be square.
            /// </summary>
            public T Determinant
            {
                get
                {
                    // Check if the matrix is square
                    if (!IsSquare())
                    {
                        string message = String.Format(
                            "This is not a square matrix.{0}" +
                            "Determinants are only defined for square matrices.",
                            Environment.NewLine);
                        throw new ArithmeticException(message);
                    }

                    // Matrix is a scalar
                    if (m == 1)
                    {
                        return array[1, 1];
                    }

                    // Rule of Sarrus; used for the determinant of 2x2 or 3x3 matrices
                    if (m == 2)
                    {
                        return RuleOfSarrus2();
                    }
                    if (m == 3)
                    {
                        return RuleOfSarrus3();
                    }

                    // TODO: better Possibilities to calculate Determinants than the Laplace Expansion
                    // * QR decomposition (Requirement: Elements are real; needs gram-schmidt-orthogonalisation)
                    // * Cholesky decomposition (Requirement: Hermition and positive-definite)
                    // * Leibniz formula (Requirement: needs permutations of (1,...,n))
                    // * LU decompositon (Runtime: n³)

                    // if nothing else is working, than this one will solve it
                    return LaplaceExpansion();
                }
            }

            /// <summary>
            /// Check if the matrix is a square Matrix.
            /// This is true if the number of columns is equal to the number of rows.
            /// </summary>
            /// <returns>
            /// Returns true if the Matrix is square.
            /// Otherwise it returns false.
            /// </returns>
            private bool IsSquare()
            {
                return this.m == this.n;
            }

            /// <summary>Calculates Rule of Sarrus for a 2x2 Matrix.</summary>
            private T RuleOfSarrus2()
            {
                return array[1, 1] * array[2, 2] - array[2, 1] * array[1, 2];
            }

            /// <summary>Calculates Rule of Sarrus for a 2x2 Matrix.</summary>
            private T RuleOfSarrus3()
            {
                return array[1, 1] * array[2, 2] * array[3, 3]
                     + array[1, 2] * array[2, 3] * array[3, 1]
                     + array[1, 3] * array[2, 1] * array[3, 2]
                     - array[1, 3] * array[2, 2] * array[3, 1]
                     - array[1, 1] * array[2, 3] * array[3, 2]
                     - array[1, 2] * array[2, 1] * array[3, 3];
            }

            /// <summary>
            /// Calculates the Determinant using the Laplace Expansion Algorithm.
            /// Notice, that this algorithm has a Runtime of n! and therefore
            /// should only get used for small matrices.
            /// </summary>
            private T LaplaceExpansion()
            {
                throw new NotImplementedException();

                T sum = 0; // initialzize sum
                uint j = 1; // the first Vector is given everytime; others may not

                // evolve by j-th Column (Vector)
                for (i = 1; i <= n; i++) // the loop of the Sum
                {
                    // check if value is zero
                    // if so there is no need to calculate the cofacotor,
                    // because the product will always be zero
                    T value = array[i, j];
                    if (value == (T)0) continue;

                    // calculate the Laplace
                    sum += value * Cofactor(i, j);
                }
            }

            /// <summary>Calculates the cofactor (minor) of the matrix.</summary>
            /// <param name="i">The Column that is not needed.</param>
            /// <param name="j">The Row that is not needed.</param>
            /// <returns>The Cofactor a<sub>i,j</sub>.</returns>
            public T Cofactor(uint i, uint j)
            {
                return ((-1) ^ (i + j)) // the sign
                    * CramersMatrix(i, j);
            }

            /// <summary>Alias for Cofactor.</summary>
            /// <param name="i">The Column that is not needed.</param>
            /// <param name="j">The Row that is not needed.</param>
            /// <returns>The Minor.</returns>
            public T Minor(uint i, uint j)
            {
                return Cofactor(i, j);
            }

            /// <summary>This returns the Matrix without the i-th Column and without the j-th Row.</summary>
            /// <param name="i">The Column which gets removed.</param>
            /// <param name="j">The Row which gets removed.</param>
            /// <returns>Cramer's Matrix (C<sub>i,j</sub>)</returns>
            public NMatrix<T> CramersMatrix(uint i, uint j)
            {
                CheckPosition(i, j);

                if (this.m == 1 || this.n == 1)
                {
                    throw new ArithmeticException("Matrix is to small to calculate a Cramers Matrix.");
                }

                NMatrix<T> CramerMatrix = new NMatrix<T>(this.m - 1, this.n - 1);

                uint c, r;

                for (uint col = 0; col < this.m; col++)
                {
                    if (col = i) continue; // skip columns
                    (col > i) ? c = col - 1 : c = col; // c is next column

                    for (uint row = 0; row < this.m; row++)
                    {
                        if (row == j) continue; // skip row
                        (row > j) ? r = row - 1 : r = row; // r is next row

                        CramerMatrix[col, row] = array[c, r]; // copy value into the new array
                    }
                }

                return CramerMatrix;
            }
        }
    }

  • 05-04-2007 13:27 In reply to

    • MovGP0
    • Top 10 Contributor
      Male
    • Joined on 11-21-2005
    • Austria
    • Posts 49

    Matrix Class Interface Implementation

    /*
     * Generic Matrix Class for Math.NET
     * written by MovGP0 (mailto:movgp0@gmail.com)
     * licensed under the terms of the CC GNU LGPL 2.1
     * (http://creativecommons.org/licenses/LGPL/2.1/)
     */
    using System;
    using System.Collections.Generic;
    using System.Linq;
    using System.Text;

    namespace MathNet.Numerics.LinearAlgebra
    {
        // implementation of interfaces
        partial class NMatrix<T> : ICloneable //, IEnumerable<T>, ICollection<T>, IList<T>
        {
            #region ICloneable Members

            //object ICloneable.Clone()
            //{
            //    uint m, n;
            //    this.Dimension(out m, out n);
            //    NMatrix<T> A = new NMatrix<T>(m, n);

            //    for (i = 1; i <= m; i++)
            //    {
            //        for (j = 1; j <= n; j++)
            //        {
            //            A[i, j] = this[i, j];
            //        }
            //    }
            //    return (object)A;
            //}

            /// <summary>Creates a memberwise Deep-Copy. The returned Object is of type NMatrix<T>.</summary>
            /// <remarks>Because of extensive value checking this method is not very efficient.</remarks>
            /// <returns>A copy of the Matrix.</returns>
            object ICloneable.Clone()
            {
                return new NMatrix<T>(array);
            }

            #endregion
        }
    }
  • 05-05-2007 12:51 In reply to

    Re: Rewrite of Matrix Class

    Interesting, does Linq* somehow solve the problem with Generics vs. ValueTypes vs. Calculations? How does it do that?

    *or rather: the new v3.5 CLR (and C# compiler) versions

  • 05-05-2007 21:02 In reply to

    • Andrii
    • Top 10 Contributor
      Male
    • Joined on 05-05-2007
    • Ukraine, Odesa
    • Posts 8

    Re: Rewrite of Matrix Class

    I have some idea.

    Could you write Complex Matrix Class for complex number or maybe template class for supporting both number's types.

  • 05-06-2007 20:35 In reply to

    Re: Rewrite of Matrix Class

    Christoph Rüegg:

    Interesting, does Linq* somehow solve the problem with Generics vs. ValueTypes vs. Calculations? How does it do that?

    *or rather: the new v3.5 CLR (and C# compiler) versions

     

    Just to clarify what I meant: How come something like "array[i, j] += A[i, j];" actually compiles when array is of type T[,] of a T with no restrictions?

Page 1 of 1 (7 items)
Powered by Community Server (Non-Commercial Edition), by Telligent Systems