Shark Machine Learning Library

• About Shark
• Sourceforge
  • Project Summary
  • Downloads
  • Subversion Repository
• Getting Started
• Tutorials
• FAQ
• Main Modules
  • ReClaM
  • EALib
  • MOO-EALib
  • Fuzzy
• Tools
  • Mixture
  • Array
  • Rng
  • LinAlg
  • FileUtil

LinAlg.cpp

Go to the documentation of this file.

//===========================================================================
//======================================================================


#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <SharkDefs.h>
#include <Array/ArrayOp.h>
#include <LinAlg/LinAlg.h>


//===========================================================================
Array< double > mean(const Array< double >& x)
{
    SIZE_CHECK(x.ndim() > 1)

    Array< double > m(x[ 0 ]);
    for (unsigned i = 1; i < x.dim(0); ++i)
        m += x[ i ];
    m /= double(x.dim(0));
    return m;
}



//===========================================================================
Array< double > variance(const Array< double >& x)
{
    SIZE_CHECK(x.ndim() > 1 && x.dim(0) > 0)

    unsigned i, j;
    unsigned size = x.nelem() / x.dim(0);
    Array< double > m(size);   // vector of mean values.
    Array< double > v(size);   // vector of variance values

    m = 0.;
    v = 0.;

    for (i = 0; i < x.dim(0); ++i)
        for (j = 0; j < size; ++j) {
            m(j) += x.elem(i * size + j);// sum of elements in j's row
            v(j) += Shark::sqr(x.elem(i * size + j));//sum of squares od elements in j's row
        }

    for (j = 0; j < size; ++j) {
        m(j) /= x.dim(0);//normalize mean
        v(j)  = Shark::max(0., v(j) / x.dim(0) - m(j) * m(j));
    }

    // re-arrange variance array
    if (x.ndim() > 1)
        v.resize(x[ 0 ], true);//ensure v is a column vector

    return v;
}



//===========================================================================
double angle(const Array< double >& x, const Array< double >& y)
{
    //return scalar_product( x, y ) / sqrt( sumOfSqr( x ) * sumOfSqr( y ) );

    SIZE_CHECK(x.ndim() == 1 && x.samedim(y))

    double sumxx = 0;
    double sumyy = 0;
    double sumxy = 0;

    for (unsigned i = 0; i < x.nelem(); ++i) {
        sumxx += x(i) * x(i);
        sumyy += y(i) * y(i);
        sumxy += x(i) * y(i);
    }

    return (sumxx == 0 || sumyy == 0) ? 0. : sumxy / sqrt(sumxx * sumyy);
}


//===========================================================================
void meanvar
(
    const Array< double >& x,
    Array< double >&       m,
    Array< double >&       v
)
{
    SIZE_CHECK(x.ndim() > 1)

    unsigned i, j;

    // number of elements for each "row":
    unsigned size = x.nelem() / x.dim(0);

    m.resize(size, false);
    v.resize(size, false);
    m = 0.;
    v = 0.;

    for (i = 0; i < x.dim(0); ++i) {
        for (j = 0; j < size; ++j) {
            // Sum of elements of each "row":
            m(j) += x.elem(i * size + j);
            v(j) += Shark::sqr(x.elem(i * size + j));
        }
    }

    for (j = 0; j < size; ++j) {
        // Calculate mean value by dividing
        // by no. of "rows":
        m(j) /= x.dim(0);
        v(j)  = v(j) / x.dim(0) - m(j) * m(j);
    }

    // re-arrange arrays
    m.resize(x[ 0 ], true);
    v.resize(x[ 0 ], true);
}

/*
    /brief Calculates the mean and variance values of 1d-arrays p(x)

    calculates the first two moments of a discrete value array and
    his corresponding x-values
    /param pxA      1d array with functionvalues
    /param xA       1d array with corresponding x values
    /param mA       retunsvalue meanval
    /param vA       returnvalue variance
    /param startA   start of a calculation window
    /param endA     end of a calculation window
*/
void meanvar
(
    const Array< double >& pxA,
    const Array< double >& xA,
    double &mA,
    double &vA,
    const int startA,
    const int endA
)
{
//  SIZE_CHECK( xA.ndim( ) != 1 )
//  SIZE_CHECK( pxA.ndim( ) != xA.ndim( ) )
//  SIZE_CHECK( pxA.nelem( ) != xA.nelem( ) )

    int size = pxA.nelem();
    int startL  = (startA < 0) ? 0 : startA;
    int endL    = (endA == -1) ? size : endA;
    if (startL == endL) {
        mA = startA;
        vA = 0;
    }
    double ew = 0;
    double sum = 0;
    int i;
    // calculate first moment
    for (i = startL; i < endL ;++i) {
        ew  += pxA.elem(i) * xA.elem(i);
        sum += pxA.elem(i);
    }
    mA = (ew == 0) ? 0 : ew / sum;
    // calculate variance
    // variance is the second central moment
    const int ii = 2;
    double vc = 0;
    for (i = startL; i < endL ;++i) {
        vc += pow((xA.elem(i) - mA), (double)ii) * pxA.elem(i);
    }
    vA = sqrt((vc == 0) ? 0 : vc / sum);

}




