/**
\file     cmatrix.c
\brief    'c' functions for vector and matrix operations.
\author   Glenn D. MacGougan (GDM)
\date     2008-09-22
\version  0.06 Beta

\b Version \b Information \n
This is the open source version (BSD license). The Professional Version
is avaiable via http://www.zenautics.com. The Professional Version
is highly optimized using SIMD and includes optimization for multi-core 
processors.

\b License \b Information \n
Copyright (c) 2008, Glenn D. MacGougan, Zenautics Technologies Inc. \n

Redistribution and use in source and binary forms, with or without
modification, of the specified files is permitted provided the following 
conditions are met: \n

- Redistributions of source code must retain the above copyright
  notice, this list of conditions and the following disclaimer. \n
- 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. \n
- The name(s) of the contributor(s) may not be used to endorse or promote 
  products derived from this software without specific prior written 
  permission. \n

THIS SOFTWARE IS PROVIDED BY THE 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.

\b NOTES: \n
This code was developed using rigourous unit testing for every function 
and operation. Despite any rigorous development process, bugs are
inevitable. Please report bugs and suggested fixes to glenn @ zenautics.com.\n
*/

#include <stdio.h>  // for FILE*
#include <stdlib.h> // for calloc, malloc, free
#include <string.h> // for strlen, sprintf, strstr, strcmp, and others
#include <ctype.h>  // for isalpha
#include <math.h>
#include <float.h>
#include <errno.h>

#include "cmatrix.h"

#ifndef _MATRIX_NO_PLOTTING
#include "cplot.h"   // for CPLOT - plotting capabilities directly to an image file.
#endif

// deal with msvc empty projects
#ifndef WIN32
  #ifdef _WIN32
    #define WIN32
  #endif
#endif

#if defined _MSC_VER && _MSC_VER < 1400
#define _CRT_SECURE_NO_DEPRECATE
#endif

#ifndef _MSC_VER
#define _CRT_SECURE_NO_DEPRECATE
#endif


#include "kiss_fft.h" // Use kiss FFT, when INTEL_IPPS is disabled

//#define MTX_DEBUG
#ifdef MTX_DEBUG
#include <time.h>
#endif


#ifndef PI
#define PI (3.1415926535897932384626433832795) //!< better value
#endif

#ifndef TWOPI
#define TWOPI (6.283185307179586476925286766559) //!< 2.0*PI
#endif

#ifndef HALFPI
#define HALFPI (1.5707963267948966192313216916398) //!< PI/2.0
#endif


#define MTX_MAX_READFROMFILE_BUFFER (65536)
#define MATRIX_MIN_RLE_TOLERANCE (1.0e-16)
#define MTX_MAX_COMMENT_LENGTH (1024*6)

// binary version identifiers
#define MTX_ID_SIZE (8)
#define MTX_ID_COMPRESSED_01 ("MTX01\n") //!< identifier used to indicate a file stored using SaveCompressed (version 1)
#define MTX_ID_LEGACY_V01 ("Matrix\n") //!< legacy identifier used to indicate a file stored using Save with imprecise double RLE
#define MTX_ID_LEGACY_V02 ("MTXV02\n") //!< legacy identifier used to indicate a file stored using Save with imprecise double RLE

#define MTX_VERSION_NR_DEFAULT (101) //!< identifier used to indicate a file stored using basic Save
#define MTX_VERSION_NR_COMPRESSED_01 (102) //!< identifier used to indicate a file stored using SaveCompressed (version 1)
#define MTX_VERSION_NR_LEGACY_V01 (1) //!< legacy identifier used to indicate a file stored using Save with imprecise double RLE
#define MTX_VERSION_NR_LEGACY_V02 (2) //!< legacy identifier used to indicate a file stored using Save with imprecise double RLE

#define MTX_NK (8) //!< the number of byte columns used to represent a column of doubles when using MTX_ID_COMPRESSED_01

#define MTX_NAN     (sqrt(-1.0))
#define MTX_POS_INF (-log(0.0))
#define MTX_NEG_INF ( log(0.0))


/// \brief This static global variable indicates whether matrix
/// operations for single elements are treated as scalar
/// operations.
/// e.g. A = B*C, B is 1x1 C is 10x10. If this was disabled, this
/// operation would return FALSE as an error. When enabled, A, is
/// treated as a scalar and is multiplied into every element of C.
BOOL MTX_static_global_treat_1x1_as_scalar = TRUE;


typedef struct
{
  char id[8];
  unsigned headersize;
  unsigned isReal;
  unsigned nrows;
  unsigned ncols;
  unsigned filesize;
  unsigned crc;
  char *comment;
}_MTX_STRUCT_FileHeader;

