vtkPCAStatistics.cxx

#include "vtkPCAStatistics.h"

#include "vtkDoubleArray.h"
#include "vtkIdTypeArray.h"
#include "vtkInformation.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkMultiCorrelativeStatisticsAssessFunctor.h"
#include "vtkObjectFactory.h"
#include "vtkSmartPointer.h"
#include "vtkStringArray.h"
#include "vtkTable.h"
#include "vtkVariantArray.h"

#include <map>
#include <vector>
#include <sstream>

#include "alglib/svd.h"

// To Do:
// - Add option to pre-multiply EigenVectors by normalization coeffs
// - In vtkPCAAssessFunctor, pre-multiply EigenVectors by normalization coeffs (if req)
// -

#define VTK_PCA_NORMCOLUMN "PCA Cov Norm"
#define VTK_PCA_COMPCOLUMN "PCA"

vtkObjectFactoryNewMacro(vtkPCAStatistics)

const char* vtkPCAStatistics::NormalizationSchemeEnumNames[NUM_NORMALIZATION_SCHEMES + 1] =
{
  "None",
  "TriangleSpecified",
  "DiagonalSpecified",
  "DiagonalVariance",
  "InvalidNormalizationScheme"
};

const char* vtkPCAStatistics::BasisSchemeEnumNames[NUM_BASIS_SCHEMES + 1] =
{
  "FullBasis",
  "FixedBasisSize",
  "FixedBasisEnergy",
  "InvalidBasisScheme"
};

// ----------------------------------------------------------------------
void vtkPCAStatistics::GetEigenvalues(int request, vtkDoubleArray* eigenvalues)
{
  vtkSmartPointer<vtkMultiBlockDataSet> outputMetaDS =
      vtkMultiBlockDataSet::SafeDownCast(
        this->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );

  if(!outputMetaDS)
  {
    vtkErrorMacro(<<"NULL dataset pointer!");
  }

  vtkSmartPointer<vtkTable> outputMeta =
    vtkTable::SafeDownCast( outputMetaDS->GetBlock( request + 1 ) );

  if(!outputMetaDS)
  {
    vtkErrorMacro(<<"NULL table pointer!");
  }

  vtkDoubleArray* meanCol = vtkArrayDownCast<vtkDoubleArray>(outputMeta->GetColumnByName("Mean"));
  vtkStringArray* rowNames = vtkArrayDownCast<vtkStringArray>(outputMeta->GetColumnByName("Column"));

  eigenvalues->SetNumberOfComponents(1);

  // Get values
  int eval = 0;
  for(vtkIdType i = 0; i < meanCol->GetNumberOfTuples(); i++)
  {
    std::stringstream ss;
    ss << "PCA " << eval;

    std::string rowName = rowNames->GetValue(i);
    if(rowName.compare(ss.str()) == 0)
    {
      eigenvalues->InsertNextValue(meanCol->GetValue(i));
      eval++;
    }
  }

}

// ----------------------------------------------------------------------
double vtkPCAStatistics::GetEigenvalue(int request, int i)
{
  vtkSmartPointer<vtkDoubleArray> eigenvalues =
    vtkSmartPointer<vtkDoubleArray>::New();
  this->GetEigenvalues(request, eigenvalues);
  return eigenvalues->GetValue(i);
}

// ----------------------------------------------------------------------
void vtkPCAStatistics::GetEigenvalues(vtkDoubleArray* eigenvalues)
{
  this->GetEigenvalues(0, eigenvalues);
}

// ----------------------------------------------------------------------
double vtkPCAStatistics::GetEigenvalue(int i)
{
  return this->GetEigenvalue(0,i);
}

// ----------------------------------------------------------------------
void vtkPCAStatistics::GetEigenvectors(int request, vtkDoubleArray* eigenvectors)
{
  // Count eigenvalues
  vtkSmartPointer<vtkDoubleArray> eigenvalues =
    vtkSmartPointer<vtkDoubleArray>::New();
  this->GetEigenvalues(request, eigenvalues);
  vtkIdType numberOfEigenvalues = eigenvalues->GetNumberOfTuples();

  vtkSmartPointer<vtkMultiBlockDataSet> outputMetaDS =
      vtkMultiBlockDataSet::SafeDownCast(
        this->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );

  if(!outputMetaDS)
  {
    vtkErrorMacro(<<"NULL dataset pointer!");
  }

  vtkSmartPointer<vtkTable> outputMeta =
    vtkTable::SafeDownCast( outputMetaDS->GetBlock( request + 1 ) );

  if(!outputMeta)
  {
    vtkErrorMacro(<<"NULL table pointer!");
  }

  vtkDoubleArray* meanCol = vtkArrayDownCast<vtkDoubleArray>(outputMeta->GetColumnByName("Mean"));
  vtkStringArray* rowNames = vtkArrayDownCast<vtkStringArray>(outputMeta->GetColumnByName("Column"));

  eigenvectors->SetNumberOfComponents(numberOfEigenvalues);

  // Get vectors
  int eval = 0;
  for(vtkIdType i = 0; i < meanCol->GetNumberOfTuples(); i++)
  {
    std::stringstream ss;
    ss << "PCA " << eval;

    std::string rowName = rowNames->GetValue(i);
    if(rowName.compare(ss.str()) == 0)
    {
      std::vector<double> eigenvector;
      for(int val = 0; val < numberOfEigenvalues; val++)
      {
        // The first two columns will always be "Column" and "Mean", so start with the next one
        vtkDoubleArray* currentCol = vtkArrayDownCast<vtkDoubleArray>(outputMeta->GetColumn(val+2));
        eigenvector.push_back(currentCol->GetValue(i));
      }

      eigenvectors->InsertNextTypedTuple(&eigenvector.front());
      eval++;
    }
  }

}