//===========================================================================
double corrcoef(const Array< double >& x, const Array< double >& y)
{
    SIZE_CHECK(x.ndim() == 1 && x.samedim(y))

    double sumxx = 0;
    double sumyy = 0;
    double sumxy = 0;
    double sumx  = 0;
    double sumy  = 0;

    if (x.nelem() > 0) {
        for (unsigned i = 0; i < x.nelem(); ++i) {
            sumxx += x(i) * x(i);
            sumyy += y(i) * y(i);
            sumxy += x(i) * y(i);
            sumx  += x(i);
            sumy  += y(i);
        }

        sumxx /= x.nelem();
        sumyy /= x.nelem();
        sumxy /= x.nelem();
        sumx  /= x.nelem();
        sumy  /= x.nelem();

        sumxx -= sumx * sumx;
        sumyy -= sumy * sumy;
        sumxy -= sumx * sumy;
    }

    return (sumxx == 0 || sumyy == 0) ? 0. : sumxy / sqrt(sumxx * sumyy);
}



//===========================================================================
Array< double > corrcoef(const Array< double >& x)
{
    SIZE_CHECK(x.ndim() > 1)

    unsigned i, j;
    Array< double > C(covariance(x));

    for (i = 0; i < C.dim(0); ++i)
        for (j = 0; j < i; ++j)
            if (C(i, i) == 0 || C(j, j) == 0)
                C(i, j) = C(j, i) = 0;
            else
                C(i, j) = C(j , i) = C(i, j) / sqrt(C(i, i) * C(j, j));

    for (i = 0; i < C.dim(0); ++i)
        C(i, i) = 1;

    return C;
}


//===========================================================================
double covariance(const Array< double >& x, const Array< double >& y)
{
    SIZE_CHECK(x.ndim() == 1 && x.samedim(y))

    double sumxy = 0;
    double sumx  = 0;
    double sumy  = 0;

    if (x.nelem() > 0) {
        for (unsigned i = 0; i < x.nelem(); ++i) {
            sumxy += x(i) * y(i);
            sumx  += x(i);
            sumy  += y(i);
        }

        sumxy /= x.nelem();
        sumx  /= x.nelem();
        sumy  /= x.nelem();
    }

    return sumxy - sumx * sumy;
}

//===========================================================================
Array< double > covariance(const Array< double >& x)
{
    SIZE_CHECK(x.ndim() == 2)

    unsigned num = x.dim(0);
    unsigned dim = x.dim(1);

    Array< double > mean(dim);
    Array< double > covar(dim, dim);

    mean  = 0.;
    covar = 0.;

    for (unsigned i = num; i--;) {
        mean  += x[ i ];
        covar += outerProduct(x[ i ], x[ i ]);
    }

    mean  /= double(num);
    covar /= double(num);
    covar -= outerProduct(mean, mean);

    return covar;
}

//===========================================================================
void matMat(Array2D<double> &A, const Array2D<double> &B, const Array2D<double> &C) {
  unsigned i, j, k, n = B.dim(0), m = C.dim(1), l = C.dim(0);
  SIZE_CHECK( B.dim(1) == C.dim(0) );  
  A.resize(n, m, false);
  for(i=0; i<n; i++) {
    for(j=0; j<m; j++) {
      A(i,j)=B(i, 0) * C(0, j);
      for(k=1; k<l; k++) {
                A(i,j)+=B(i, k) * C(k, j);
      }
    }
  }
}

//===========================================================================
void matColVec(Array<double> &A, const Array2D<double> &B, const Array<double> &C)
{
    SIZE_CHECK(B.ndim() == 2);
    SIZE_CHECK(C.ndim() == 1);

    unsigned int v, vc = B.dim(0);
    unsigned int h, hc = B.dim(1);

    SIZE_CHECK(C.dim(0) == hc);
    A.resize(vc, false);

    for (v=0; v<vc; v++)
    {
        double value = 0.0;
        for (h=0; h<hc; h++) value += B(v, h) * C(h);
        A(v) = value;
    }
}