typedef struct
{
  unsigned length[MTX_NK]; // length of compressed data column
  unsigned char isCompressed[MTX_NK]; // indicates which columns are compressed using RLE
  unsigned totalLength; // total length of MTX_NK data columns
}_MTX_STRUCT_CompressedColumnHeader;


/// struct specific for MTX_ReadFromFile and related functions (a simple linked list)
typedef struct _MTX_listItemCplx
{
  BOOL isReal; //!< A boolean to indicate if real data or complex data is attached to this item.
  double *rowptr; //!< The pointer to real data.
  stComplex *rowptr_cplx; //!< The pointer to complex data.
  struct _MTX_listItemCplx *next; //!< The pointer to the next item in the list.
}_MTX_STRUCT_ReadFromFileListElem;


/// static function for matrix memory allocation
static BOOL MTX_static_alloc( MTX *M, const unsigned nrows, const unsigned ncols, const BOOL setToZero, const BOOL isReal );

/// static function for converting a complex stored matrix to a real matrix (either all real component or all imaginary component)
static BOOL MTX_static_ConvertComplexTo( MTX *M, BOOL useReal );

static BOOL MTX_static_get_row_array_from_string( char *datastr, _MTX_STRUCT_ReadFromFileListElem *L, const unsigned ncols );

/// This function gets the next valid line of data. Whitespace lines are skipped.
static BOOL MTX_static_get_next_valid_data_line(
  FILE *in, //!< The input file pointer (input).
  char *linebuf, //!< A exisiting buffer to store the input line (input/output).
  unsigned *line_length, //!< The length of the line read (output).
  BOOL *atEOF //!< A boolean to indicate if EOF has been reached.
  );

/// This function gets the next valid line of data from a matrix string. Whitespace lines are skipped.
static BOOL MTX_static_get_next_valid_data_line_from_matrix_string(
  const char *strMatrix, //!< The matrix string pointer (input).
  const unsigned strLength, //!< The length of the matrix string (input).
  unsigned *index, //!< The starting/(next line) index into the matrix string pointer (input/output).
  char *linebuf, //!< A exisiting buffer to store the input line (input/output).
  unsigned *line_length, //!< The length of the line read (output).
  BOOL *atEndOfString //!< A boolean to indicate if the end of the strMatrix string has been reached.
  );


/// Extract a complex value from a string with a leading digit.
/// The string is either bi or a+bi.
static BOOL MTX_static_extract_cplx_from_string_with_leading_digit(
  char *datastr, //!< The entire input data string.
  const unsigned indexS, //!< The start index of the complex element.
  const unsigned indexE, //!< The inclusive end index of the complex element.
  double *re, //!< The extracted real component.
  double *im //!< The extracted imag component.
  );


static BOOL MTX_static_extract_real_into_cplx_from_string(
  char *datastr, //!< The entire input data string.
  const unsigned indexS, //!< The start index of the complex element.
  double *re, //!< The extracted real component.
  double *im //!< The extracted imag component.
  );

/// This static function looks for complex data in a line string.
static BOOL MTX_static_look_for_complex_data(
  char *linebuf, //!< A string containing a line of data (input).
  const unsigned line_length, //!< The length of the string (input).
  BOOL *hasComplex //!< A boolean indicating if there is any complex data (output).
  );



/// static function, rounds a value to an integer
static BOOL MTX_static_round_value_to_integer( double *value );

/// static function, rounds a value at the specified precision
static BOOL MTX_static_round_value( double *value, const unsigned precision );

////
// functions for quicksorting
static void MTX_static_quicksort( double *a, unsigned start, unsigned end ); //!< The normal quicksort function
static void MTX_static_swap_doubles( double *a, double *b ); //!< swap two doubles a and b
static int MTX_static_partition( double *a, unsigned start, unsigned end ); //!< partition the vector
static void MTX_static_quicksort_indexed( double *a, double *index, unsigned start, unsigned end ); //!< quicksort that also returns a sorted indexing vector
static void MTX_static_swap_doubles_indexed( double *a, double *b, double *index_a, double *index_b ); //!< swap the doubles and indexes
static int MTX_static_partition_indexed( double *a, double *index, unsigned start, unsigned end ); //!< partition the vectors


/// Compute the sqrt of a complex value.
static void MTX_static_quick_sqrt( const double *a_re, const double *a_im, double *re, double *im );


/// A static function to multiply a*b complex values
static void MTX_static_quick_complex_mult_ab(
  const double* a_re,
  const double* a_im,
  const double* b_re,
  const double* b_im,
  double *re,
  double *im );