// ----------------------------------------------------------------------
void vtkPCAStatistics::GetEigenvectors(vtkDoubleArray* eigenvectors)
{
  this->GetEigenvectors(0, eigenvectors);
}

// ----------------------------------------------------------------------
void vtkPCAStatistics::GetEigenvector(int request, int i, vtkDoubleArray* eigenvector)
{
  vtkSmartPointer<vtkDoubleArray> eigenvectors =
    vtkSmartPointer<vtkDoubleArray>::New();
  this->GetEigenvectors(request, eigenvectors);

  double* evec = new double[eigenvectors->GetNumberOfComponents()];
  eigenvectors->GetTypedTuple(i, evec);

  eigenvector->Reset();
  eigenvector->Squeeze();
  eigenvector->SetNumberOfComponents(eigenvectors->GetNumberOfComponents());
  eigenvector->InsertNextTypedTuple(evec);
  delete[] evec;
}

// ----------------------------------------------------------------------
void vtkPCAStatistics::GetEigenvector(int i, vtkDoubleArray* eigenvector)
{
  this->GetEigenvector(0, i, eigenvector);
}

// ======================================================== vtkPCAAssessFunctor
class vtkPCAAssessFunctor : public vtkMultiCorrelativeAssessFunctor
{
public:
  static vtkPCAAssessFunctor* New();

  vtkPCAAssessFunctor() { }
  ~vtkPCAAssessFunctor() VTK_OVERRIDE { }
  virtual bool InitializePCA(
                             vtkTable* inData, vtkTable* reqModel,
                             int normScheme, int basisScheme, int basisSize, double basisEnergy );

  void operator () ( vtkDoubleArray* result, vtkIdType row ) VTK_OVERRIDE;

  std::vector<double> EigenValues;
  std::vector<std::vector<double> > EigenVectors;
  vtkIdType BasisSize;
};

// ----------------------------------------------------------------------
vtkPCAAssessFunctor* vtkPCAAssessFunctor::New()
{
  return new vtkPCAAssessFunctor;
}

