Core.inl

Go to the documentation of this file.

//Copyright (c) 2004-2005, Baris Sumengen
//All rights reserved.
//
// CIMPL Matrix Performance Library
//
//Redistribution and use in source and binary
//forms, with or without modification, are
//permitted provided that the following
//conditions are met:
//
//    * No commercial use is allowed. 
//    This software can only be used
//    for non-commercial purposes. This 
//    distribution is mainly intended for
//    academic research and teaching.
//    * Redistributions of source code must
//    retain the above copyright notice, this
//    list of conditions and the following
//    disclaimer.
//    * Redistributions of binary form must
//    mention the above copyright notice, this
//    list of conditions and the following
//    disclaimer in a clearly visible part 
//    in associated product manual, 
//    readme, and web site of the redistributed 
//    software.
//    * Redistributions in binary form must
//    reproduce the above copyright notice,
//    this list of conditions and the
//    following disclaimer in the
//    documentation and/or other materials
//    provided with the distribution.
//    * The name of Baris Sumengen may not be
//    used to endorse or promote products
//    derived from this software without
//    specific prior written permission.
//
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
//HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
//EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
//NOT LIMITED TO, THE IMPLIED WARRANTIES OF
//MERCHANTABILITY AND FITNESS FOR A PARTICULAR
//PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
//CONTRIBUTORS BE LIABLE FOR ANY
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
//EXEMPLARY, OR CONSEQUENTIAL DAMAGES
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
//OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
//DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
//HOWEVER CAUSED AND ON ANY THEORY OF
//LIABILITY, WHETHER IN CONTRACT, STRICT
//LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
//OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
//OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//POSSIBILITY OF SUCH DAMAGE.


namespace MathCore
{



        template <class T>
        Vector<int> Find(Vector<T>& m)
        {
                std::vector<int> ind;
                for(int i=0; i<m.Length(); i++)
                {
                        if(m.ElemNC(i) != 0)
                        {
                                ind.push_back(i);
                        }
                }
                Vector<int> temp((int)ind.size());
                memcpy(temp.Data(), &ind[0], sizeof(int)*ind.size() );
                //std::vector<int>::const_iterator constIterator;
                //int j = 0;
                //for(constIterator = ind.begin(); constIterator != ind.end(); constIterator++)
                //{
                //      temp[j] = *constIterator;
                //      j++;
                //}
                return temp;
        }
        
        template <class T>
        Vector<int> Find(Matrix<T>& m)
        {
                std::vector<int> ind;
                for(int i=0; i<m.Length(); i++)
                {
                        if(m.ElemNC(i) != 0)
                        {
                                ind.push_back(i);
                        }
                }
                Vector<int> temp((int)ind.size());
                memcpy(temp.Data(), &ind[0], sizeof(int)*ind.size() );
                return temp;
        }
        
        template <class T>
        void Find(Matrix<T>& m, Vector<int>& I, Vector<int>& J)
        {
                std::vector<int> indI;
                std::vector<int> indJ;
                for(int j=0; j<m.Columns(); j++)
                {
                        for(int i=0; i<m.Rows(); i++)
                        {
                                if(m.ElemNC(i,j) != 0)
                                {
                                        indI.push_back(i);
                                        indJ.push_back(j);
                                }
                        }
                }
                I = Vector<int>((int)indI.size());
                memcpy(I.Data(), &indI[0], sizeof(int)*indI.size() );
                J = Vector<int>((int)indJ.size());
                memcpy(J.Data(), &indJ[0], sizeof(int)*indJ.size() );
        }


        
        // Vectors need to be of length 3
        template <class T>
        Vector<T> Cross(Vector<T>& x, Vector<T>& y)
        {
                Vector<T> temp(3);
                temp[0] = x[1]*y[2]-x[2]*y[1];
                temp[1] = x[2]*y[0]-x[0]*y[2];
                temp[2] = x[0]*y[1]-x[1]*y[0];
                
                return temp;
        }

        
        template <class T>
        Matrix<T> Cross(Matrix<T>& x, Matrix<T>& y)
        {
                if(x.Rows() == 3)
                {
                        return Cross(x,y,1);
                }
                else if(x.Columns() == 3)
                {
                        return Cross(x,y,2);
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("One of the dimension of the matrices x and y should be of length 3!");
                }
        }

        
        template <class T>
        Matrix<T> Cross(Matrix<T>& x, Matrix<T>& y, int dimension)
        {
                if(x.Rows() != y.Rows() || x.Columns() != y.Columns())
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Matrices are not of the same size!");
                }
                
                Matrix<T> temp(x.Rows(), x.Columns());
                