/// A static function to multiply a*b*c complex values
static void MTX_static_quick_complex_mult_abc(
  const double* a_re,
  const double* a_im,
  const double* b_re,
  const double* b_im,
  const double* c_re,
  const double* c_im,
  double *re,
  double *im );


/// A static function to compute the complex result of a/b.
static void MTX_static_quick_complex_divide(
                                     const double* a_re, //!< The real part of a (input).
                                     const double* a_im, //!< The imag part of a (input).
                                     const double* b_re, //!< The real part of b (input).
                                     const double* b_im, //!< The imag part of b (input).
                                     double *re, //!< The real part of the result.
                                     double *im ); //!< The imag part of the result.



/// Calculate a CRC value to be used by CRC calculation functions.
static unsigned MTX_static_CRC32(unsigned ulCRC);

/// Updates the 32 bit CRC with a block of data.
/// This function can be called once (with uiCRC initialized to zero) to get the crc
/// for a single byte vector or multiple times to apply the crc calculation to multiple
/// bytes vectors.
static void MTX_static_updateCRC( unsigned char *pBytes, const unsigned nBytes, unsigned *uiCRC );

/// closes the file and frees memory used by MTX_SaveCompressed and MTX_Load
static void MTX_static_SaveAndLoadCleanUp( FILE *fid, unsigned char **bytes, unsigned char **compressed, const unsigned nk );

/// loads a legacy verison of a .mtx binary matrix file
static BOOL MTX_static_ReadCompressed_LegacyVersion( MTX* M, const char *path );


/// Performs factorization by gaussian elimination with scaled parital pivoting
/// /b Reference /n
/// [1] Chaney, Ward & David Kincaid, "Numerical Mathematics and Computing, 3rd Edition",
/// Cole Publishing Co., 1994, Belmont, CA, p.237)
static BOOL MTX_static_Factorize( BOOL *isFullRank, const unsigned n, unsigned* index, MTX *A );

/// Solve AX=b
/// factorized A is obtained from MTX_static_Factorize
static BOOL MTX_static_SolveByGaussianElimination(
  const MTX *b,
  MTX *X,
  const MTX *A, // factorized A
  unsigned *index );

/// Clean up dynamic memory used in MTX_Det.
static void MTX_static_Det_cleanup( unsigned *index, double *scale, MTX *U, MTX *magMtx );



/// Perform the FFT or IFFT of the columns in the src matrix and
/// store the result in the dst matrix. If the number of rows in the
/// src matrix is not a power of two, the DFT or IDFT is performed.
static BOOL MTX_static_fft(
                           const MTX *src, //!< The source matrix.
                           MTX *dst, //!< The result matrix (always complex).
                           BOOL isFwd //!< A boolean to indicate if this is a fwd transform or the inverse transform
                           );

/// Perform the FFT or IFFT of the columns in the src matrix inplace.
/// If the number of rows in the src matrix is not a power of two, the DFT or IDFT is performed.
static BOOL MTX_static_fft_inplace(
                                   MTX *src, //!< The source matrix.
                                   BOOL isFwd //!< A boolean to indicate if this is a fwd transform or the inverse transform
                                   );




/// \brief  Get a value from the uniform distribution [0,1].
/// \pre srand(seed) has been called.
static double MTX_static_get_rand_value();

/// \brief  Get a value from the standard normal gaussian distribution.
///
/// \pre srand(seed) has been called.
///
/// REFERENCE: \n
/// Scheinerman, E. R (2006). "C++ for Mathematicians: An Introduction for Students and Professionals."
/// Chapman and Hall/CRC, Taylor and Francis Group. pp 61-63.
static double MTX_static_get_randn_value();



BOOL MTX_Initialize_MTXEngine()
{
  return TRUE;
}

BOOL MTX_Enable1x1MatricesForTreatmentAsScalars( BOOL enable )
{
  MTX_static_global_treat_1x1_as_scalar = enable;
  return TRUE;
}

BOOL MTX_isNull( const MTX *M )
{
  if( !M )
    return TRUE;

  if( M->data == NULL && M->cplx == NULL )
    return TRUE;

  return FALSE;
}

BOOL MTX_isConformalForMultiplication( const MTX *A, const MTX *B )
{
  if( MTX_isNull( A ) )
  {
    return FALSE;
  }
  if( MTX_isNull( B ) )
  {
    return FALSE;
  }

  return( A->ncols == B->nrows );
}

BOOL MTX_isConformalForAddition( const MTX *A, const MTX *B )
{
  if( MTX_isNull( A ) )
  {
    return FALSE;
  }
  if( MTX_isNull( B ) )
  {
    return FALSE;
  }

  return( A->nrows == B->nrows && A->ncols == B->ncols );
}


