linalg.cpp

#include "linalg.h"
#include "util.h"
#include "qscr.h"
#include <sstream>
#include <iomanip>
#include <cassert>
#include <cstring>
#include <cstdlib>
#include <cmath>


matrix matrix::addcolumns()
{
  long i,j;
  matrix B(rows);
  for(i=0;i<=rows-1;i++)
    {   
      for(j=0;j<=columns-1;j++)
        B.A[i][0] += A[i][j];
    }
  return(B);
}

matrix matrix::getrow(long i) const
{
  long j;
  matrix B(1,columns);
  for (j=0; j<=columns-1;j++)
    B.A[0][j] = A[i][j];
  return(B);
}

matrix matrix::getcolumn(long i) const
{
  long j;
  matrix B(rows);
  for (j=0; j<=rows-1; j++)
    B.A[j][0] = A[j][i];
  return(B);
}

void matrix::setrow(const matrix &B, long i)
{
  long j;
  if (columns != B.columns)
    error("Wrong dimensions in matrix setrow \n");
  for (j=0; j<=columns-1; j++)
    A[i][j] = B.A[0][j];
}
  
void matrix::setcolumn(const matrix &B, long i)
{
  long j;
  if (rows != B.rows)
    error("Wrong dimensions in matrix setcolumn \n");
  for (j=0; j<=rows-1; j++)
    A[j][i] = B.A[j][0];
}

long matrix::operator==(const matrix &B)
{
  long i,j,eq=1;
  if ((rows != B.rows) & (columns != B.columns))
    error("Wrong dimensions in matrix equal \n");                           
  else 
    {
      for(i=0; (i<=rows-1) & eq; i++)
        {
          for(j=0; (j<=columns-1) & eq; j++)
            eq = (A[i][j] == B.A[i][j]);
        }
    }
  return(eq);
}
        
void matrix::print() const
{
  long i,j;
  //  std::cerr << "\n";
  for(i=0;i<=rows-1;i++)
    {
      for(j=0;j<=columns-1;j++) {
        std::string ss = numStr(A[i][j]) + "    ";
        qaddstr(ss.c_str());    
      }
        qinsertln();
        //std::cerr << A[i][j] << " ";
        //      std::cerr << "\n";
    }
}

void matrix::stdprint() const
{
  long i,j;
  std::cerr  << "[";
  for(i=0;i<=rows-1;i++)
    {
      for(j=0;j<=columns-1;j++) {
        std::cerr << A[i][j];
        if (j<columns-1)
          std::cerr << " ";
      }
      if (i<rows-1)
        std::cerr << ";" << std::endl;
    }
  std::cerr << "]";
}



matrix matrix::operator*(const double &x)
{
  long i,j;
  matrix C(rows,columns);
  for(i=0;i <= C.rows-1; i++)
    {
      for(j=0; j<=C.columns-1;j++)
        {
          C.A[i][j] = x*A[i][j];
        }
    }
  return(C);
}


matrix matrix::operator*(const matrix &B)
{
  long i,j,k;
  matrix C(rows,B.columns);
  if (columns != B.rows) 
    error("Wrong dimensions in matrix multiplication \n");
  else 
    { 
      for(i=0;i <= C.rows-1; i++)
        {
          for(j=0; j<=C.columns-1;j++)
            {
              for(k=0;k<=columns-1;k++)
                C.A[i][j] += A[i][k]*B.A[k][j];
            }
        }
    }
  return(C);
}

matrix matrix::trans() 
{
  long i,j;
  matrix B(columns,rows);
  for (i=0;i<=columns-1;i++)
    {
      for(j=0;j<=rows-1;j++)
        B.A[i][j] = A[j][i];
    }
  return(B);
}
    
matrix matrix::operator+(const matrix &B)
{
  long i,j;
  matrix C(rows,columns);
  if ((rows != B.rows) & (columns != B.columns))
    error("Wrong dimensions in matrix plus \n");                            
  else 
    {
      for(i=0;i <= rows-1; i++)
        {
          for(j=0; j<=columns-1;j++)
            C.A[i][j] = A[i][j] + B.A[i][j];
        }
    }
  return(C);
}

