CorrelationMatrix.cpp

#include "CorrelationMatrix.h"

#include <QDataStream>
#include <QDebug>
#include <QFile>
#include <QList>
#include <QString>
#include <QStringList>
#include <QtCore/qmath.h>

#include "IException.h"
#include "LinearAlgebra.h"
#include "Pvl.h"
#include "PvlObject.h"
#include "SparseBlockMatrix.h"

namespace Isis {
  CorrelationMatrix::CorrelationMatrix() {
    m_covarianceFileName = new FileName(""); 
    m_correlationFileName = new FileName(""); 
    m_visibleBlocks = new QList<SparseBlockColumnMatrix>();
    m_imagesAndParameters = new QMap<QString, QStringList>();
    m_diagonals = new QList<double>();
  }


  CorrelationMatrix::CorrelationMatrix(PvlObject storedMatrixData) {
    //m_imagesAndParameters = NULL;
    m_imagesAndParameters = new QMap<QString, QStringList>();
    m_covarianceFileName = new FileName("");
    m_correlationFileName = new FileName("");
    m_diagonals = NULL;
    m_visibleBlocks = NULL;

    if (storedMatrixData.name() != "CorrelationMatrixData") {
      QString msg = "This Pvl Object does not have the correct correlation information. The Object "
                    "you are looking for is called CorrelationMatrixData.";
      throw IException(IException::User, msg, _FILEINFO_);
    }

    try {
      m_covarianceFileName =
          new FileName(storedMatrixData.findKeyword("CovarianceMatrixFileName")[0]);
    }
    catch (IException &e) {
      QString msg = "Could not find the Covariance Matrix .dat file name.";
      throw IException(e, IException::User, msg, _FILEINFO_);
    }
    
    try {
      QString corrFileName = storedMatrixData.findKeyword("CorrelationMatrixFileName")[0];
      if (corrFileName == "NULL") {
        m_correlationFileName = new FileName;
      }
      else {
        m_correlationFileName = new FileName(corrFileName);
      }
    }
    catch (IException &e) {
      QString msg = "Could not find the Correlation Matrix .dat file name.";
      throw IException(e, IException::User, msg, _FILEINFO_);
    }

    try {
      PvlObject::PvlKeywordIterator
          imgsIt = storedMatrixData.findGroup("ImagesAndParameters").begin();
      while ( imgsIt != storedMatrixData.findGroup("ImagesAndParameters").end() ) {
        QStringList params = (*imgsIt)[0].split(",");
        m_imagesAndParameters->insert(imgsIt->name(), params);
        imgsIt++;
      }
    }
    catch (IException &e) {
      QString msg = "Could not get Images and Parameters from ImagesAndParameters group.";
      throw IException(e, IException::User, msg, _FILEINFO_);
    }
  }


  CorrelationMatrix::CorrelationMatrix(const CorrelationMatrix &other) {
    m_imagesAndParameters = new QMap<QString, QStringList>(*other.m_imagesAndParameters);
    m_covarianceFileName = new FileName(*other.m_covarianceFileName);
    m_correlationFileName = new FileName(*other.m_correlationFileName);
    m_diagonals = new QList<double>(*other.m_diagonals);
    m_visibleBlocks = new QList<SparseBlockColumnMatrix>(*other.m_visibleBlocks);
  }


  CorrelationMatrix::~CorrelationMatrix() {
    delete m_imagesAndParameters;
    m_imagesAndParameters = NULL;

    delete m_covarianceFileName;
    m_covarianceFileName = NULL;

    delete m_correlationFileName;
    m_correlationFileName = NULL;

    delete m_diagonals;
    m_diagonals = NULL;

    delete m_visibleBlocks;
    m_visibleBlocks = NULL;
  }


  CorrelationMatrix &CorrelationMatrix::operator=(const CorrelationMatrix &other) {

    if (&other != this) {

      delete m_imagesAndParameters;
      m_imagesAndParameters = NULL;
      m_imagesAndParameters = new QMap<QString, QStringList>(*other.m_imagesAndParameters);
  
      delete m_covarianceFileName;
      m_covarianceFileName = NULL;
      m_covarianceFileName = new FileName(*other.m_covarianceFileName);
  
      delete m_correlationFileName;
      m_correlationFileName = NULL;
      m_correlationFileName = new FileName(*other.m_correlationFileName);
  
      delete m_diagonals;
      m_diagonals = NULL;
      m_diagonals = new QList<double>(*other.m_diagonals);
  
      delete m_visibleBlocks;
      m_visibleBlocks = NULL;
      m_visibleBlocks = new QList<SparseBlockColumnMatrix>(*other.m_visibleBlocks);

    }

    return *this;
  }