                if(dimension == 1)
                {
                        if(temp.Rows() != 3)
                        {
                                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                                Utility::RunTimeError("Both matrices x and y should be of length 3 along the dimension on which the cross product is taken!");
                        }
                        for(int z=0; z<temp.Columns(); z++)
                        {
                                temp.ElemNC(0,z) = x.ElemNC(1,z)*y.ElemNC(2,z) - x.ElemNC(2,z)*y.ElemNC(1,z);
                                temp.ElemNC(1,z) = x.ElemNC(2,z)*y.ElemNC(0,z) - x.ElemNC(0,z)*y.ElemNC(2,z);
                                temp.ElemNC(2,z) = x.ElemNC(0,z)*y.ElemNC(1,z) - x.ElemNC(1,z)*y.ElemNC(0,z);
                        }
                }
                else if(dimension == 2)
                {
                        if(temp.Columns() != 3)
                        {
                                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                                Utility::RunTimeError("Both matrices x and y should be of length 3 along the dimension on which the cross product is taken!");
                        }
                        for(int z=0; z<temp.Rows(); z++)
                        {
                                temp.ElemNC(z,0) = x.ElemNC(z,1)*y.ElemNC(z,2) - x.ElemNC(z,2)*y.ElemNC(z,1);
                                temp.ElemNC(z,1) = x.ElemNC(z,2)*y.ElemNC(z,0) - x.ElemNC(z,0)*y.ElemNC(z,2);
                                temp.ElemNC(z,2) = x.ElemNC(z,0)*y.ElemNC(z,1) - x.ElemNC(z,1)*y.ElemNC(z,0);
                        }
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
                }
                return temp;
        }


        
        template <class T>
        Vector<T> CumSum(Vector<T>& m)
        {
                Vector<T> temp(m.Length());
                temp.ElemNC(0) = m.ElemNC(0);
                for(int i=1; i<m.Length(); i++)
                {
                        temp.ElemNC(i) = temp.ElemNC(i-1) + m.ElemNC(i);
                }
                return temp;
        }
        
        template <class T>
        Matrix<T> CumSum(Matrix<T>& m)
        {
                return CumSum(m,1);
        }
        
        template <class T>
        Matrix<T> CumSum(Matrix<T>& m, int dimension)
        {
                Matrix<T> temp(m.Rows(), m.Columns());
                if(dimension == 1)
                {
                        for(int z=0; z<m.Columns(); z++)
                        {
                                temp.ElemNC(0,z) = m.ElemNC(0,z);
                                for(int i=1; i<m.Rows(); i++)
                                {
                                        temp.ElemNC(i,z) = temp.ElemNC(i-1,z) + m.ElemNC(i,z);
                                }
                        }
                }
                else if(dimension == 2)
                {
                        for(int z=0; z<m.Rows(); z++)
                        {
                                temp.ElemNC(z,0) = m.ElemNC(z,0);
                                for(int i=1; i<m.Columns(); i++)
                                {
                                        temp.ElemNC(z,i) = temp.ElemNC(z,i-1) + m.ElemNC(z,i);
                                }
                        }
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
                }       
                return temp;
        }


        
        
        template <class T>
        Vector<T>& CumSumI(Vector<T>& m)
        {
                for(int i=1; i<m.Length(); i++)
                {
                        m.ElemNC(i) = m.ElemNC(i-1) + m.ElemNC(i);
                }
                return m;
        }
        
        template <class T>
        Matrix<T>& CumSumI(Matrix<T>& m)
        {
                return CumSum(m,1);
        }
        
        template <class T>
        Matrix<T>& CumSumI(Matrix<T>& m, int dimension)
        {
                if(dimension == 1)
                {
                        for(int z=0; z<m.Columns(); z++)
                        {
                                for(int i=1; i<m.Rows(); i++)
                                {
                                        m.ElemNC(i,z) = m.ElemNC(i-1,z) + m.ElemNC(i,z);
                                }
                        }
                }
                else if(dimension == 2)
                {
                        for(int z=0; z<m.Rows(); z++)
                        {
                                for(int i=1; i<m.Columns(); i++)
                                {
                                        m.ElemNC(z,i) = m.ElemNC(z,i-1) + m.ElemNC(z,i);
                                }
                        }
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
                }       
                return m;
        }


        
        template <class T>
        Vector<T> CumProd(Vector<T>& m)
        {
                Vector<T> temp(m.Length());
                temp.ElemNC(0) = m.ElemNC(0);
                for(int i=1; i<m.Length(); i++)
                {
                        temp.ElemNC(i) = temp.ElemNC(i-1) * m.ElemNC(i);
                }
                return temp;
        }
        
        template <class T>
        Matrix<T> CumProd(Matrix<T>& m)
        {
                return CumProd(m,1);
        }
        
        template <class T>
        Matrix<T> CumProd(Matrix<T>& m, int dimension)
        {
                Matrix<T> temp(m.Rows(), m.Columns());
                if(dimension == 1)
                {
                        for(int z=0; z<m.Columns(); z++)
                        {
                                temp.ElemNC(0,z) = m.ElemNC(0,z);
                                for(int i=1; i<m.Rows(); i++)
                                {
                                        temp.ElemNC(i,z) = temp.ElemNC(i-1,z) * m.ElemNC(i,z);
                                }
                        }
                }
                else if(dimension == 2)
                {
                        for(int z=0; z<m.Rows(); z++)
                        {
                                temp.ElemNC(z,0) = m.ElemNC(z,0);
                                for(int i=1; i<m.Columns(); i++)
                                {
                                        temp.ElemNC(z,i) = temp.ElemNC(z,i-1) * m.ElemNC(z,i);
                                }
                        }
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
                }       
                return temp;
        }

        
        template <class T>
        Vector<T>& CumProdI(Vector<T>& m)
        {
                for(int i=1; i<m.Length(); i++)
                {
                        m.ElemNC(i) = m.ElemNC(i-1) * m.ElemNC(i);
                }
                return m;
        }
        
        template <class T>
        Matrix<T>& CumProdI(Matrix<T>& m)
        {
                return CumProd(m,1);
        }
        
        template <class T>
        Matrix<T>& CumProdI(Matrix<T>& m, int dimension)
        {
                if(dimension == 1)
                {
                        for(int z=0; z<m.Columns(); z++)
                        {
                                for(int i=1; i<m.Rows(); i++)
                                {
                                        m.ElemNC(i,z) = m.ElemNC(i-1,z) * m.ElemNC(i,z);
                                }
                        }
                }
                else if(dimension == 2)
                {
                        for(int z=0; z<m.Rows(); z++)
                        {
                                for(int i=1; i<m.Columns(); i++)
                                {
                                        m.ElemNC(z,i) = m.ElemNC(z,i-1) * m.ElemNC(z,i);
                                }
                        }
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
                }       
                return m;
        }

        
        template <class T>
        Matrix<T> Diag(Vector<T>& m)
        {
                return Diag(m,0);
        }
        
        template <class T>
        Matrix<T> Diag(Vector<T>& m, int k)
        {
                int absk = abs(k);
                int dim = m.Length() + absk;
                Matrix<T> temp(dim,dim,0);
                if(k>=0)
                {
                        for(int i=0; i<m.Length(); i++)
                        {
                                temp.ElemNC(i,i+absk) = m[i];
                        }
                }
                else
                {
                        for(int i=0; i<m.Length(); i++)
                        {
                                temp.ElemNC(i+absk,i) = m[i];
                        }
                }
                return temp;
        }

        
        template <class T>
        Vector<T> Diag(Matrix<T>& m)
        {
                return Diag(m,0);
        }

        
        template <class T>
        Vector<T> Diag(Matrix<T>& m, int k)
        {
                Vector<T> temp;
                int absk = abs(k);
                int dim;
                if(k>=0)
                {
                        if(absk>=m.Columns())
                        {
                                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                                Utility::RunTimeError("diagonal index k is outside the matrix dimensions!");
                        }
                        dim = m.Columns()-absk > m.Rows() ? m.Rows() : m.Columns()-absk;
                        temp = Vector<T>(dim);
                        for(int i=0; i<temp.Length(); i++)
                        {
                                temp[i] = m.ElemNC(i,i+absk);
                        }
                }
                else
                {
                        if(absk>=m.Rows())
                        {
                                cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                                Utility::RunTimeError("diagonal index k is outside the matrix dimensions!");
                        }
                        dim = m.Rows()-absk > m.Columns() ? m.Columns() : m.Rows()-absk;
                        temp = Vector<T>(dim);
                        for(int i=0; i<temp.Length(); i++)
                        {
                                temp[i] = m.ElemNC(i+absk,i);
                        }
                }
                
                return temp;
        }


        
        template <class T>
        Matrix<T> FlipLR(Matrix<T>& m)
        {
                Matrix<T> temp(m.Rows(), m.Columns());
                for(int j=0; j<m.Columns(); j++)
                {
                        for(int i=0; i<m.Rows(); i++)
                        {
                                temp.ElemNC(i,m.Columns()-j-1) = m.ElemNC(i,j);
                        }
                }
                return temp;
        }

        
        template <class T>
        Matrix<T>& FlipLRI(Matrix<T>& m)
        {
                for(int j=0; j<m.Columns()/2; j++)
                {
                        for(int i=0; i<m.Rows(); i++)
                        {
                                 T dummy = m.ElemNC(i,j);
                                 m.ElemNC(i,j) = m.ElemNC(i,m.Columns()-j-1);
                                 m.ElemNC(i,m.Columns()-j-1) = dummy;
                        }
                }
                return m;
        }


        
        template <class T>
        Matrix<T> FlipUD(Matrix<T>& m)
        {
                Matrix<T> temp(m.Rows(), m.Columns());
                for(int j=0; j<m.Columns(); j++)
                {
                        for(int i=0; i<m.Rows(); i++)
                        {
                                temp.ElemNC(m.Rows()-i-1,j) = m.ElemNC(i,j);
                        }
                }
                return temp;
        }

        
        template <class T>
        Matrix<T>& FlipUDI(Matrix<T>& m)
        {
                for(int j=0; j<m.Columns(); j++)
                {
                        for(int i=0; i<m.Rows()/2; i++)
                        {
                                T dummy = m.ElemNC(i,j);
                                m.ElemNC(i,j) = m.ElemNC(m.Rows()-i-1,j);
                                m.ElemNC(m.Rows()-i-1,j) = dummy;
                        }
                }
                return m;
        }


        
        template <class T>
        Vector<T> Reverse(Vector<T>& m)
        {
                Vector<T> temp(m.Length());
                for(int i=0; i<m.Length(); i++)
                {
                        temp.ElemNC(m.Length()-i-1) = m.ElemNC(i);
                }
                return temp;
        }

        
        template <class T>
        Vector<T>& ReverseI(Vector<T>& m)
        {
                for(int i=0; i<(m.Length()+1)/2; i++)
                {
                        T dummy = m.ElemNC(i);
                        m.ElemNC(i) = m.ElemNC(m.Length()-i-1);
                        m.ElemNC(m.Length()-i-1) = dummy;
                }
                return m;
        }


        template <class T>
        Matrix<T> Rot90(Matrix<T>& m)
        {
                return Rot90(m,1);
        }

        
        template <class T>
        Matrix<T> Rot90(Matrix<T>& m, int k)
        {
                int rot = k % 4;
                if(rot<0)
                {
                        rot += 4;
                }
                
                if(rot==0)
                {
                        return m.Clone();
                }
                else if(rot==1)
                {
                        return FlipUD(~m);
                }
                else if(rot==2)
                {
                        return FlipLR(FlipUD(m));
                }
                else if(rot==3)
                {
                        return FlipLR(~m);
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Something went wrong when trying to rotate the matrix!");
                }
        }
        
        template <class T>
        Vector<T> Maximum(Matrix<T>& m)
        {
                Vector<T> maxes(m.Columns());
                for(int j=0; j<m.Columns(); j++)
                {
                        maxes[j] = Maximum(m[j]);
                }
                return maxes;
        }
        
        template <class T>
        T Maximum(Vector<T>& m)
        {
                T max = m.ElemNC(0);
                for(int i=1; i<m.Length(); i++)
                {
                        max = max >= m.ElemNC(i) ? max : m.ElemNC(i);
                }
                return max;
        }
        
        template <class T>
        Vector<T> Minimum(Matrix<T>& m)
        {
                Vector<T> mins(m.Columns());
                for(int j=0; j<m.Columns(); j++)
                {
                        mins[j] = Minimum(m[j]);
                }
                return mins;
        }
        
        template <class T>
        T Minimum(Vector<T>& m)
        {
                T min = m.ElemNC(0);
                for(int i=1; i<m.Length(); i++)
                {
                        min = min <= m.ElemNC(i) ? min : m.ElemNC(i);
                }
                return min;
        }

        
        template <class T>
        Vector<double> Mean(Matrix<T>& m)
        {
                Vector<double> means(m.Columns());
                for(int j=0; j<m.Columns(); j++)
                {
                        means[j] = Mean(m[j]);
                }
                return means;
        }
        
        template <class T>
        double Mean(Vector<T>& m)
        {
                double mean = 0;
                for(int i=0; i<m.Length(); i++)
                {
                        mean += m.ElemNC(i);
                }
                mean /= m.Length();
                return mean;
        }
        
        template <class T>
        Vector<T> Median(Matrix<T>& m)
        {
                Vector<T> medians(m.Columns());
                for(int j=0; j<m.Columns(); j++)
                {
                        medians[j] = Median(m[j]);
                }
                return medians;
        }
        
        template <class T>
        T Median(Vector<T>& m)
        {
                Vector<T> temp = m.Clone();
                T *arr = temp.Data();
                int n = temp.Length();
                int low, high;
                int median;
                int middle, ll, hh;

                low = 0;
                high = m.Length()-1;
                median = (low + high)/2;
                for(;;)
                {
                        if(high <= low)
                        {
                                return arr[median];
                        }
                        if(high == low + 1)
                        {
                                if(arr[low] > arr[high])
                                {
                                        Swap(arr[low], arr[high]);
                                }
                                return arr[median];
                        }

                        middle = (low + high)/2;
                        if(arr[middle] > arr[high])
                        {
                                Swap(arr[middle], arr[high]);
                        }
                        if(arr[low] > arr[high])
                        {
                                Swap(arr[low], arr[high]);
                        }
                        if(arr[middle] > arr[low])
                        {
                                Swap(arr[middle], arr[low]);
                        }

                        Swap(arr[middle], arr[low+1]) ;

                        ll = low + 1;
                        hh = high;
                        for(;;)
                        {
                                do ll++; while(arr[low] > arr[ll]);
                                do hh--; while(arr[hh]  > arr[low]);

                                if (hh < ll)
                                {
                                        break;
                                }

                                Swap(arr[ll], arr[hh]);
                        }

                        Swap(arr[low], arr[hh]);

                        if (hh <= median)
                        {
                                low = ll;
                        }
                        if (hh >= median)
                        {
                                high = hh - 1;
                        }
                }
        }

        
        template <class T>
        Vector<double> Var(Matrix<T>& m)
        {
                return Var(m,true);
        }
        
        template <class T>
        double Var(Vector<T>& m)
        {
                return Var(m,true);
        }

        
        template <class T>
        Vector<double> Var(Matrix<T>& m, bool unbiased)
        {
                Vector<double> vars(m.Columns());
                for(int j=0; j<m.Columns(); j++)
                {
                        vars[j] = Var(m[j], unbiased);
                }
                return vars;
        }
        
        template <class T>
        double Var(Vector<T>& m, bool unbiased)
        {
                double mean = Mean(m);
                double var = 0;
                for(int i=0; i<m.Length(); i++)
                {
                        var += (m.ElemNC(i)-mean)*(m.ElemNC(i)-mean);
                }
                
                if(unbiased)
                {
                        return var/(m.Length()-1);
                }
                else
                {
                        return var/m.Length();
                }
        }


        
        template <class T>
        Vector<double> Std(Matrix<T>& m)
        {
                return Std(m,true);
        }
        
        template <class T>
        double Std(Vector<T>& m)
        {
                return Std(m,true);
        }

        
        template <class T>
        Vector<double> Std(Matrix<T>& m, bool unbiased)
        {
                Vector<double> stds(m.Columns());
                for(int j=0; j<m.Columns(); j++)
                {
                        stds[j] = Std(m[j], unbiased);
                }
                return stds;
        }
        
        template <class T>
        double Std(Vector<T>& m, bool unbiased)
        {
                return sqrt(Var(m,unbiased));
        }



        
        template <class T>
        Vector<T> Sum(Matrix<T>& m)
        {
                return Sum(m,1);
        }

        
        template <class T>
        Vector<T> Sum(Matrix<T>& m, int dimension)
        {
                Vector<T> sums;
                if(dimension == 1)
                {
                        sums = Vector<T>(m.Columns());
                        for(int j=0; j<m.Columns(); j++)
                        {
                                sums[j] = Sum(m[j]);
                        }
                }
                else if(dimension == 2)
                {
                        sums = Vector<T>(m.Rows());
                        for(int j=0; j<m.Rows(); j++)
                        {
                                T sum = 0;
                                for(int i=0; i<m.Columns(); i++)
                                {
                                        sum += m.ElemNC(j,i);
                                }
                                sums[j] = sum;
                        }
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
                }
                return sums;
        }
        
        template <class T>
        T Sum(Vector<T>& m)
        {
                T sum = 0;
                for(int i=0; i<m.Length(); i++)
                {
                        sum += m.ElemNC(i);
                }
                return sum;
        }



        
        template <class T>
        Vector<T> Prod(Matrix<T>& m)
        {
                return Prod(m,1);
        }

        
        template <class T>
        Vector<T> Prod(Matrix<T>& m, int dimension)
        {
                Vector<T> prods;
                if(dimension == 1)
                {
                        prods = Vector<T>(m.Columns());
                        for(int j=0; j<m.Columns(); j++)
                        {
                                prods[j] = Prod(m[j]);
                        }
                }
                else if(dimension == 2)
                {
                        prods = Vector<T>(m.Rows());
                        for(int j=0; j<m.Rows(); j++)
                        {
                                T prod = 1;
                                for(int i=0; i<m.Columns(); i++)
                                {
                                        prod *= m.ElemNC(j,i);
                                }
                                prods[j] = prod;
                        }
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
                }
                return prods;
        }
        
        template <class T>
        T Prod(Vector<T>& m)
        {
                T prod = 1;
                for(int i=0; i<m.Length(); i++)
                {
                        prod *= m.ElemNC(i);
                }
                return prod;
        }

        
        template <class T>
        Vector<T> Diff(Vector<T>& m)
        {
                if(m.Length() == 1)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Diff can only be called on vectors of length 2 or more!");
                }
                Vector<T> temp(m.Length()-1);
                for(int i=0;i<temp.Length(); i++)
                {
                        temp.ElemNC(i) = m.ElemNC(i+1) - m.ElemNC(i);
                }
                return temp;
        }
        
        template <class T>
        Matrix<T> Diff(Matrix<T>& m)
        {
                if(m.Rows() == 1)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Diff can only be called on matrices with 2 or more rows!");
                }
                Matrix<T> temp(m.Rows()-1, m.Columns());
                for(int j=0; j<temp.Columns(); j++)
                {
                        for(int i=0; i<temp.Rows(); i++)
                        {
                                
                                temp.ElemNC(i,j) = m.ElemNC(i+1,j) - m.ElemNC(i,j);
                        }
                }
                return temp;
        }



        
        template <class T>
        Matrix<T> RepMat(Matrix<T>& m, int M, int N)
        {
                Matrix<T> temp(M*m.Rows(),N*m.Columns());
                for(int k=0; k<M; k++)
                {
                        for(int l=0; l<N; l++)
                        {
                                int s = k*m.Rows();
                                int t = l*m.Columns();
                                for(int i=0; i<m.Rows(); i++)
                                {
                                        for(int j=0; j<m.Columns(); j++)
                                        {
                                                temp.ElemNC(s+i,t+j) = m.ElemNC(i,j);
                                        }
                                }               
                        }
                }
                return temp;
        }


        
        template <class T>
        Matrix<T> Reshape(Matrix<T>& m, int M, int N)
        {
                if(m.Length() != M*N)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Number of elements should be the same when using Reshape()!");
                }
                Matrix<T> temp(M,N);
                memcpy(temp.Data(), m.Data(), sizeof(T)*m.Length());
                return temp;
        }

        template <class T>
        Matrix<T> Reshape(Vector<T>& m, int M, int N)
        {
                if(m.Length() != M*N)
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("Number of elements should be the same when using Reshape()!");
                }
                Matrix<T> temp(M,N);
                memcpy(temp.Data(), m.Data(), sizeof(T)*m.Length());
                return temp;
        }
        
        template <class T>
        Vector<T>& QuickSort(Vector<T>& m)
        {
                int M = 7;
                int NSTACK = 50;
                
                int i, j, k; 
                int l=1;
                int ir = m.Length();
                Vector<int> stack(m.Length());
                int *istack = stack.Data()-1;
                int jstack = 0;
                
                T a;
                T* arr = m.Data()-1;
                
                for (;;) 
                {
                        if (ir-l < M) 
                        {
                                for (j=l+1;j<=ir;j++) 
                                {
                                        a=arr[j];
                                        for (i=j-1;i>=l;i--) 
                                        {
                                                if (arr[i] <= a) 
                                                {
                                                        break;
                                                }
                                                arr[i+1]=arr[i];
                                        }
                                        arr[i+1]=a;
                                }
                                if (jstack == 0) 
                                {
                                        break;
                                }
                                ir=istack[jstack--];
                                l=istack[jstack--];
                        } 
                        else 
                        {
                                k=(l+ir) >> 1;
                                Swap(arr[k],arr[l+1]);
                                if (arr[l] > arr[ir]) 
                                {
                                        Swap(arr[l],arr[ir]);
                                }
                                if (arr[l+1] > arr[ir]) 
                                {
                                        Swap(arr[l+1],arr[ir]);
                                }
                                if (arr[l] > arr[l+1]) 
                                {
                                        Swap(arr[l],arr[l+1]);
                                }
                                i=l+1; 
                                j=ir;
                                a=arr[l+1];
                                for (;;) 
                                {
                                        do i++; while (arr[i] < a); 
                                        do j--; while (arr[j] > a); 
                                        if (j < i) 
                                        {
                                                break; 
                                        }
                                        Swap(arr[i],arr[j]);
                                }
                                arr[l+1]=arr[j];
                                arr[j]=a;
                                jstack += 2;
                                if (jstack > NSTACK) 
                                {
                                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                                        Utility::RunTimeError("Stack size is too small in quick sort!");
                                }
                                if (ir-i+1 >= j-l)
                                {
                                        istack[jstack]=ir;
                                        istack[jstack-1]=i;
                                        ir=j-1;
                                } 
                                else
                                {
                                        istack[jstack]=j-1;
                                        istack[jstack-1]=l;
                                        l=i;
                                }
                        }
                }
                return m;
        }
        
                
        template <class T>
        Vector<T> Sort(Vector<T>& m)
        {
                Vector<T> sorted = m.Clone();
                QuickSort(sorted);
                return sorted;
        }
        
        template <class T>
        Matrix<T> Sort(Matrix<T>& m, int dimension)
        {
                Matrix<T> temp;
                if(dimension == 1)
                {
                        temp = m.Clone();
                        for(int j=0; j<temp.Columns(); j++)
                        {
                                QuickSort(temp[j]);
                        }
                }
                else if(dimension == 2)
                {
                        temp = Matrix<T>(m.Rows(), m.Columns());
                        for(int i=0; i<temp.Rows(); i++)
                        {
                                Vector<T> dummy(temp.Columns());
                                for(int j=0; j<temp.Columns(); j++)
                                {
                                        dummy.ElemNC(j) = m.ElemNC(i,j);
                                }
                                QuickSort(dummy);
                                for(int j=0; j<temp.Columns(); j++)
                                {
                                        temp.ElemNC(i,j) = dummy.ElemNC(j);
                                }
                        }
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
                }
                return temp;
        }
        
        template <class T>
        Matrix<T> Sort(Matrix<T>& m)
        {
                return Sort(m,1);
        }




        
        template <class T>
        Vector<T>& SortI(Vector<T>& m)
        {
                return QuickSort(m);
        }
        
        template <class T>
        Matrix<T>& SortI(Matrix<T>& m, int dimension)
        {
                if(dimension == 1)
                {
                        for(int j=0; j<m.Columns(); j++)
                        {
                                QuickSort(m[j]);
                        }
                }
                else if(dimension == 2)
                {
                        Vector<T> dummy(m.Columns());
                        for(int i=0; i<m.Rows(); i++)
                        {
                                for(int j=0; j<m.Columns(); j++)
                                {
                                        dummy.ElemNC(j) = m.ElemNC(i,j);
                                }
                                QuickSort(dummy);
                                for(int j=0; j<m.Columns(); j++)
                                {
                                        m.ElemNC(i,j) = dummy.ElemNC(j);
                                }
                        }
                }
                else
                {
                        cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
                        Utility::RunTimeError("dimension should be either 1 or 2 for matrices!");
                }
                return m;
        }
        
        template <class T>
        Matrix<T>& SortI(Matrix<T>& m)
        {
                return Sort(m,1);
        }





        
        template <class T>
        void MeshGrid(Vector<T>& x, Vector<T>& y, Matrix<T>& X, Matrix<T>& Y)
        {
                X = Matrix<T>(y.Length(), x.Length());
                Y = Matrix<T>(y.Length(), x.Length());
                for(int i=0; i<y.Length(); i++)
                {
                        for(int j=0; j<x.Length(); j++)
                        {
                                X.ElemNC(i,j) = x(j);
                                Y.ElemNC(i,j) = y(i);
                        }
                }
        }





        
        
        template <class T>
        Vector<int> Sign(Vector<T>& m)
        {
                Vector<int> temp(m.Length());
                for(int i=0; i<m.Length(); i++)
                {
                        if(m.ElemNC(i) > 0)
                        {
                                temp.ElemNC(i) = 1;
                        }
                        else if(m.ElemNC(i) < 0)
                        {
                                temp.ElemNC(i) = -1;
                        }
                        else
                        {
                                temp.ElemNC(i) = 0;
                        }
                }
                return temp;
        }
        
        template <class T>
        Matrix<int> Sign(Matrix<T>& m)
        {
                Matrix<int> temp(m.Rows(), m.Columns());
                for(int i=0; i<m.Length(); i++)
                {
                        if(m.ElemNC(i) > 0)
                        {
                                temp.ElemNC(i) = 1;
                        }
                        else if(m.ElemNC(i) < 0)
                        {
                                temp.ElemNC(i) = -1;
                        }
                        else
                        {
                                temp.ElemNC(i) = 0;
                        }
                }
                return temp;
        }

        
        template <class T>
        Vector<double> Gradient(Vector<T>& m)
        {
                Vector<double> temp(m.Length());
                for(int i=1; i<m.Length()-1; i++)
                {
                        temp.ElemNC(i) = (m.ElemNC(i+1)-m.ElemNC(i-1))/2.0;
                }
                temp.ElemNC(0) = m.ElemNC(1)-m.ElemNC(0);
                temp.ElemNC(m.Length()-1) = m.ElemNC(m.Length()-1)-m.ElemNC(m.Length()-2);
                return temp;
        }

        
        template <class T>
        void Gradient(Matrix<T>& m, Matrix<double>& X, Matrix<double>& Y)
        {
                Matrix<double> X(m.Rows(), m.Columns());
                Matrix<double> Y(m.Rows(), m.Columns());
                
                for(int j=1; j<m.Columns()-1; j++)
                {
                        for(int i=0; i<m.Rows(); i++)
                        {
                                X.ElemNC(i,j) = (m.ElemNC(i,j+1)-m.ElemNC(i,j-1))/2.0;
                        }
                }
                for(int i=0; i<m.Rows(); i++)
                {
                        X.ElemNC(i,0) = m.ElemNC(i,1)-m.ElemNC(i,0);
                        X.ElemNC(i,m.Columns()-1) = m.ElemNC(i,Columns()-1)-m.ElemNC(i,Columns()-2);
                }
                
                for(int j=0; j<m.Columns(); j++)
                {
                        for(int i=1; i<m.Rows()-1; i++)
                        {
                                Y.ElemNC(i,j) = (m.ElemNC(i+1,j)-m.ElemNC(i-1,j))/2.0;
                        }
                }
                for(int j=0; j<m.Columns(); j++)
                {
                        Y.ElemNC(0,j) = m.ElemNC(1,j)-m.ElemNC(0,j);
                        Y.ElemNC(m.Rows()-1,j) = m.ElemNC(m.Rows()-1,j))-m.ElemNC(m.Rows()-2,j));
                }
        }

        
        template <class T>
        Vector<double> Del2(Vector<T>& m)
        {
                Vector<double> temp(m.Length());
                if(m.Length()>3)
                {
                        for(int i=1; i<m.Length()-1; i++)
                        {
                                temp.ElemNC(i) = (m.ElemNC(i+1) - 2*m.ElemNC(i) + m.ElemNC(i-1))/4.0;
                        }
                        temp.ElemNC(0) = 2*temp.ElemNC(1)-temp.ElemNC(2);
                        temp.ElemNC(m.Length()-1) = 2*temp.ElemNC(m.Length()-2)-temp.ElemNC(m.Length()-3);
                }
                else if(m.Length()==3)
                {
                        temp.ElemNC(1) = (m.ElemNC(2) - 2*m.ElemNC(1) + m.ElemNC(0))/4.0;
                        temp.ElemNC(0) = temp.ElemNC(1);
                        temp.ElemNC(2) = temp.ElemNC(1);
                }
                else
                {
                        temp.ElemNC(0) = 0;
                        temp.ElemNC(m.Length()-1) = 0;
                }
                return temp;
        }

        
        template <class T>
        Matrix<double> Del2(Matrix<T>& m)
        {
                Matrix<double> temp(m.Rows(), m.Columns(),0.0);
                
                Vector<double> temp2(m.Columns());
                for(int i=0; i<m.Rows(); i++)
                {
                        if(m.Columns()>3)
                        {
                                for(int j=1; j<m.Columns()-1; j++)
                                {
                                        temp2.ElemNC(j) = (m.ElemNC(i,j+1) - 2*m.ElemNC(i,j) + m.ElemNC(i,j-1))/4.0;
                                }
                                temp2.ElemNC(0) = 2*temp2.ElemNC(1)-temp2.ElemNC(2);
                                temp2.ElemNC(m.Columns()-1) = 2*temp2.ElemNC(m.Columns()-2)-temp2.ElemNC(m.Columns()-3);
                                for(int j=0; j<m.Columns(); j++)
                                {
                                        temp.ElemNC(i,j) += temp2.ElemNC(j);
                                }
                        }
                        else if(m.Columns()==3)
                        {
                                double dummy = (m.ElemNC(i,2) - 2*m.ElemNC(i,1) + m.ElemNC(i,0))/4.0;
                                temp.ElemNC(i,0) += dummy;
                                temp.ElemNC(i,1) += dummy;
                                temp.ElemNC(i,2) += dummy;
                        }
                }
                
                Vector<double> temp3(m.Rows());
                for(int j=0; j<m.Columns(); j++)
                {
                        if(m.Rows()>3)
                        {
                                for(int i=1; i<m.Rows()-1; i++)
                                {
                                        temp3.ElemNC(i) = (m.ElemNC(i+1,j) - 2*m.ElemNC(i,j) + m.ElemNC(i-1,j))/4.0;
                                }
                                temp3.ElemNC(0) = 2*temp3.ElemNC(1)-temp3.ElemNC(2);
                                temp3.ElemNC(m.Rows()-1) = 2*temp3.ElemNC(m.Rows()-2)-temp3.ElemNC(m.Rows()-3);
                                for(int i=0; i<m.Rows(); i++)
                                {
                                        temp.ElemNC(i,j) += temp3.ElemNC(i);
                                }
                        }
                        else if(m.Rows()==3)
                        {
                                double dummy = (m.ElemNC(2,j) - 2*m.ElemNC(1,j) + m.ElemNC(0,j))/4.0;
                                temp.ElemNC(0,j) += dummy;
                                temp.ElemNC(1,j) += dummy;
                                temp.ElemNC(2,j) += dummy;
                        }
                }

                return temp;
        }



        
        template <class T>
        Vector<T> CircShift(Vector<T>& m, int shift)
        {
                Vector<T> temp(m.Length());
                for(int i=0; i<m.Length(); i++)
                {
                        int ind = (i-shift) % m.Length();
                        ind = ind >=0 ? ind : ind+m.Length();
                        temp.ElemNC(i) = m.ElemNC(ind);
                }
                return temp;
        }

        
        template <class T>
        Matrix<T> CircShift(Matrix<T>& m, int rshift, int cshift)
        {
                Matrix<T> temp(m.Rows(), m.Columns());
                for(int i=0; i<m.Rows(); i++)
                {
                        int rind = (i-rshift) % m.Rows();
                        rind = rind >=0 ? rind : rind+m.Rows();
                        for(int j=0; j<m.Columns(); j++)
                        {
                                int cind = (j-cshift) % m.Columns();
                                cind = cind >=0 ? cind : cind+m.Columns();
                                temp.ElemNC(i,j) = m.ElemNC(rind,cind);
                        }
                }
                
                return temp;
        }

        
        
        template <class T>
        Vector<T> FFTShift(Vector<T>& m)
        {
                return CircShift(m, m.Length()/2);
        }

        
        template <class T>
        Matrix<T> FFTShift(Matrix<T>& m)
        {
                return CircShift(m, -m.Rows()/2, -m.Columns()/2);
        }

        template <class T>
        Vector<T> IFFTShift(Vector<T>& m)
        {
                return CircShift(m, m.Length()/2);
        }

        
        template <class T>
        Matrix<T> IFFTShift(Matrix<T>& m)
        {
                return CircShift(m, -m.Rows()/2, -m.Columns()/2);
        }
        
        
        
        template <class T>
        Vector<int> Floor(Vector<T>& m)
        {
                Vector<int> temp(m.Length());
                for(int i=0; i<m.Length(); i++)
                {
                        temp.ElemNC(i) = (int)floor(m.ElemNC(i));
                }
                return temp;
        }


        template <class T>
        Matrix<int> Floor(Matrix<T>& m)
        {
                Matrix<int> temp(m.Length());
                for(int i=0; i<m.Length(); i++)
                {
                        temp.ElemNC(i) = (int)floor(m.ElemNC(i));
                }
                return temp;
        }


        template <class T>
        Vector<int> Ceil(Vector<T>& m)
        {
                Vector<int> temp(m.Length());
                for(int i=0; i<m.Length(); i++)
                {
                        temp.ElemNC(i) = (int)ceil(m.ElemNC(i));
                }
                return temp;
        }

        
        template <class T>
        Matrix<int> Ceil(Matrix<T>& m)
        {
                Matrix<int> temp(m.Length());
                for(int i=0; i<m.Length(); i++)
                {
                        temp.ElemNC(i) = (int)ceil(m.ElemNC(i));
                }
                return temp;
        }

        
        template <class T>
        Vector<int> Fix(Vector<T>& m)
        {
                Vector<int> temp(m.Length());
                for(int i=0; i<m.Length(); i++)
                {
                        temp.ElemNC(i) = (m.ElemNC(i) >= 0) ? (int)floor(m.ElemNC(i)) : (int)ceil(m.ElemNC(i));
                }
                return temp;
        }

        
        template <class T>
        Matrix<int> Fix(Matrix<T>& m)
        {
                Matrix<int> temp(m.Length());
                for(int i=0; i<m.Length(); i++)
                {
                        temp.ElemNC(i) = (m.ElemNC(i) >= 0) ? (int)floor(m.ElemNC(i)) : (int)ceil(m.ElemNC(i));
                }
                return temp;
        }

        
        template <class T>
        Vector<int> Round(Vector<T>& m)
        {
                Vector<int> temp(m.Length());
                for(int i=0; i<m.Length(); i++)
                {
                        temp.ElemNC(i) = (int)floor(m.ElemNC(i)+0.5);
                }
                return temp;
        }

        
        template <class T>
        Matrix<int> Round(Matrix<T>& m)
        {
                Matrix<int> temp(m.Length());
                for(int i=0; i<m.Length(); i++)
                {
                        temp.ElemNC(i) = (int)floor(m.ElemNC(i)+0.5);
                }
                return temp;
        }



        
        
        
};

---