/*
 * Mtx.cpp
 *
 *  Created on: Feb 15, 2012
 *      Author: godin
 */

#include "data/Mtx.h"
#include "io/MtxIOStream.h"
#include "utils/AuvMath.h"

Mtx::Mtx( Logger& logger, int m, int n, float initVal, Syslog::Severity severity )
    : m_( m ),
      n_( n ),
      value_( NULL ),
      logger_( logger ),
      severity_( severity ),
      unlimitedN_( false )
{
    setSize( m, n );
    if( m * n > 0 )
    {
        for( int i = 0; i < m; ++i )
        {
            for( int j = 0; j < n; ++j )
            {
                value_[i][j] = initVal;
            }
        }
    }
}

Mtx::Mtx( Logger& logger, int m, int n, float* ptr, Syslog::Severity severity )
    : m_( m ),
      n_( n ),
      value_( NULL ),
      logger_( logger ),
      severity_( severity ),
      unlimitedN_( false )
{
    setSize( m, n );
    memcpy( value_, ptr, m * n * sizeof( float ) ); // Don't usually like working at the level of memcpy, but it seems the most straightforward way to do this.
}

Mtx::Mtx( const char* filename, Logger& logger, Syslog::Severity severity )
    : m_( 0 ),
      n_( 0 ),
      value_( NULL ),
      logger_( logger ),
      severity_( severity ),
      unlimitedN_( false )
{
    MtxIOStream stream( filename, false );
    if( stream.isReadable() )
    {
        *this = stream.read( logger, 0, 0, Syslog::FAULT );
    }
}

Mtx::Mtx( const char* filename, Logger& logger, int minN, int maxN, Syslog::Severity severity )
    : m_( 0 ),
      n_( 0 ),
      value_( NULL ),
      logger_( logger ),
      severity_( severity ),
      unlimitedN_( false )
{
    MtxIOStream stream( filename, false );
    if( stream.isReadable() )
    {
        *this = stream.read( logger, minN, maxN, severity );
    }
}

Mtx::Mtx( const Mtx& rhs )
    : m_( 0 ),
      n_( 0 ),
      value_( NULL ),
      logger_( rhs.logger_ ),
      severity_( rhs.severity_ )
{
    *this = rhs;
}

Mtx::~Mtx()
{
    clear();
}

void Mtx::setSize( int m, int n )
{
    clear();
    value_ = New2D( m, n, true );
    m_ = m;
    n_ = n;
    if( NULL != value_ )
    {
        value_[m_][0] = 1;
    }
    /*value_ = m * n > 0 ? new float*[m + 1] : NULL;
    if( m * n > 0 )
    {
        value_[0] = new float[m_ * n_ + 1];
        for( int i = 1; i <= m_; ++i )
        {
            value_[i] = value_[i - 1] + n_;
        }
        value_[m_][0] = 1;
    }*/
}

Mtx& Mtx::operator=( const Mtx & rhs )
{
    if( m_ == rhs.m_ && n_ == rhs.n_ && m_ * n_ > 0 && value_[m_][0] == 1.0f )
    {
        memcpy( value_[0], rhs.value_[0], sizeof( float )*m_ * n_ );
    }
    else
    {
        clear();
        m_ = rhs.m_;
        n_ = rhs.n_;
        if( m_ * n_ > 0 )
        {
            value_ = rhs.value_;
            value_[m_][0] = value_[m_][0] + 1;
        }
    }
    return *this;
}

Mtx& Mtx::operator*=( const float rhs )
{
    prepareValueForModification();
    for( int i = 0; i < m_; ++i )
    {
        for( int j = 0; j < n_; ++j )
        {
            value_[i][j] *= rhs;
        }
    }
    return *this;
}