BOOL MTX_isSquare( const MTX *A )
{
  if( MTX_isNull( A ) )
  {
    return FALSE;
  }

  return( A->nrows == A->ncols );
}

BOOL MTX_isSameSize( const MTX *A, const MTX *B )
{
  return MTX_isConformalForAddition( A, B );
}


BOOL MTX_Init( MTX *M )
{
  if( !M )
  {
    MTX_ERROR_MSG( "Cannot initialize NULL pointer." )
      return FALSE;
  }

  M->ncols = 0;
  M->nrows = 0;
  M->isReal = TRUE;
  M->cplx = NULL;
  M->data = NULL;
  M->comment = NULL;

  return TRUE;
}

BOOL MTX_SetComment( MTX *M, const char *comment )
{
  unsigned length;

  if( !M )
  {
    MTX_ERROR_MSG( "Cannot initialize NULL pointer." );
    return FALSE;
  }

  if( !comment )
  {
    MTX_ERROR_MSG( "if( !comment )" );
    return FALSE;
  }

  if( M->comment )
    free( M->comment );

  length = (unsigned int)strlen(comment);
  if( length == 0 )
  {
    MTX_ERROR_MSG( "strlen returned 0." );
    return FALSE;
  }

  M->comment = (char*)malloc( (length+1)*sizeof(char) ); // +1 for the null terminator
  if( !M->comment )
  {
    // memory allocation failure
    MTX_ERROR_MSG( "malloc returned NULL." );
    return FALSE;
  }

#ifndef _CRT_SECURE_NO_DEPRECATE
  if( strcpy_s( M->comment, length+1, comment ) != 0 )
  {
    MTX_ERROR_MSG( "strcpy_s returned 0." );
    free(M->comment);
    M->comment = NULL;
    return FALSE;
  }
#else
  strcpy( M->comment, comment );
#endif

  return TRUE;
}

BOOL MTX_Free( MTX *M )
{
  unsigned j = 0;

  if( !M )
  {
    MTX_ERROR_MSG( "Cannot free NULL pointer." );
    return FALSE;
  }

  if( M->isReal )
  {
    if( M->data == NULL )
    {
      if( M->comment )
        free( M->comment );

      M->comment = NULL;
      M->nrows = 0;
      M->ncols = 0;
      return TRUE;
    }
  }
  else
  {
    if( M->cplx == NULL )
    {
      if( M->comment )
        free( M->comment );

      M->comment = NULL;
      M->nrows = 0;
      M->ncols = 0;
      return TRUE;
    }
  }


  if( M->isReal )
  {
    for( j = 0; j < M->ncols; j++ )
    {
      free( M->data[j] );
    }
  }
  else
  {
    for( j = 0; j < M->ncols; j++ )
    {
      free( M->cplx[j] );
    }
  }

  // free the array of pointers
  if( M->isReal )
    free( M->data );
  else
    free( M->cplx );

  M->nrows = 0;
  M->ncols = 0;
  M->isReal = TRUE;
  M->cplx = NULL;
  M->data = NULL;

  if( M->comment )
    free( M->comment );
  M->comment = NULL;
  return TRUE;
}


BOOL MTX_Malloc( MTX *M, const unsigned nrows, const unsigned ncols, const BOOL isReal )
{
  return MTX_static_alloc( M, nrows, ncols, FALSE, isReal );
}

BOOL MTX_Calloc( MTX *M, const unsigned nrows, const unsigned ncols, const BOOL isReal )
{
  return MTX_static_alloc( M, nrows, ncols, TRUE, isReal );
}

BOOL MTX_static_alloc( MTX *M, const unsigned nrows, const unsigned ncols, const BOOL setToZero, const BOOL isReal )
{
  unsigned i = 0;
  unsigned j = 0;

  // invalid call
  if( nrows == 0 || ncols == 0 )
  {
    MTX_ERROR_MSG( "if( nrows == 0 || ncols == 0 )" );
    return FALSE;
  }
  if( !M )
  {
    MTX_ERROR_MSG( "Cannot set a NULL pointer." );
    return FALSE;
  }

  // Check if the matrix is already the right size and type.
  if( M->isReal == isReal )
  {
    if( M->nrows > 0 && M->ncols > 0 )
    {
      if( M->nrows == nrows && M->ncols == ncols )
      {
        // already the right size and type
        if( setToZero )
        {
          if( !MTX_Zero( M ) )
          {
            MTX_ERROR_MSG( "MTX_Zero returned FALSE." );
            return FALSE;