matrix matrix::operator-(const matrix &B)
{
  long i,j;
  matrix C(rows,columns);
  if ((rows != B.rows) & (columns != B.columns))
    error("Wrong dimensions in matrix minus \n");                           
  else 
    {
      for(i=0;i <= rows-1; i++)
        {
          for(j=0; j<=columns-1;j++)
            C.A[i][j] = A[i][j] - B.A[i][j];
        }
    }
  return(C);
}

double matrix::twonorm()
{
  long i,j;
  double norm=0;
  for(i=0;i<=rows-1;i++)
    {
      for(j=0;j<=columns-1;j++)
        norm += A[i][j]*A[i][j];
    }
  norm = sqrt(norm);
  return(norm);
}

double matrix::onenorm()
{
  long i,j;
  double norm=0;
  for(i=0;i<=rows-1;i++)
    {
      for(j=0;j<=columns-1;j++)
        norm += num(A[i][j]);
    }
  return(norm);
}

safevector<double> matrix::mean() const {
  safevector<double> m(columns);
  for (long s=0; s<columns; s++) {
    double sum = 0;
    for (long r=0; r<rows; r++) {
      sum += A[r][s];
    }
    m[s] = sum/(double)rows;
  }
  return m;
}

safevector<double> matrix::var() const {
  safevector<double> mu = mean();
  safevector<double> v(columns);
  for (long s=0; s<columns; s++) {
    double sum = 0;
    for (long r=0; r<rows; r++) {
      double xe = A[r][s]-mu[s];
      sum +=  xe*xe;
    }
    v[s] = sum/((double)rows-1);
  }
  return v;
}

matrix matrix::covar() const {
  safevector <double> m = mean();
  matrix res(columns, columns);
  for (int x=0; x<columns; x++) {
    for (int y=x; y<columns; y++) {
      double sum = 0;
      for (int i=0; i<rows; i++) {
        sum += (A[i][x]-m[x])*(A[i][y]-m[y]);
      }
      double val = sum/(double)(rows-1);
      res.set(val,x,y);
      res.set(val,y,x);
    }
  }
  return res;
}


matrix matrix::inv()
{
  long *indxc =new long[rows];
  long *indxr = new long[rows];
  long *ipiv = new long[rows];
  long i, icol = 0, irow = 0, j, k, l, ll;
  double big, dum, pivinv;
  matrix B(rows,columns);

  for(i=0;i<=rows-1;i++)
    {
      indxc[i]=0;
      indxr[i]=0;
      ipiv[i]=0;
    }

  if (columns != rows) 
    error("Cannot invert non-square matrix \n");
  else
    for(i=0;i<=rows-1;i++)
      {
        for(j=0;j<=columns-1;j++)
          B.A[i][j] = A[i][j];
      }
 
      
  for (i=0;i<=columns-1;i++)
    {
      big = 0.0;
      for (j=0;j<=rows-1;j++)
        {
          if (ipiv[j] != 1)
            for (k=0;k<=columns-1;k++)
              {
                if (ipiv[k] == 0)
                  {
                    if (fabs(B.A[j][k]) >=big)
                      {
                        big = fabs(B.A[j][k]);
                        irow = j;
                        icol = k;
                      }
                  } 
                else if (ipiv[k] > 1) error("inv: Singular Matrix-1");
              }
        }
      ++(ipiv[icol]);
      if (irow != icol)
        {
          for (l=0;l<=columns-1;l++) B.swap(irow,icol);
        }
      indxr[i] = irow;
      indxc[i] = icol;
      if (B.A[icol][icol] == 0.0) error("inv: Singular Matrix-2");
      pivinv = 1.0/B.A[icol][icol]; 
      B.A[icol][icol] = 1.0;
      for (l=0;l<=columns-1;l++) B.A[icol][l] *= pivinv;
      for (ll=0;ll<=rows-1;ll++)
        if (ll != icol)
          {
            dum = B.A[ll][icol];
            B.A[ll][icol] = 0.0;
            for (l=0;l<=columns-1;l++) 
              B.A[ll][l] -= B.A[icol][l]*dum;
      
          }
    }
  for(l=columns-1;l>=0;l--)
    {
      if (indxr[l] != indxc[l])
        B.swapc(indxr[l],indxc[l]);
    }
  delete(ipiv);
  delete(indxr);
  delete(indxc);
  return(B);
}