  void CorrelationMatrix::computeCorrelationMatrix() {

    if ( !isValid() ) {
      QString msg = "Cannot compute correlation matrix without a specified file name. Use "
                    "setCorrelationFileName(FileName) before calling computeCorrelationMatrix().";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    delete m_visibleBlocks;
    m_visibleBlocks = NULL;
    m_visibleBlocks = new QList<SparseBlockColumnMatrix>;

    // Create file handle
    QFile matrixInput( m_covarianceFileName->expanded() );
    QFile matrixOutput( m_correlationFileName->expanded() );

    // Open file to write to
    matrixInput.open(QIODevice::ReadOnly);
    matrixOutput.open(QIODevice::WriteOnly);

    // Open Stream
    QDataStream inStream(&matrixInput);
    QDataStream outStream(&matrixOutput);

    double firstParam1 = 0, //starting param for each iteration
            firstParam2 = 0;
    double param1 = 0, // current param for each iteration
            param2 = 0;
    // Read one column at a time
    SparseBlockColumnMatrix sbcm;
    while ( !inStream.atEnd() ) {
      inStream >> sbcm;

      // Store diagonal
      int numOfBlocks = sbcm.size();
      int lastBlock = numOfBlocks - 1;
      int numDiagonals = sbcm[lastBlock]->size1();

      // Get Diagonals
      for (int i = 0; i < numDiagonals; i++) {
        double val = ( *(sbcm[lastBlock]) )(i, i);
        m_diagonals->append(val);
      }

      // compute correlations
      QMapIterator<int, LinearAlgebra::Matrix *> block(sbcm);

      while ( block.hasNext() ) { // each block in the column
        block.next();
        for (int row = 0; row < (int)block.value()->size1(); row++) { // each row in the block
          for (int column = 0; column < (int)block.value()->size2(); column++) { // each column 
            // correlation = covariance / (variance1 * variance2)
            ( *block.value() )(row, column) = ( *block.value() )(row, column) /
                                              sqrt( (*m_diagonals)[param1] *
                                                    (*m_diagonals)[param2] );
            param2++;
          }
          param1++;
          param2 = firstParam2;
        }
        firstParam1 += block.value()->size1();
        param1 = firstParam1;
      }
      firstParam1 = 0;  // start each column at first element of diagonal list
      param1 = firstParam1;
      firstParam2 += block.value()->size2();
      param2 = firstParam2;

      outStream << sbcm;
      m_visibleBlocks->append(sbcm);
    }
    
    // close file
    matrixInput.close();
    matrixOutput.close();
  }



  void CorrelationMatrix::retrieveVisibleElements(int x, int y) {

//     if ( !correlationMatrixExists() ) {
      // call computeCorrelationMatrix
//     }
    // read the values we want from the correlation matrix file.

    // store values by column?
    // return list of values in m_visibleElements
  }



//   QList<MatrixElement*> CorrelationMatrix::visibleElements() {
//     return *m_visibleElements;
//   }



  bool CorrelationMatrix::isValid() {

      return !(m_correlationFileName->name() == "" || m_covarianceFileName->name() == "");
  }


  bool CorrelationMatrix::hasCovMat() {
    return !(m_covarianceFileName->name() == "");
  }


  // Set Methods
  void CorrelationMatrix::setCorrelationFileName(FileName correlationFileName) {
    if (m_correlationFileName == NULL) {
      m_correlationFileName = new FileName(correlationFileName);
    }
    else {
      *m_correlationFileName = correlationFileName;
    }
  }


  void CorrelationMatrix::setCovarianceFileName(FileName covarianceFileName) {
    if (m_covarianceFileName == NULL) {
      m_covarianceFileName = new FileName(covarianceFileName);
    }
    else {
      *m_covarianceFileName = covarianceFileName;
    }
    //Make the correlation matrix file name match the covariance matrix file name.
    if (!isValid()) {
      QString fName = covarianceFileName.expanded().replace( QString("inverse"),
                                                            QString("correlation") );
      setCorrelationFileName( FileName(fName) );
    }
  }


  void CorrelationMatrix::setImagesAndParameters(QMap<QString, QStringList> imagesAndParameters) {
    if (m_imagesAndParameters == NULL) {
      m_imagesAndParameters = new QMap<QString, QStringList>(imagesAndParameters);
    }
    else {
      *m_imagesAndParameters = imagesAndParameters;
    }      
  }


  FileName CorrelationMatrix::correlationFileName() {
    return *m_correlationFileName;
  }


  FileName CorrelationMatrix::covarianceFileName() {
    return *m_covarianceFileName;
  }


  QMap<QString, QStringList> *CorrelationMatrix::imagesAndParameters() {
    return m_imagesAndParameters;
  }


  void CorrelationMatrix::retrieveWholeMatrix() {
//     SparseBlockColumnMatrix sbcm;
//     QFile matrixInput( m_correlationFileName->expanded() );
//     matrixInput.open(QIODevice::ReadOnly);
//     QDataStream inStream(&matrixInput);
//
//     while( !inStream.atEnd() ) {
//       inStream >> sbcm;
//       m_visibleBlocks->append(&sbcm);
//     }
  }


  void CorrelationMatrix::retrieveThreeVisibleBlocks() {
  }


  QList<SparseBlockColumnMatrix> *CorrelationMatrix::visibleBlocks() {
    return m_visibleBlocks;
  }


  PvlObject CorrelationMatrix::pvlObject() {
    PvlObject corrMatInfo("CorrelationMatrixData");
    
    corrMatInfo += PvlKeyword( "CovarianceMatrixFileName", m_covarianceFileName->expanded() );
    corrMatInfo += PvlKeyword( "CorrelationMatrixFileName", m_correlationFileName->expanded() );

    PvlGroup imgsAndParams("ImagesAndParameters");
    QMapIterator<QString, QStringList> imgParamIt(*m_imagesAndParameters);
    while ( imgParamIt.hasNext() ) {
      imgParamIt.next();
      imgsAndParams += PvlKeyword( imgParamIt.key(), imgParamIt.value().join(",") );
    }
    corrMatInfo += imgsAndParams;
    
    return corrMatInfo;
  }
}