#include "RMat.h"   // For RMat class
#include <limits>   // For NaN
#include <iostream> // For cout, cerr, endl, and flush
#include <assert.h> // For assert

#if COMPILE_WITH_R
#include "R.h"      // For Rprintf
#endif

using std::cout;
using std::cerr;
using std::endl;
using std::flush;

//- Constructor #1
//
RMat::RMat() : mValues(0), mNumRows(0),
mNumCols(0), mNaN(std::numeric_limits<double>::quiet_NaN())
{
}

//- Constructor #2
//
RMat::RMat(int x, int y) : mValues(0), mNumRows(0),
mNumCols(0), mNaN(std::numeric_limits<double>::quiet_NaN())
{
    AllocateMemory(x,y);
}

//- Constructor #3
//  This constructor will be used to point an RMat object to the memory
//  location provided by R, via the '.C()' function.
//
RMat::RMat(double *vElements, int x, int y) :
mNaN(std::numeric_limits<double>::quiet_NaN())
{
    mValues  = vElements;
    mNumRows = x;
    mNumCols = y;
}

//- Constructor #4 (copy constructor)
//
RMat::RMat(const RMat &vRhs) : mValues(0), mNumRows(0),
mNumCols(0), mNaN(std::numeric_limits<double>::quiet_NaN())
{
    AllocateMemory(vRhs.NumRows(),vRhs.NumCols());
    for (int i=1; i<=mNumRows; i++)
    {
        for (int j=1; j<=mNumCols; j++)
        {
            (*this)(i,j) = vRhs(i,j);
        }
    }
}

//- Destructor
//
RMat::~RMat()
{
    DeallocateMemory();
}

//- Assignment operator #1
//
RMat&
RMat::operator= (const RMat &vRhs)
{
    DeallocateMemory();
    AllocateMemory(vRhs.NumRows(),vRhs.NumCols());
    for (int i=1; i<=mNumRows; i++)
    {
        for (int j=1; j<=mNumCols; j++)
        {
            (*this)(i,j) = vRhs(i,j);
        }
    }
    return (*this);
}

//- Assignment operator #2
//
RMat&
RMat::operator= (double vValue)
{
    for (int i=0; i<mNumRows*mNumCols; i++)
    {
        mValues[i] = vValue;
    }
    return (*this);
}

//- Multiplication operator
//
RMat
RMat::operator* (const RMat &vRhs)
{
    assert(mNumCols==vRhs.mNumRows);
    RMat vProduct(mNumRows,vRhs.NumCols());
    for (int i=1; i<=mNumRows; i++)
    {
        for (int j=1; j<=vRhs.NumCols(); j++)
        {
            for (int k=1;k<=mNumCols;k++)
            { 
                vProduct(i,j)+=(*this)(i,k)*vRhs(k,j); 
            } 
        }
    }
    return vProduct;
}

//- Element access
//
double&
RMat::operator()(int x, int y) const
{
    //- Basic range checks.
    //
    if ( ( x < 1 ) || ( x > mNumRows )
      || ( y < 1 ) || ( y > mNumCols ) )
    {
         cerr << "Rmat::operator(): range error!" << endl;
         return (double &) mNaN;
    }
    
    //- Offset is ((y-1)*mNumRows)+(x-1)
    //  rather than ((x-1)*mNumCols)+(y-1)
    //  because we want to organize memory the same
    //  way as it is in an R matrix.  R is column-major.
    //  Reference:
    //  http://en.wikipedia.org/wiki/Row-major_order#Column-major_order
    //
    return mValues[((y-1)*mNumRows)+(x-1)];
}


//- Transpose
//
RMat
RMat::Transpose()
{
    RMat vTranspose(mNumCols,mNumRows);
    for (int i=1; i<=mNumRows; i++)
    {
        for (int j=1; j<=mNumCols; j++)
        {
            vTranspose(j,i) = (*this)(i,j); 
        }
    }
    return vTranspose;
}

//- Print #1
//
void
RMat::Print()
{
    for (int i=1; i<=mNumRows; i++)
    {
        for (int j=1; j<=mNumCols; j++)
        {
            cout << (*this)(i,j) << "  ";
        }
        cout << endl;
    }
}

//- Print #2
//
void
RMat::Print(const char *vString)
{
    cout << vString;
    Print();
}

#if COMPILE_WITH_R
//- RPrint #1
//
void
RMat::RPrint()
{
    for (int i=1; i<=mNumRows; i++)
    {
        for (int j=1; j<=mNumCols; j++)
        {
            Rprintf("%g  ",(*this)(i,j));
        }
        Rprintf("\n");
        R_FlushConsole();
        R_ProcessEvents();
    }
}

//- RPrint #2
//
void
RMat::RPrint(const char *vString)
{
    Rprintf("%s",vString);
    RPrint();
}
#endif

//- Get number of rows
//
int
RMat::NumRows() const
{
    return mNumRows;
}

//- Get number of columns
//
int
RMat::NumCols() const
{
    return mNumCols;
}

//- Allocate memory
//  Initialize all entries to 0.
//
bool             
RMat::AllocateMemory(int x, int y)
{
    DeallocateMemory();
    try {
        mValues  = new double [x*y];
    }
    catch (...) {
        cerr << "Rmat::AllocateMemory(): unable to allocate memory for "
             << x << " by " << y << " matrix!" << endl;
        return false;
    }

    mNumRows = x;
    mNumCols = y;
    for (int i=1; i<=mNumRows; i++)
    {
        for (int j=1; j<=mNumCols; j++)
        {
            (*this)(i,j) = 0.0;
        }
    }
    return true;
}

//- Deallocate memory
//
void             
RMat::DeallocateMemory()
{
    if ( mValues != 0 )
    {
        delete [] mValues;
        mValues = 0;
    }
}