matrix matrix::map(double(*f)(double)) {
  matrix res(rows,columns);
  for (int r=0; r<rows; r++)
    for (int s=0; s<columns; s++)
      res.set(f(A[r][s]),r,s);

  return res;
}

double matrix::maxvalue(){

  double max=A[0][0];  
  for(int i=0; i<rows;i++){
    for(int j=0; j<columns; j++)
      if (A[i][j]>max) max=A[i][j];
  }
  return(max);

}



const char* matrix::row2string(int r) {
  assert(r>=0 && r<rows);

  std::ostringstream mstr;
  mstr.setf(std::ios::fixed);
  std::string res;
  for (int i=0; i<columns; i++) {
    mstr << std::setw(9) << A[r][i] << " ";
  }
  
  std::string ss = mstr.str();    

  return ss.c_str();
}


matrix covar2(matrix const& X, matrix const& Y) {
  assert(X.row()==Y.row() && (X.column()==Y.column() || Y.column()==1));

  safevector <double> mx = X.mean();
  safevector <double> my = Y.mean();
  matrix cv(X.column());


  // If Y is a column-vector we calculate [cov(X_1,Y),...,cov(X_n,Y)].
  if (Y.column()==1) {
    for (int s=0; s<X.column(); s++) {
      double sum=0;
      for (int i=0; i<X.row(); i++) {
        sum += (X.get(i,s)-mx[s])*(Y.get(i)-my[0]);
        //      std::cerr << "sum+=" << sum << std::endl;
      }
      cv.set(sum/(double)(X.row()-1),s);
    }
    return cv; 
  }
 

  // Else we calculate [cov(X_1,Y_1),....,cov(X_n,Y_n)].
  for (int s=0; s<X.column(); s++) {
    double sum=0;
    for (int i=0; i<X.row(); i++) {
      sum += (X.get(i,s)-mx[s])*(Y.get(i,s)-my[s]);
    }
    cv.set(sum/(double)(X.row()-1),s);
  }
  return cv;
}


//####################################################################

#ifdef TEST


double cubic(double x) {
  return x*x*x;
}

double square(double x) {
  return x*x;
}

double identity(double x) {
  return x;
}

int main() {
  try {

    matrix A(4,4);
    for (unsigned i=0; i<4; i++)
      for (unsigned j=0; j<4; j++)
        A.set(i+j, i,j);
    
    std::cerr << "A= "; A.stdprint();
    matrix B = A.map(cubic);
    std::cerr << "cubic=\n "; B.stdprint();
    B = A.map(square);
    std::cerr << "square=\n "; B.stdprint();
    B = A.map(identity);
    std::cerr << "identity=\n "; B.stdprint();
    

    matrix C(4,4);
    C.set(3,0,0);    C.set(2,0,1);    C.set(1,0,2);    C.set(5,0,3);
    C.set(4,1,0);    C.set(2,1,1);    C.set(1,1,2);    C.set(5,1,3);
    C.set(3,2,0);    C.set(4,2,1);    C.set(1,2,2);    C.set(4,2,3);
    C.set(1,3,0);    C.set(3,3,1);    C.set(1,3,2);    C.set(3,3,3);

    matrix vec(4);
    vec.set(1,0);    vec.set(4,1);    vec.set(2,2);    vec.set(9,3);
    std::cerr << "C="; C.stdprint(); 
    std::cerr << "x="; vec.stdprint(); 
    matrix covCx = covar2(C,vec);
    std::cerr << "cov(C,x)= "; covCx.stdprint();


    std::cerr << "Nothing to test in linalg.cpp\n";
  } catch(std::string const& err) {
    std::cerr << err << std::endl;
    return 1;
  }
  return 0;
}


#endif//TEST

---