Mtx& Mtx::operator+=( const Mtx & rhs )
{
    prepareValueForModification();
    if( rhs.m_ > m_ || rhs.n_ > n_ )
    {
        logger_.syslog( Str( "Attempting [" ) + m_ + "*" + n_ + "]+=[" + rhs.m_ + "*" + rhs.n_ + "], truncating", severity_ );
    }
    for( int i = 0; i < m_ && i < rhs.m_; ++i )
    {
        for( int j = 0; j < n_ && j < rhs.n_; ++j )
        {
            value_[i][j] += rhs.value_[i][j];
        }
    }
    return *this;
}

Mtx& Mtx::operator+=( const float rhs )
{
    prepareValueForModification();
    for( int i = 0; i < m_; ++i )
    {
        for( int j = 0; j < n_; ++j )
        {
            value_[i][j] += rhs;
        }
    }
    return *this;
}

Mtx& Mtx::operator-=( const Mtx & rhs )
{
    prepareValueForModification();
    if( rhs.m_ > m_ || rhs.n_ > n_ )
    {
        logger_.syslog( Str( "Attempting [" ) + m_ + "*" + n_ + "]-=[" + rhs.m_ + "*" + rhs.n_ + "], truncating.", severity_ );
    }
    for( int i = 0; i < m_ && i < rhs.m_; ++i )
    {
        for( int j = 0; j < n_ && j < rhs.n_; ++j )
        {
            value_[i][j] -= rhs.value_[i][j];
        }
    }
    return *this;
}

Mtx Mtx::subset( int minM, int maxM, int minN, int maxN ) const
{
    minM = minM >= 0 ? minM : m_ + minM;
    maxM = maxM >= 0 ? maxM : m_ + maxM;
    minN = minN >= 0 ? minN : n_ + minN;
    maxN = maxN >= 0 ? maxN : n_ + maxN;
    int m = maxM - minM + 1;
    int n = maxN - minN + 1;
    Mtx mtx( logger_, m, n, nanf( "" ), severity_ );
    for( int i = 0; i < m; ++i )
    {
        for( int j = 0; j < n; ++j )
        {
            mtx.value_[i][j] = value_[i + minM][j + minN];
        }
    }
    return mtx;
}

bool Mtx::write( const char* filename, Logger& logger,
                 int minM, int maxM, int minN, int maxN,
                 Syslog::Severity severity )
{
    MtxIOStream stream( filename, true );
    if( stream.isWritable() )
    {
        return stream.write( *this, logger, minM, maxM, minN, maxN, severity );
    }
    return false;
}

bool Mtx::append( const char* filename, Logger& logger,
                  Syslog::Severity severity )
{
    MtxIOStream stream( filename, false );
    if( stream.isWritable() && stream.isReadable() )
    {
        return stream.append( *this, logger, severity );
    }
    return false;
}

void Mtx::clear()
{
    if( NULL != value_ )
    {
        if( value_[m_][0] > 1.0f )
        {
            value_[m_][0] = value_[m_][0] - 1.0f;
        }
        else
        {
            Delete2D( value_ );
        }
        value_ = NULL;
    }
    m_ = 0;
    n_ = 0;
}

void Mtx::prepareValueForModification()
{
    if( m_ * n_ > 0 )
    {
        if( value_[m_][0] > 1.0f )
        {
            value_[m_][0] = value_[m_][0] - 1.0f;
            float** oldValue = value_;
            value_ = NULL;
            setSize( m_, n_ );
            memcpy( value_[0], oldValue[0], sizeof( float )*m_ * n_ );
        }
    }
}

Mtx Mtx::operator*( const Mtx& rhs ) const
{
    Mtx retVal( logger_, 0, 0, nanf( "" ), severity_ );
    if( rhs.m_ != n_ )
    {
        logger_.syslog( Str( "Attempting [a*" ) + n_ + "]x[" + rhs.m_ + "*b], returning empty", severity_ );
        return retVal;
    }
    retVal.setSize( m_, rhs.n_ );
    for( int i = 0; i < m_; ++i )
    {
        for( int j = 0; j < rhs.n_; ++j )
        {
            float value( 0.0f );
            for( int k = 0; k < n_; ++k )
            {
                value += value_[i][k] * rhs.value_[k][j];
            }
            retVal.value_[i][j] = value;
        }
    }
    return retVal;
}

