Commit 0cdf6f37 authored by George Zagaris's avatar George Zagaris
Browse files

Merge remote-tracking branch 'origin/master' into AMR-Refactoring

parents 5ac7dc34 0661c8ba
/*=========================================================================
Program: Visualization Toolkit
Module: TestRandomMomentStatisticsMPI.cxx
Program: Visualization Toolkit
Module: TestRandomMomentStatisticsMPI.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
/*
......@@ -48,7 +48,11 @@
struct RandomSampleStatisticsArgs
{
int nVals;
bool descOnly;
bool skipDescriptive;
bool skipPDescriptive;
bool skipPCorrelative;
bool skipPMultiCorrelative;
bool skipPPCA;
int* retVal;
int ioRank;
};
......@@ -69,7 +73,7 @@ void RandomSampleStatistics( vtkMultiProcessController* controller, void* arg )
// Seed random number generator
vtkMath::RandomSeed( static_cast<int>( vtkTimerLog::GetUniversalTime() ) * ( myRank + 1 ) );
// Generate an input table that contains samples of mutually independent random variables over [0, 1]
// Generate an input table that contains samples of mutually independent random variables
int nUniform = 2;
int nNormal = 2;
int nVariables = nUniform + nNormal;
......@@ -117,563 +121,603 @@ void RandomSampleStatistics( vtkMultiProcessController* controller, void* arg )
doubleArray[c]->Delete();
}
// "68-95-99.7 rule"
// Actually testing for 1, ..., numRuleVa standard deviations
int numRuleVal = 6;
// Reference values
double sigmaRuleVal[] = { 68.2689492137,
95.4499736104,
99.7300203937,
99.9936657516,
99.9999426697,
99.9999998027 };
// Tolerances
double sigmaRuleTol[] = { 1.,
.5,
.1,
.05,
.01,
.005 };
// ************************** Descriptive Statistics **************************
// Synchronize and start clock
com->Barrier();
// Create timer to be used by all tests
vtkTimerLog *timer=vtkTimerLog::New();
timer->StartTimer();
// For verification, instantiate a serial descriptive statistics engine and set its ports
vtkDescriptiveStatistics* ds = vtkDescriptiveStatistics::New();
ds->SetInput( vtkStatisticsAlgorithm::INPUT_DATA, inputData );
// ************************** Serial descriptive Statistics **************************
// Select all columns
for ( int c = 0; c < nVariables; ++ c )
// Skip serial descriptive statistics if requested
if ( ! args->skipDescriptive )
{
ds->AddColumn( columnNames[c] );
}
// Synchronize and start clock
com->Barrier();
timer->StartTimer();
// Test (serially) with Learn and Derive options only
ds->SetLearnOption( true );
ds->SetDeriveOption( true );
ds->SetAssessOption( false );
ds->SetTestOption( false );
ds->Update();
// For verification, instantiate a serial descriptive statistics engine and set its ports
vtkDescriptiveStatistics* ds = vtkDescriptiveStatistics::New();
ds->SetInput( vtkStatisticsAlgorithm::INPUT_DATA, inputData );
// Get output data and meta tables
vtkMultiBlockDataSet* outputMetaDS = vtkMultiBlockDataSet::SafeDownCast( ds->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );
vtkTable* outputPrimary = vtkTable::SafeDownCast( outputMetaDS->GetBlock( 0 ) );
vtkTable* outputDerived = vtkTable::SafeDownCast( outputMetaDS->GetBlock( 1 ) );
// Select all columns
for ( int c = 0; c < nVariables; ++ c )
{
ds->AddColumn( columnNames[c] );
}
// Collect (local) cardinalities, extrema, and means
int nRows = outputPrimary->GetNumberOfRows();
int np = com->GetNumberOfProcesses();
int n2Rows = 2 * nRows;
// Test (serially) with Learn operation only (this is only to verify parallel statistics)
ds->SetLearnOption( true );
ds->SetDeriveOption( false );
ds->SetAssessOption( false );
ds->SetTestOption( false );
ds->Update();
double* extrema_l = new double[n2Rows];
double* extrema_g = new double[n2Rows];
// Get output data and meta tables
vtkMultiBlockDataSet* outputMetaDS = vtkMultiBlockDataSet::SafeDownCast( ds->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );
vtkTable* outputPrimary = vtkTable::SafeDownCast( outputMetaDS->GetBlock( 0 ) );
double* cardsAndMeans_l = new double[n2Rows];
double* cardsAndMeans_g = new double[n2Rows];
// Collect (local) cardinalities, extrema, and means
int nRows = outputPrimary->GetNumberOfRows();
int np = com->GetNumberOfProcesses();
int n2Rows = 2 * nRows;
double dn;
for ( vtkIdType r = 0; r < nRows; ++ r )
{
dn = outputPrimary->GetValueByName( r, "Cardinality" ).ToDouble();
cardsAndMeans_l[2 * r] = dn;
cardsAndMeans_l[2 * r + 1] = dn * outputPrimary->GetValueByName( r, "Mean" ).ToDouble();
double* extrema_l = new double[n2Rows];
double* extrema_g = new double[n2Rows];
extrema_l[2 * r] = outputPrimary->GetValueByName( r, "Minimum" ).ToDouble();
// Collect -max instead of max so a single reduce op. (minimum) can process both extrema at a time
extrema_l[2 * r + 1] = - outputPrimary->GetValueByName( r, "Maximum" ).ToDouble();
}
double* cardsAndMeans_l = new double[n2Rows];
double* cardsAndMeans_g = new double[n2Rows];
// Reduce all extremal values, and gather all cardinalities and means, on process calcProc
int calcProc = np - 1;
double dn;
for ( vtkIdType r = 0; r < nRows; ++ r )
{
dn = outputPrimary->GetValueByName( r, "Cardinality" ).ToDouble();
cardsAndMeans_l[2 * r] = dn;
cardsAndMeans_l[2 * r + 1] = dn * outputPrimary->GetValueByName( r, "Mean" ).ToDouble();
com->Reduce( extrema_l,
extrema_g,
n2Rows,
vtkCommunicator::MIN_OP,
calcProc );
extrema_l[2 * r] = outputPrimary->GetValueByName( r, "Minimum" ).ToDouble();
// Collect -max instead of max so a single reduce op. (minimum) can process both extrema at a time
extrema_l[2 * r + 1] = - outputPrimary->GetValueByName( r, "Maximum" ).ToDouble();
}
com->Reduce( cardsAndMeans_l,
cardsAndMeans_g,
n2Rows,
vtkCommunicator::SUM_OP,
calcProc );
// Reduce all extremal values, and gather all cardinalities and means, on process calcProc
int calcProc = np - 1;
// Have process calcProc calculate global cardinality and mean, and send all results to I/O process
if ( myRank == calcProc )
{
if ( ! com->Send( extrema_g,
n2Rows,
args->ioRank,
65 ) )
{
vtkGenericWarningMacro("MPI error: process "<<myRank<< "could not send global results. Serial/parallel sanity check will be meaningless.");
*(args->retVal) = 1;
}
com->Reduce( extrema_l,
extrema_g,
n2Rows,
vtkCommunicator::MIN_OP,
calcProc );
if ( ! com->Send( cardsAndMeans_g,
n2Rows,
args->ioRank,
66 ) )
{
vtkGenericWarningMacro("MPI error: process "<<myRank<< "could not send global results. Serial/parallel sanity check will be meaningless.");
*(args->retVal) = 1;
}
}
com->Reduce( cardsAndMeans_l,
cardsAndMeans_g,
n2Rows,
vtkCommunicator::SUM_OP,
calcProc );
// Have I/O process receive results from process calcProc
if ( myRank == args->ioRank )
{
if ( ! com->Receive( extrema_g,
n2Rows,
calcProc,
65 ) )
// Have process calcProc calculate global cardinality and mean, and send all results to I/O process
if ( myRank == calcProc )
{
vtkGenericWarningMacro("MPI error: I/O process "<<args->ioRank<<" could not receive global results. Serial/parallel sanity check will be meaningless.");
*(args->retVal) = 1;
if ( ! com->Send( extrema_g,
n2Rows,
args->ioRank,
65 ) )
{
vtkGenericWarningMacro("MPI error: process "<<myRank<< "could not send global results. Serial/parallel sanity check will be meaningless.");
*(args->retVal) = 1;
}
if ( ! com->Send( cardsAndMeans_g,
n2Rows,
args->ioRank,
66 ) )
{
vtkGenericWarningMacro("MPI error: process "<<myRank<< "could not send global results. Serial/parallel sanity check will be meaningless.");
*(args->retVal) = 1;
}
}
if ( ! com->Receive( cardsAndMeans_g,
n2Rows,
calcProc,
66 ) )
// Have I/O process receive results from process calcProc
if ( myRank == args->ioRank )
{
vtkGenericWarningMacro("MPI error: I/O process "<<args->ioRank<<" could not receive global results. Serial/parallel sanity check will be meaningless.");
*(args->retVal) = 1;
}
}
if ( ! com->Receive( extrema_g,
n2Rows,
calcProc,
65 ) )
{
vtkGenericWarningMacro("MPI error: I/O process "<<args->ioRank<<" could not receive global results. Serial/parallel sanity check will be meaningless.");
*(args->retVal) = 1;
}
// Synchronize and stop clock
com->Barrier();
timer->StopTimer();
if ( ! com->Receive( cardsAndMeans_g,
n2Rows,
calcProc,
66 ) )
{
vtkGenericWarningMacro("MPI error: I/O process "<<args->ioRank<<" could not receive global results. Serial/parallel sanity check will be meaningless.");
*(args->retVal) = 1;
}
}
if ( com->GetLocalProcessId() == args->ioRank )
{
cout << "\n## Completed serial calculations of descriptive statistics (with assessment):\n"
<< " With partial aggregation calculated on process "
<< calcProc
<< "\n"
<< " Wall time: "
<< timer->GetElapsedTime()
<< " sec.\n";
// Synchronize and stop clock
com->Barrier();
timer->StopTimer();
for ( vtkIdType r = 0; r < nRows; ++ r )
if ( com->GetLocalProcessId() == args->ioRank )
{
cout << " "
<< outputPrimary->GetColumnName( 0 )
<< "="
<< outputPrimary->GetValue( r, 0 ).ToString()
<< " "
<< "Cardinality"
<< "="
<< cardsAndMeans_g[2 * r]
<< " "
<< "Minimum"
<< "="
<< extrema_g[2 * r]
<< " "
<< "Maximum"
<< "="
<< - extrema_g[2 * r + 1]
<< " "
<< "Mean"
<< "="
<< cardsAndMeans_g[2 * r + 1] / cardsAndMeans_g[2 * r]
<< "\n";
cout << "\n## Completed serial calculations of descriptive statistics:\n"
<< " With partial aggregation calculated on process "
<< calcProc
<< "\n"
<< " Wall time: "
<< timer->GetElapsedTime()
<< " sec.\n";
cout << " Calculated the following primary statistics:\n";
for ( vtkIdType r = 0; r < nRows; ++ r )
{
cout << " "
<< outputPrimary->GetColumnName( 0 )
<< "="
<< outputPrimary->GetValue( r, 0 ).ToString()
<< " "
<< "Cardinality"
<< "="
<< cardsAndMeans_g[2 * r]
<< " "
<< "Minimum"
<< "="
<< extrema_g[2 * r]
<< " "
<< "Maximum"
<< "="
<< - extrema_g[2 * r + 1]
<< " "
<< "Mean"
<< "="
<< cardsAndMeans_g[2 * r + 1] / cardsAndMeans_g[2 * r]
<< "\n";
}
}
}
// Clean up
delete [] cardsAndMeans_l;
delete [] cardsAndMeans_g;
// Clean up
delete [] cardsAndMeans_l;
delete [] cardsAndMeans_g;
delete [] extrema_l;
delete [] extrema_g;
delete [] extrema_l;
delete [] extrema_g;
ds->Delete();
ds->Delete();
} // if ( ! args->skipDescriptive )
else if ( com->GetLocalProcessId() == args->ioRank )
{
cout << "\n## Skipped serial calculations of descriptive statistics.\n";
}
// Now on to the actual parallel descriptive engine
// ************************** Parallel Descriptive Statistics **************************
// Skip parallel descriptive statistics if requested
if ( ! args->skipPDescriptive )
{
// Now on to the actual parallel descriptive engine
// "68-95-99.7 rule" for 1 up to numRuleVal standard deviations
int numRuleVal = 6;
// Reference values
double sigmaRuleVal[] = { 68.2689492137,
95.4499736104,
99.7300203937,
99.9936657516,
99.9999426697,
99.9999998027 };
// Tolerances
double sigmaRuleTol[] = { 1.,
.5,
.1,
.05,
.01,
.005 };
// Synchronize and start clock
com->Barrier();
timer->StartTimer();
// Instantiate a parallel descriptive statistics engine and set its input data
vtkPDescriptiveStatistics* pds = vtkPDescriptiveStatistics::New();
pds->SetInput( vtkStatisticsAlgorithm::INPUT_DATA, inputData );
// Select all columns
for ( int c = 0; c < nVariables; ++ c )
{
pds->AddColumn( columnNames[c] );
}
// Synchronize and start clock
com->Barrier();
timer->StartTimer();
// Test (in parallel) with Learn, Derive, and Assess operations turned on
pds->SetLearnOption( true );
pds->SetDeriveOption( true );
pds->SetAssessOption( true );
pds->SetTestOption( false );
pds->SignedDeviationsOff(); // Use unsigned deviations
pds->Update();
// Instantiate a parallel descriptive statistics engine and set its input data
vtkPDescriptiveStatistics* pds = vtkPDescriptiveStatistics::New();
pds->SetInput( vtkStatisticsAlgorithm::INPUT_DATA, inputData );
// Synchronize and stop clock
com->Barrier();
timer->StopTimer();
// Select all columns
for ( int c = 0; c < nVariables; ++ c )
{
pds->AddColumn( columnNames[c] );
}
// Get output data and meta tables
vtkMultiBlockDataSet* outputMetaDS = vtkMultiBlockDataSet::SafeDownCast( pds->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );
vtkTable* outputPrimary = vtkTable::SafeDownCast( outputMetaDS->GetBlock( 0 ) );
vtkTable* outputDerived = vtkTable::SafeDownCast( outputMetaDS->GetBlock( 1 ) );
vtkTable* outputData = pds->GetOutput( vtkStatisticsAlgorithm::OUTPUT_DATA );
// Test (in parallel) with Learn, Derive, and Assess options turned on
pds->SetLearnOption( true );
pds->SetDeriveOption( true );
pds->SetAssessOption( true );
pds->SetTestOption( false );
pds->SignedDeviationsOff(); // Use unsigned deviations
pds->Update();
// Synchronize and stop clock
com->Barrier();
timer->StopTimer();
// Get output data and meta tables
outputMetaDS = vtkMultiBlockDataSet::SafeDownCast( pds->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );
outputPrimary = vtkTable::SafeDownCast( outputMetaDS->GetBlock( 0 ) );
outputDerived = vtkTable::SafeDownCast( outputMetaDS->GetBlock( 1 ) );
vtkTable* outputData = pds->GetOutput( vtkStatisticsAlgorithm::OUTPUT_DATA );
if ( com->GetLocalProcessId() == args->ioRank )
{
cout << "\n## Completed parallel calculation of descriptive statistics (with assessment):\n"
<< " Total sample size: "
<< outputPrimary->GetValueByName( 0, "Cardinality" ).ToInt()
<< " \n"
<< " Wall time: "
<< timer->GetElapsedTime()
<< " sec.\n";
cout << " Calculated the following primary statistics:\n";
for ( vtkIdType r = 0; r < outputPrimary->GetNumberOfRows(); ++ r )
if ( com->GetLocalProcessId() == args->ioRank )
{
cout << " ";
for ( int i = 0; i < outputPrimary->GetNumberOfColumns(); ++ i )
cout << "\n## Completed parallel calculation of descriptive statistics (with assessment):\n"
<< " Total sample size: "
<< outputPrimary->GetValueByName( 0, "Cardinality" ).ToInt()
<< " \n"
<< " Wall time: "
<< timer->GetElapsedTime()
<< " sec.\n";
cout << " Calculated the following primary statistics:\n";
for ( vtkIdType r = 0; r < outputPrimary->GetNumberOfRows(); ++ r )
{
cout << outputPrimary->GetColumnName( i )
<< "="
<< outputPrimary->GetValue( r, i ).ToString()
<< " ";
cout << " ";
for ( int i = 0; i < outputPrimary->GetNumberOfColumns(); ++ i )
{
cout << outputPrimary->GetColumnName( i )
<< "="
<< outputPrimary->GetValue( r, i ).ToString()
<< " ";
}
cout << "\n";
}
cout << "\n";
}
cout << " Calculated the following derived statistics:\n";
for ( vtkIdType r = 0; r < outputDerived->GetNumberOfRows(); ++ r )
{
cout << " ";
for ( int i = 0; i < outputDerived->GetNumberOfColumns(); ++ i )
cout << " Calculated the following derived statistics:\n";
for ( vtkIdType r = 0; r < outputDerived->GetNumberOfRows(); ++ r )
{
cout << outputDerived->GetColumnName( i )
<< "="
<< outputDerived->GetValue( r, i ).ToString()
<< " ";
cout << " ";
for ( int i = 0; i < outputDerived->GetNumberOfColumns(); ++ i )
{
cout << outputDerived->GetColumnName( i )
<< "="
<< outputDerived->GetValue( r, i ).ToString()
<< " ";
}
cout << "\n";
}
cout << "\n";
}
}
// Verify that the DISTRIBUTED standard normal samples indeed statisfy the 68-95-99.7 rule
if ( com->GetLocalProcessId() == args->ioRank )
{
cout << "\n## Verifying whether the distributed standard normal samples satisfy the 68-95-99.7 rule:\n";
}
vtkDoubleArray* relDev[2];
relDev[0] = vtkDoubleArray::SafeDownCast(
outputData->GetColumnByName( "d(Standard Normal 0)" ) );
relDev[1] = vtkDoubleArray::SafeDownCast(
outputData->GetColumnByName( "d(Standard Normal 1)" ) );
// Verify that the DISTRIBUTED standard normal samples indeed statisfy the 68-95-99.7 rule
if ( com->GetLocalProcessId() == args->ioRank )
{
cout << "\n## Verifying whether the distributed standard normal samples satisfy the 68-95-99.7 rule:\n";
}
if ( !relDev[0] || ! relDev[1] )
{
cout << "*** Error: "
<< "Empty output column(s) on process "
<< com->GetLocalProcessId()
<< ".\n";
}
vtkDoubleArray* relDev[2];
relDev[0] = vtkDoubleArray::SafeDownCast(
outputData->GetColumnByName( "d(Standard Normal 0)" ) );
relDev[1] = vtkDoubleArray::SafeDownCast(
outputData->GetColumnByName( "d(Standard Normal 1)" ) );
// For each normal variable, count deviations of more than 1, ..., numRuleVal standard deviations from the mean
for ( int c = 0; c < nNormal; ++ c )
{
// Allocate and initialize counters
int* outsideStdv_l = new int[numRuleVal];
for ( int i = 0; i < numRuleVal; ++ i )
if ( !relDev[0] || ! relDev[1] )
{
outsideStdv_l[i] = 0;
cout << "*** Error: "
<< "Empty output column(s) on process "
<< com->GetLocalProcessId()
<< ".\n";
}
// Count outliers
double dev;
int n = outputData->GetNumberOfRows();
for ( vtkIdType r = 0; r < n; ++ r )
// For each normal variable, count deviations of more than 1, ..., numRuleVal standard deviations from the mean
for ( int c = 0; c < nNormal; ++ c )
{
dev = relDev[c]->GetValue( r );
// Allocate and initialize counters
int* outsideStdv_l = new int[numRuleVal];
for ( int i = 0; i < numRuleVal; ++ i )
{
outsideStdv_l[i] = 0;