// ----------------------------------------------------------------------
bool vtkPCAAssessFunctor::InitializePCA( vtkTable* inData,
                                         vtkTable* reqModel,
                                         int normScheme,
                                         int basisScheme,
                                         int fixedBasisSize,
                                         double fixedBasisEnergy )
{
  if ( ! this->vtkMultiCorrelativeAssessFunctor::Initialize( inData,
                                                             reqModel,
                                                             false /* no Cholesky decomp */ ) )
  {
    return false;
  }

  // Put the PCA basis into a matrix form we can use.
  vtkIdType m = reqModel->GetNumberOfColumns() - 2;
  vtkDoubleArray* evalm = vtkArrayDownCast<vtkDoubleArray>( reqModel->GetColumnByName( VTK_MULTICORRELATIVE_AVERAGECOL ) );
  if ( ! evalm )
  {
    vtkGenericWarningMacro( "No \"" VTK_MULTICORRELATIVE_AVERAGECOL "\" column in request." );
    return false;
  }

  // Check that the derived model includes additional rows specifying the
  // normalization as required.
  vtkIdType mrmr = reqModel->GetNumberOfRows(); // actual number of rows
  vtkIdType ermr; // expected number of rows
  switch ( normScheme )
  {
    case vtkPCAStatistics::NONE:
      // m+1 covariance/Cholesky rows, m eigenvector rows, no normalization factors
      ermr = 2 * m + 1;
      break;
    case vtkPCAStatistics::DIAGONAL_SPECIFIED:
    case vtkPCAStatistics::DIAGONAL_VARIANCE:
      // m+1 covariance/Cholesky rows, m eigenvector rows, 1 normalization factor row
      ermr = 2 * m + 2;
      break;
    case vtkPCAStatistics::TRIANGLE_SPECIFIED:
      // m+1 covariance/Cholesky rows, m eigenvector rows, m normalization factor rows
      ermr = 3 * m + 1;
      break;
    case vtkPCAStatistics::NUM_NORMALIZATION_SCHEMES:
    default:
    {
      vtkGenericWarningMacro( "The normalization scheme specified (" << normScheme << ") is invalid." );
      return false;
    }
  }

  // Allow derived classes to add rows, but never allow fewer than required.
  if ( mrmr < ermr )
  {
    vtkGenericWarningMacro(
                           "Expected " << ( 2 * m + 1 ) << " or more columns in request but found "
                           << reqModel->GetNumberOfRows() << "." );
    return false;
  }

  // OK, we got this far; we should succeed.
  vtkIdType i, j;
  double eigSum = 0.;
  for ( i = 0; i < m; ++ i )
  {
    double eigVal = evalm->GetValue( m + 1 + i );
    eigSum += eigVal;
    this->EigenValues.push_back( eigVal );
  }
  this->BasisSize = -1;
  switch ( basisScheme )
  {
    case vtkPCAStatistics::NUM_BASIS_SCHEMES:
    default:
      vtkGenericWarningMacro( "Unknown basis scheme " << basisScheme << ". Using FULL_BASIS." );
      VTK_FALLTHROUGH;
    case vtkPCAStatistics::FULL_BASIS:
      this->BasisSize = m;
      break;
    case vtkPCAStatistics::FIXED_BASIS_SIZE:
      this->BasisSize = fixedBasisSize;
      break;
    case vtkPCAStatistics::FIXED_BASIS_ENERGY:
    {
      double frac = 0.;
      for ( i = 0; i < m; ++ i )
      {
        frac += this->EigenValues[i] / eigSum;
        if ( frac > fixedBasisEnergy )
        {
          this->BasisSize = i + 1;
          break;
        }
      }
      if ( this->BasisSize < 0 )
      { // OK, it takes all the eigenvectors to approximate that well...
        this->BasisSize = m;
      }
    }
      break;
  }

  // FIXME: Offer mode to include normalization factors (none,diag,triang)?
  // Could be done here by pre-multiplying this->EigenVectors by factors.
  for ( i = 0; i < this->BasisSize; ++ i )
  {
    std::vector<double> evec;
    for ( j = 0; j < m; ++ j )
    {
      evec.push_back( reqModel->GetValue( m + 1 + i, j + 2 ).ToDouble() );
    }
    this->EigenVectors.push_back( evec );
  }
  return true;
}

// ----------------------------------------------------------------------
void vtkPCAAssessFunctor::operator () ( vtkDoubleArray* result, vtkIdType row )
{
  vtkIdType i;
  result->SetNumberOfValues( this->BasisSize );
  std::vector<std::vector<double> >::iterator it;
  vtkIdType m = this->GetNumberOfColumns();
  for ( i = 0; i < m; ++ i )
  {
    this->Tuple[i] = this->Columns[i]->GetTuple( row )[0] - this->Center[i];
  }
  i = 0;
  for ( it = this->EigenVectors.begin(); it != this->EigenVectors.end(); ++ it, ++ i )
  {
    double cv = 0.;
    std::vector<double>::iterator tvit;
    std::vector<double>::iterator evit = this->Tuple.begin();
    for ( tvit = it->begin(); tvit != it->end(); ++ tvit, ++ evit )
    {
      cv += (*evit) * (*tvit);
    }
    result->SetValue( i, cv );
  }
}

// ======================================================== vtkPCAStatistics
vtkPCAStatistics::vtkPCAStatistics()
{
  this->SetNumberOfInputPorts( 4 ); // last port is for normalization coefficients.
  this->NormalizationScheme = NONE;
  this->BasisScheme = FULL_BASIS;
  this->FixedBasisSize = -1;
  this->FixedBasisEnergy = 1.;
}

// ----------------------------------------------------------------------
vtkPCAStatistics::~vtkPCAStatistics()
{
}

// ----------------------------------------------------------------------
void vtkPCAStatistics::PrintSelf( ostream& os, vtkIndent indent )
{
  this->Superclass::PrintSelf( os, indent );
  os << indent << "NormalizationScheme: " << this->GetNormalizationSchemeName( this->NormalizationScheme ) << "\n";
  os << indent << "BasisScheme: " << this->GetBasisSchemeName( this->BasisScheme ) << "\n";
  os << indent << "FixedBasisSize: " << this->FixedBasisSize << "\n";
  os << indent << "FixedBasisEnergy: " << this->FixedBasisEnergy << "\n";
}