/// Simple sort routine, invented by Donald Shell 1959.
/// Returns true if the data is modified.
/// Improved with gap sequence by Ciura. It is better than O(n^[4/3]).
/// See Ciura, Marcin (2001). "Best Increments for the Average Case of Shellsort".
///  In Freiwalds, Rusins. Proceedings of the 13th International Symposium on
///  Fundamentals of Computation Theory. London: Springer-Verlag. pp. 106–117.
///  ISBN 3-540-42487-3.
bool Mtx::shellSort( int ( &comp )( const Mtx&, int, const Mtx&, int ) )
{
    Mtx temp( logger_, m_, 1, nanf( "" ), severity_ );

    const int ngaps = 12;
    const int gaps[ngaps] = { 1750, 701, 301, 132, 57, 23, 10, 4, 1};
    bool modified = false;
    for( int g = 0; g < ngaps; ++g )
    {
        int gap = gaps[g];
        for( int i = gap; i < n_; ++i )
        {
            CopyValues( *this, 0, i, temp, 0, 0, m_, 1 );
            int j = i;
            for( ; j >= gap && comp( *this, j - gap, temp, 0 ) > 0; j -= gap )
            {
                CopyValues( *this, 0, j - gap, *this, 0, j, m_, 1 );
                modified = true;
            }
            if( j != i )
            {
                CopyValues( temp, 0, 0, *this, 0, j, m_, 1 );
            }
        }
    }
    return modified;
}

bool Mtx::CopyValues( const Mtx& src, int srcStartM, int srcStartN, const Mtx& dst, int dstStartM, int dstStartN, int numM, int numN )
{
    if( src.m_ < srcStartM + numM || src.n_ < srcStartN + numN
            || dst.m_ < dstStartM + numM || dst.n_ < dstStartN + numN )
    {
        return false;
    }

    for( int i = 0; i < numM; ++i )
    {
        float* srcVal = src.value_[i + srcStartM];
        float* dstVal = dst.value_[i + dstStartM];
        for( int j = 0; j < numN; ++j )
        {
            dstVal[j + dstStartN] = srcVal[j + srcStartN];
        }
    }
    return true;
}

/// Comparison function for indexM = 0;
int Mtx::Comp0( const Mtx& lhs, int lhsJ, const Mtx& rhs, int rhsJ )
{
    return AuvMath::Sign( lhs( 0, lhsJ ) - rhs( 0, rhsJ ) );
}

/// Comparison function for indexM = 1;
int Mtx::Comp1( const Mtx& lhs, int lhsJ, const Mtx& rhs, int rhsJ )
{
    return AuvMath::Sign( lhs( 1, lhsJ ) - rhs( 1, rhsJ ) );
}

/// Comparison function for indexM = 2;
int Mtx::Comp2( const Mtx& lhs, int lhsJ, const Mtx& rhs, int rhsJ )
{
    return AuvMath::Sign( lhs( 2, lhsJ ) - rhs( 2, rhsJ ) );
}

/// Comparison function for indexM = 0 AND indexM = 1;
int Mtx::Comp01( const Mtx& lhs, int lhsJ, const Mtx& rhs, int rhsJ )
{
    int retVal = Mtx::Comp0( lhs, lhsJ, rhs, rhsJ );
    if( retVal == 0 )
    {
        retVal = Mtx::Comp1( lhs, lhsJ, rhs, rhsJ );
    }
    return retVal;
}

/// Comparison function for indexM = 0 AND indexM = 2;
int Mtx::Comp02( const Mtx& lhs, int lhsJ, const Mtx& rhs, int rhsJ )
{
    int retVal = Mtx::Comp0( lhs, lhsJ, rhs, rhsJ );
    if( retVal == 0 )
    {
        retVal = Mtx::Comp2( lhs, lhsJ, rhs, rhsJ );
    }
    return retVal;
}