//===========================================================================
void matColVec(ArrayReference<double> A, const Array2D<double>& B, const ArrayReference<double> C)
{
    SIZE_CHECK(B.ndim() == 2);
    SIZE_CHECK(C.ndim() == 1);

    unsigned int v, vc = B.dim(0);
    unsigned int h, hc = B.dim(1);

    SIZE_CHECK(C.dim(0) == hc);
    A.resize(vc, false);

    for (v=0; v<vc; v++)
    {
        double value = 0.0;
        for (h=0; h<hc; h++) value += B(v, h) * C(h);
        A(v) = value;
    }
}

//===========================================================================
void matColVec(Array<double> &A, const Array2D<double> &B, const Array<double> &C, unsigned int index)
{
    unsigned i, j, n, m;
    SIZE_CHECK( B.ndim() == 2 );
    SIZE_CHECK( C.ndim() == 2 );
    RANGE_CHECK(index < C.dim(1));
    SIZE_CHECK( B.dim(1) == C.dim(0) );
    n = B.dim(0);
    m = B.dim(1);
    A.resize(n, false);
    for(i=0; i<n; i++) {
        A(i)=B(i, 0) * C(0, index);
        for(j=1; j<m; j++) {
            A(i)+=B(i, j) * C(j, index);
        }
    }
}

//===========================================================================
double vecMatVec(const Array<double> &A, const Array2D<double> &B, const Array<double> &C) {
  unsigned i, j, n, m;
  double sum =  0, help;
  SIZE_CHECK( A.ndim() == 1 );  
  SIZE_CHECK( B.ndim() == 2 );  
  SIZE_CHECK( C.ndim() == 1 );  
  SIZE_CHECK( B.dim(1) == C.dim(0) );  
  SIZE_CHECK( B.dim(0) == A.dim(0) );  
  n = B.dim(0);
  m = B.dim(1);  

  for(i=0; i<n; i++) {
    help=B(i, 0) * C(0);
    for(j=1; j<m; j++) {
    help+=B(i, j) * C(j);
    }
    sum+=A(i)*help;
  }
  return sum;
}

//===========================================================================
double vecMatVec(const Array<double> &A, unsigned int i, const Array2D<double> &B, const Array<double> &C, unsigned int j) {
    unsigned a, c, n, m;
    double sum =  0, help;
    SIZE_CHECK( A.ndim() == 2 );
    SIZE_CHECK( B.ndim() == 2 );
    SIZE_CHECK( C.ndim() == 2 );
    RANGE_CHECK(i < A.dim(1));
    RANGE_CHECK(j < C.dim(1));
    SIZE_CHECK( B.dim(1) == C.dim(0) );
    SIZE_CHECK( B.dim(0) == A.dim(0) );
    n = B.dim(0);
    m = B.dim(1);

    for(a=0; a<n; a++) {
        help=B(a, 0) * C(0, j);
        for(c=1; c<m; c++) {
            help+=B(a, c) * C(c, j);
        }
        sum+=A(a, i)*help;
    }
    return sum;
}

//===========================================================================
void invertSymm( Array2D<double> &I, const Array2D< double >& A) {
  SIZE_CHECK( A.ndim() == 2 );  
  SIZE_CHECK( A.dim(1) == A.dim(0) );  

  Array<double> lambda;
  Array2D<double> eigenVectors;
  unsigned n = A.dim(0);
  
  I.resize(A, false);
  lambda.resize(n, false);
  eigenVectors.resize(A, false);
  lambda = 0;
  eigenVectors = 0.;
  
  eigensymm(A, eigenVectors, lambda);
  I = 0.;
  
  for(unsigned i = 0; i < n; ++i) I(i,i) = 1./lambda(i);

  I = innerProduct(I , transpose(eigenVectors));
  I = innerProduct(eigenVectors , I);
}


//===========================================================================
void CholeskyDecomposition(const Array2D< double >& M,
                                                     Array2D< double >& C) 
{
    int i, j, k, m;
    double s;
    SIZE_CHECK(M.dim(0) == M.dim(1));
    
    m = M.dim(0);
    C.resize(m, m, false);
    for(j = 0; j < m; j++) {
        for(i = j; i < m; i++) {
            s = M(i, j);
            for(k = 0; k < j; k++) {
                s -= C(i, k) * C(j, k);
            }
            if (i == j) {
                C(i, j) = sqrt(s);
            }
            else {
                C(i, j) = s/C(j , j);
                C(j, i) = 0;
            }
        }
    }
}