// ----------------------------------------------------------------------
bool vtkPCAStatistics::SetParameter( const char* parameter,
                                     int vtkNotUsed( index ),
                                     vtkVariant value )
{
  if ( ! strcmp( parameter, "NormalizationScheme" ) )
  {
    this->SetNormalizationScheme( value.ToInt() );

    return true;
  }

  if ( ! strcmp( parameter, "BasisScheme" ) )
  {
    this->SetBasisScheme( value.ToInt() );

    return true;
  }

  if ( ! strcmp( parameter, "FixedBasisSize" ) )
  {
    this->SetFixedBasisSize( value.ToInt() );

    return true;
  }

  if ( ! strcmp( parameter, "FixedBasisEnergy" ) )
  {
    this->SetFixedBasisEnergy( value.ToDouble() );

    return true;
  }

  return false;
}

// ----------------------------------------------------------------------
const char* vtkPCAStatistics::GetNormalizationSchemeName( int schemeIndex )
{
  if ( schemeIndex < 0 || schemeIndex > NUM_NORMALIZATION_SCHEMES )
  {
    return vtkPCAStatistics::NormalizationSchemeEnumNames[NUM_NORMALIZATION_SCHEMES];
  }
  return vtkPCAStatistics::NormalizationSchemeEnumNames[schemeIndex];
}

// ----------------------------------------------------------------------
void vtkPCAStatistics::SetNormalizationSchemeByName( const char* schemeName )
{
  for ( int i = 0; i < NUM_NORMALIZATION_SCHEMES; ++ i )
  {
    if ( ! strcmp( vtkPCAStatistics::NormalizationSchemeEnumNames[i], schemeName ) )
    {
      this->SetNormalizationScheme( i );
      return;
    }
  }
  vtkErrorMacro( "Invalid normalization scheme name \"" << schemeName << "\" provided." );
}


// ----------------------------------------------------------------------
vtkTable* vtkPCAStatistics::GetSpecifiedNormalization()
{
  return vtkTable::SafeDownCast( this->GetInputDataObject( 3, 0 ) );
}

void vtkPCAStatistics::SetSpecifiedNormalization( vtkTable* normSpec )
{
  this->SetInputData( 3, normSpec );
}

// ----------------------------------------------------------------------
const char* vtkPCAStatistics::GetBasisSchemeName( int schemeIndex )
{
  if ( schemeIndex < 0 || schemeIndex > NUM_BASIS_SCHEMES )
  {
    return vtkPCAStatistics::BasisSchemeEnumNames[NUM_BASIS_SCHEMES];
  }
  return vtkPCAStatistics::BasisSchemeEnumNames[schemeIndex];
}

// ----------------------------------------------------------------------
void vtkPCAStatistics::SetBasisSchemeByName( const char* schemeName )
{
  for ( int i = 0; i < NUM_BASIS_SCHEMES; ++ i )
  {
    if ( ! strcmp( vtkPCAStatistics::BasisSchemeEnumNames[i], schemeName ) )
    {
      this->SetBasisScheme( i );
      return;
    }
  }
  vtkErrorMacro( "Invalid basis scheme name \"" << schemeName << "\" provided." );
}

// ----------------------------------------------------------------------
int vtkPCAStatistics::FillInputPortInformation( int port, vtkInformation* info )
{
  if ( port == 3 )
  {
    info->Set( vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkTable" );
    info->Set( vtkAlgorithm::INPUT_IS_OPTIONAL(), 1 );
    return 1;
  }
  return this->Superclass::FillInputPortInformation( port, info );
}

// ----------------------------------------------------------------------
static void vtkPCAStatisticsNormalizeSpec( vtkVariantArray* normData,
                                           ap::real_2d_array& cov,
                                           vtkTable* normSpec,
                                           vtkTable* reqModel,
                                           bool triangle )
{
  vtkIdType i, j;
  vtkIdType m = reqModel->GetNumberOfColumns() - 2;
  std::map<vtkStdString,vtkIdType> colNames;
  // Get a list of columns of interest for this request
  for ( i = 0; i < m; ++ i )
  {
    colNames[ reqModel->GetColumn( i + 2 )->GetName() ] = i;
  }
  // Turn normSpec into a useful array.
  std::map<std::pair<vtkIdType,vtkIdType>,double> factor;
  vtkIdType n = normSpec->GetNumberOfRows();
  for ( vtkIdType r = 0; r < n; ++ r )
  {
    std::map<vtkStdString,vtkIdType>::iterator it;
    if ( ( it = colNames.find( normSpec->GetValue( r, 0 ).ToString() ) ) == colNames.end() )
    {
      continue;
    }
    i = it->second;
    if ( ( it = colNames.find( normSpec->GetValue( r, 1 ).ToString() ) ) == colNames.end() )
    {
      continue;
    }
    j = it->second;
    if ( j < i )
    {
      vtkIdType tmp = i;
      i = j;
      j = tmp;
    }
    factor[std::pair<vtkIdType,vtkIdType>( i, j )] = normSpec->GetValue( r, 2 ).ToDouble();
  }
  // Now normalize cov, recording any missing factors along the way.
  std::ostringstream missing;
  bool gotMissing = false;
  std::map<std::pair<vtkIdType,vtkIdType>,double>::iterator fit;
  if ( triangle )
  { // Normalization factors are provided for the upper triangular portion of the covariance matrix.
    for ( i = 0; i < m; ++ i )
    {
      for ( j = i; j < m; ++ j )
      {
        double v;
        fit = factor.find( std::pair<vtkIdType,vtkIdType>( i, j ) );
        if ( fit == factor.end() )
        {
          v = 1.;
          gotMissing = true;
          missing
            << "(" << reqModel->GetColumn( i + 2 )->GetName()
            << "," << reqModel->GetColumn( j + 2 )->GetName()
            << ") ";
        }
        else
        {
          v = fit->second;
        }
        normData->InsertNextValue( v );
        cov( i, j ) /= v;
        if ( i != j )
        { // don't normalize diagonal entries twice
          cov( j, i ) /= v;
        }
      }
    }
  }
  else
  { // Only diagonal normalization factors are supplied. Off-diagonals are the product of diagonals.
    for ( i = 0; i < m; ++ i )
    {
      double v;
      double vsq;
      fit = factor.find( std::pair<vtkIdType,vtkIdType>( i, i ) );
      if ( fit == factor.end() )
      {
        vsq = v = 1.;
        gotMissing = true;
        missing
          << "(" << reqModel->GetColumn( i + 2 )->GetName()
          << "," << reqModel->GetColumn( i + 2 )->GetName()
          << ") ";
      }
      else
      {
        vsq = fit->second;
        v = sqrt( vsq );
      }
      normData->InsertNextValue( vsq );
      // normalization factor applied up and to the left.
      for ( j = 0; j < i; ++ j )
      {
        cov( i, j ) /= v;
        cov( j, i ) /= v;
      }
      // normalization factor applied down and to the right.
      for ( j = i + 1; j < m; ++ j )
      {
        cov( i, j ) /= v;
        cov( j, i ) /= v;
      }
      cov( i, i ) /= vsq;
    }
  }
  if ( gotMissing )
  {
    vtkGenericWarningMacro(
                           "The following normalization factors were expected but not provided: "
                           << missing.str().c_str() );
  }
}

// ----------------------------------------------------------------------
static void vtkPCAStatisticsNormalizeVariance( vtkVariantArray* normData,
                                               ap::real_2d_array& cov )
{
  vtkIdType i, j;
  vtkIdType m = cov.gethighbound( 0 ) - cov.getlowbound( 0 ) + 1;
  for ( i = 0; i < m; ++ i )
  {
    normData->InsertNextValue( cov( i, i ) );
    double norm = sqrt( cov( i, i ) );
    // normalization factor applied down and to the right.
    for ( j = i + 1; j < m; ++ j )
    {
      cov( i, j ) /= norm;
      cov( j, i ) /= norm;
    }
    // normalization factor applied up and to the left.
    for ( j = 0; j < i; ++ j )
    {
      cov( i, j ) /= norm;
      cov( j, i ) /= norm;
    }
    cov( i, i ) = 1.;
  }
}

// ----------------------------------------------------------------------
void vtkPCAStatistics::Derive( vtkMultiBlockDataSet* inMeta )
{
  if ( ! inMeta )
  {
    return;
  }

  // Use the parent class to compute a covariance matrix for each request.
  this->Superclass::Derive( inMeta );

  // Now that we have the covariance matrices, compute the SVD of each.
  vtkIdType nb = static_cast<vtkIdType>( inMeta->GetNumberOfBlocks() );
  for ( vtkIdType b = 1; b < nb; ++ b )
  {
    vtkTable* reqModel = vtkTable::SafeDownCast( inMeta->GetBlock( b ) );
    if ( ! reqModel )
    {
      continue;
    }
    vtkIdType m = reqModel->GetNumberOfColumns() - 2;
    ap::real_2d_array cov;
    cov.setbounds( 0, m - 1, 0, m - 1 );
    // Fill the cov array with values from the vtkTable
    vtkIdType i, j;
    for ( j = 2; j < 2 + m; ++ j )
    {
      for ( i = 0; i < j - 1; ++ i )
      {
        cov( i, j - 2 ) = reqModel->GetValue( i, j ).ToDouble();
      }
    }
    // Fill in the lower triangular portion of the matrix
    for ( j = 0; j < m - 1; ++ j )
    {
      for ( i = j; i < m; ++ i )
      {
        cov( i, j ) = cov( j, i );
      }
    }
    // If normalization of the covariance array is requested, perform it:
    vtkVariantArray* normData = vtkVariantArray::New();
    switch ( this->NormalizationScheme )
    {
      case TRIANGLE_SPECIFIED:
      case DIAGONAL_SPECIFIED:
        vtkPCAStatisticsNormalizeSpec(
                                      normData, cov, this->GetSpecifiedNormalization(), reqModel,
                                      this->NormalizationScheme == TRIANGLE_SPECIFIED );
        break;
      case DIAGONAL_VARIANCE:
        vtkPCAStatisticsNormalizeVariance( normData, cov );
        break;
      case NONE:
      case NUM_NORMALIZATION_SCHEMES:
      default:
        // do nothing
        break;
    }
    ap::real_2d_array u;
    ap::real_1d_array s;
    ap::real_2d_array vt;
    // Now that we have the covariance matrix, compute the SVD.
    // Note that vt is not computed since the VtNeeded parameter is 0.
    bool status = rmatrixsvd( cov, m, m, 2, 0, 2, s, u, vt );
    if ( ! status )
    {
      vtkWarningMacro( "Could not compute PCA for request " << b );
      continue;
    }
    vtkVariantArray* row = vtkVariantArray::New();
    //row->SetNumberOfComponents( m + 2 );
    //row->SetNumberOfTuples( 1 );
    row->SetNumberOfComponents( 1 );
    row->SetNumberOfTuples( m + 2 );
    for ( i = 0; i < m; ++ i )
    {
      std::ostringstream pcaCompName;
      pcaCompName << VTK_PCA_COMPCOLUMN << " " << i;
      row->SetValue( 0, pcaCompName.str().c_str() );
      row->SetValue( 1, s( i ) );
      for ( j = 0; j < m; ++ j )
      {
        // transpose the matrix so basis is row vectors (and thus
        // eigenvalues are to the left of their eigenvectors):
        row->SetValue( j + 2, u( j, i ) );
      }
      reqModel->InsertNextRow( row );
    }
    // Now insert the subset of the normalization data we used to
    // process this request at the bottom of the results.
    switch ( this->NormalizationScheme )
    {
      case TRIANGLE_SPECIFIED:
        for ( i = 0; i < m; ++ i )
        {
          std::ostringstream normCompName;
          normCompName << VTK_PCA_NORMCOLUMN << " " << i;
          row->SetValue( 0, normCompName.str().c_str() );
          row->SetValue( 1, 0. );
          for ( j = 0; j < i; ++ j )
          {
            row->SetValue( j + 2, 0. );
          }
          for ( ; j < m; ++ j )
          {
            row->SetValue( j + 2, normData->GetValue( j ) );
          }
          reqModel->InsertNextRow( row );
        }
        break;
      case DIAGONAL_SPECIFIED:
      case DIAGONAL_VARIANCE:
      {
        row->SetValue( 0, VTK_PCA_NORMCOLUMN );
        row->SetValue( 1, 0. );
        for ( j = 0; j < m; ++ j )
        {
          row->SetValue( j + 2, normData->GetValue( j ) );
        }
        reqModel->InsertNextRow( row );
      }
        break;
      case NONE:
      case NUM_NORMALIZATION_SCHEMES:
      default:
        // do nothing
        break;
    }
    normData->Delete();
    row->Delete();
  }
}
// Use the invalid value of -1 for p-values if R is absent
vtkDoubleArray* vtkPCAStatistics::CalculatePValues( vtkIdTypeArray* vtkNotUsed( dimCol ),
                                                    vtkDoubleArray* statCol )
{
  // A column must be created first
  vtkDoubleArray* testCol = vtkDoubleArray::New();

  // Fill this column
  vtkIdType n = statCol->GetNumberOfTuples();
  testCol->SetNumberOfTuples( n );
  for ( vtkIdType r = 0; r < n; ++ r )
  {
    testCol->SetTuple1( r, -1 );
  }

  return testCol;
}

// ----------------------------------------------------------------------
void vtkPCAStatistics::Test( vtkTable* inData,
                             vtkMultiBlockDataSet* inMeta,
                             vtkTable* outMeta )

{
  if ( ! inMeta )
  {
    return;
  }

  if ( ! outMeta )
  {
    return;
  }

  // Prepare columns for the test:
  // 0: (derived) model block index
  // 1: multivariate Srivastava skewness
  // 2: multivariate Srivastava kurtosis
  // 3: multivariate Jarque-Bera-Srivastava statistic
  // 4: multivariate Jarque-Bera-Srivastava p-value (calculated only if R is available, filled with -1 otherwise)
  // 5: number of degrees of freedom of Chi square distribution
  // NB: These are not added to the output table yet, for they will be filled individually first
  //     in order that R be invoked only once.
  vtkIdTypeArray* blockCol = vtkIdTypeArray::New();
  blockCol->SetName( "Block" );

  vtkDoubleArray* bS1Col = vtkDoubleArray::New();
  bS1Col->SetName( "Srivastava Skewness" );

  vtkDoubleArray* bS2Col = vtkDoubleArray::New();
  bS2Col->SetName( "Srivastava Kurtosis" );

  vtkDoubleArray* statCol = vtkDoubleArray::New();
  statCol->SetName( "Jarque-Bera-Srivastava" );

  vtkIdTypeArray* dimCol = vtkIdTypeArray::New();
  dimCol->SetName( "d" );

  // Retain data cardinality to check that models are applicable
  vtkIdType nRowData = inData->GetNumberOfRows();

  // Now iterate over model blocks
  unsigned int nBlocks = inMeta->GetNumberOfBlocks();
  for ( unsigned int b = 1; b < nBlocks; ++ b )
  {
    vtkTable* derivedTab = vtkTable::SafeDownCast( inMeta->GetBlock( b ) );

    // Silenty ignore empty blocks
    if ( ! derivedTab )
    {
      continue;
    }

    // Figure out dimensionality; it is assumed that the 2 first columns
    // are what they should be: namely, Column and Mean.
    int p = derivedTab->GetNumberOfColumns() - 2;

    // Return informative message when cardinalities do not match.
    if ( derivedTab->GetValueByName( p, "Mean" ).ToInt() != nRowData )
    {
      vtkWarningMacro( "Inconsistent input: input data has "
                       << nRowData
                       << " rows but primary model has cardinality "
                       << derivedTab->GetValueByName( p, "Mean" ).ToInt()
                       << " for block "
                       << b
                       <<". Cannot test." );
      continue;
    }

    // Create and fill entries of name and mean vectors
    vtkStdString *varNameX = new vtkStdString[p];
    double *mX = new double[p];
    for ( int i = 0; i < p; ++ i )
    {
      varNameX[i] = derivedTab->GetValueByName( i, "Column" ).ToString();
      mX[i] = derivedTab->GetValueByName( i, "Mean" ).ToDouble();
    }

    // Create and fill entries of eigenvalue vector and change of basis matrix
    double *wX = new double[p];
    double *P = new double[p * p];
    for ( int i = 0; i < p; ++ i )
    {
      // Skip p + 1 (Means and Cholesky) rows and 1 column (Column)
      wX[i] = derivedTab->GetValue( i + p + 1, 1).ToDouble();

      for ( int j = 0; j < p; ++ j )
      {
        // Skip p + 1 (Means and Cholesky) rows and 2 columns (Column and Mean)
        P[p * i + j] = derivedTab->GetValue( i + p + 1, j + 2).ToDouble();
      }
    }

    // Now iterate over all observations
    double tmp, t;
    double *x = new double[p];
    double *sum3 = new double[p];
    double *sum4 = new double[p];
    for ( int i = 0; i < p; ++ i )
    {
      sum3[i] = 0.;
      sum4[i] = 0.;
    }
    for ( vtkIdType r = 0; r < nRowData; ++ r )
    {
      // Read and center observation
      for ( int i = 0; i < p; ++ i )
      {
        x[i] = inData->GetValueByName( r, varNameX[i] ).ToDouble() - mX[i];
      }

      // Now calculate skewness and kurtosis per component
      for ( int i = 0; i < p; ++ i )
      {
        // Transform coordinate into eigencoordinates
        t = 0.;
        for ( int j = 0; j < p; ++ j )
        {
          // Pij = P[p*i+j]
          t += P[p * i + j] * x[j];
        }

        // Update third and fourth order sums for each eigencoordinate
        tmp = t * t;
        sum3[i] += tmp * t;
        sum4[i] += tmp * tmp;
      }
    }

    // Finally calculate moments by normalizing sums with corresponding eigenvalues and powers
    double bS1 = 0.;
    double bS2 = 0.;
    for ( int i = 0; i < p; ++ i )
    {
      tmp = wX[i] * wX[i];
      bS1 += sum3[i] * sum3[i] / ( tmp * wX[i] );
      bS2 += sum4[i] / tmp;
    }
    bS1 /= ( nRowData * nRowData * p );
    bS2 /= ( nRowData * p );

    // Finally, calculate Jarque-Bera-Srivastava statistic
    tmp = bS2 - 3.;
    double jbs = static_cast<double>( nRowData * p ) * ( bS1 / 6. + ( tmp * tmp ) / 24. );

    // Insert variable name and calculated Jarque-Bera-Srivastava statistic
    blockCol->InsertNextValue( b );
    bS1Col->InsertNextTuple1( bS1 );
    bS2Col->InsertNextTuple1( bS2 );
    statCol->InsertNextTuple1( jbs );
    dimCol->InsertNextTuple1( p + 1 );

    // Clean up
    delete [] sum3;
    delete [] sum4;
    delete [] x;
    delete [] P;
    delete [] wX;
    delete [] mX;
    delete [] varNameX;
  } // b

  // Now, add the already prepared columns to the output table
  outMeta->AddColumn( blockCol );
  outMeta->AddColumn( bS1Col );
  outMeta->AddColumn( bS2Col );
  outMeta->AddColumn( statCol );
  outMeta->AddColumn( dimCol );

  // Last phase: compute the p-values or assign invalid value if they cannot be computed
  vtkDoubleArray* testCol = this->CalculatePValues( dimCol, statCol );

  // The test column name can only be set after the column has been obtained from R
  testCol->SetName( "P" );
  // Now add the column of invalid values to the output table
  outMeta->AddColumn( testCol );
  // Clean up
  testCol->Delete();

  // Clean up
  blockCol->Delete();
  bS1Col->Delete();
  bS2Col->Delete();
  statCol->Delete();
  dimCol->Delete();
}
// ----------------------------------------------------------------------
void vtkPCAStatistics::Assess( vtkTable* inData,
                               vtkMultiBlockDataSet* inMeta,
                               vtkTable* outData )
{
  if ( ! inData )
  {
    return;
  }

  if ( ! inMeta )
  {
    return;
  }

  // For each request, add a column to the output data related to the likelihood of each input datum wrt the model in the request.
  // Column names of the metadata and input data are assumed to match.
  // The output columns will be named "RelDevSq(A,B,C)" where "A", "B", and "C" are the column names specified in the
  // per-request metadata tables.
  vtkIdType nRow = inData->GetNumberOfRows();
  int nb = static_cast<int>( inMeta->GetNumberOfBlocks() );
  AssessFunctor* dfunc = 0;
  for ( int req = 1; req < nb; ++ req )
  {
    vtkTable* reqModel = vtkTable::SafeDownCast( inMeta->GetBlock( req ) );
    if ( ! reqModel )
    { // silently skip invalid entries. Note we leave assessValues column in output data even when it's empty.
      continue;
    }

    this->SelectAssessFunctor( inData,
                               reqModel,
                               0,
                               dfunc );

    vtkPCAAssessFunctor* pcafunc = static_cast<vtkPCAAssessFunctor*>( dfunc );
    if ( ! pcafunc )
    {
      vtkWarningMacro( "Request " << req - 1 << " could not be accommodated. Skipping." );
      delete dfunc;
      continue;
    }

    // Create an array to hold the assess values for all the input data
    std::vector<double*> assessValues;
    int comp;
    for ( comp = 0; comp < pcafunc->BasisSize; ++ comp )
    {
      std::ostringstream reqNameStr;
      reqNameStr << VTK_PCA_COMPCOLUMN << "{";
      for ( int i = 0; i < pcafunc->GetNumberOfColumns(); ++ i )
      {
        if ( i > 0 )
        {
          reqNameStr << ",";
        }
        reqNameStr << pcafunc->GetColumn(i)->GetName();
      }
      reqNameStr << "}(" << comp << ")";
      vtkDoubleArray* arr = vtkDoubleArray::New();
      arr->SetName( reqNameStr.str().c_str() );
      arr->SetNumberOfTuples( nRow );
      outData->AddColumn( arr );
      arr->Delete();
      assessValues.push_back( arr->GetPointer( 0 ) );
    }

    // Something to hold assessed values for a single input datum
    vtkDoubleArray* singleResult = vtkDoubleArray::New();
    // Loop over all the input data and assess each datum:
    for ( vtkIdType row = 0; row < nRow; ++ row )
    {
      (*dfunc)( singleResult, row );
      for ( comp = 0; comp < pcafunc->BasisSize; ++ comp )
      {
        assessValues[comp][row] = singleResult->GetValue( comp );
      }
    }
    delete dfunc;
    singleResult->Delete();
  }
}

// ----------------------------------------------------------------------
void vtkPCAStatistics::SelectAssessFunctor( vtkTable* inData,
                                            vtkDataObject* inMetaDO,
                                            vtkStringArray* vtkNotUsed(rowNames),
                                            AssessFunctor*& dfunc )
{
  dfunc = 0;
  vtkTable* reqModel = vtkTable::SafeDownCast( inMetaDO );
  if ( ! reqModel )
  {
    return;
  }

  vtkPCAAssessFunctor* pcafunc = vtkPCAAssessFunctor::New();
  if ( ! pcafunc->InitializePCA(
                                inData, reqModel, this->NormalizationScheme,
                                this->BasisScheme, this->FixedBasisSize, this->FixedBasisEnergy ) )
  {
    delete pcafunc;
    return;
  }

  dfunc = pcafunc;
}