/* * Copyright 2008 Sandia Corporation. * Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive * license for use of this work by or on behalf of the * U.S. Government. Redistribution and use in source and binary forms, with * or without modification, are permitted provided that this Notice and any * statement of authorship are reproduced on all copies. */ // .SECTION Thanks // Thanks to Philippe Pebay from Sandia National Laboratories // for implementing this test. #include "vtkMultiBlockDataSet.h" #include "vtkStringArray.h" #include "vtkTable.h" #include "vtkVariantArray.h" #include "vtkContingencyStatistics.h" //============================================================================= int TestContingencyStatistics( int, char *[] ) { int testStatus = 0; vtkVariant mingledData[] = { 123,456,80,"HTTP", 123,789,80,"HTTP", 123,789,80,"HTTP", 123,456,80,"HTTP", 456,123,80,"HTTP", 456,123,80,"HTTP", 456,123,8080,"HTTP", 789,123,1122,"HTTP", 456,789,80,"HTTP", 456,789,25,"SMTP", 456,789,25,"SMTP", 456,789,25,"SMTP", 456,789,25,"SMTP", 123,789,25,"SMTP", 789,123,80,"SMTP", 123,456,20,"FTP", 789,456,20,"FTP", 789,123,20,"FTP", 789,123,122,"FTP", 789,456,20,"FTP", 789,456,20,"FTP", }; int nVals = 21; vtkVariantArray* dataset0Arr = vtkVariantArray::New(); dataset0Arr->SetNumberOfComponents( 1 ); dataset0Arr->SetName( "Source" ); vtkVariantArray* dataset1Arr = vtkVariantArray::New(); dataset1Arr->SetNumberOfComponents( 1 ); dataset1Arr->SetName( "Destination" ); vtkVariantArray* dataset2Arr = vtkVariantArray::New(); dataset2Arr->SetNumberOfComponents( 1 ); dataset2Arr->SetName( "Port" ); vtkVariantArray* dataset3Arr = vtkVariantArray::New(); dataset3Arr->SetNumberOfComponents( 1 ); dataset3Arr->SetName( "Protocol" ); for ( int i = 0; i < nVals; ++ i ) { int ti = i << 2; dataset0Arr->InsertNextValue( mingledData[ti] ); dataset1Arr->InsertNextValue( mingledData[ti + 1] ); dataset2Arr->InsertNextValue( mingledData[ti + 2] ); dataset3Arr->InsertNextValue( mingledData[ti + 3] ); } vtkTable* datasetTable = vtkTable::New(); datasetTable->AddColumn( dataset0Arr ); dataset0Arr->Delete(); datasetTable->AddColumn( dataset1Arr ); dataset1Arr->Delete(); datasetTable->AddColumn( dataset2Arr ); dataset2Arr->Delete(); datasetTable->AddColumn( dataset3Arr ); dataset3Arr->Delete(); int nMetricPairs = 2; // Entropies in the summary table should normally be retrieved as follows: // column 2: H(X,Y) // column 3: H(Y|X) // column 4: H(X|Y) int iEntropies[] = { 2, 3, 4 }; int nEntropies = 3; // correct number of entropies reported in the summary table double* H = new double[nEntropies]; // Set contingency statistics algorithm and its input data port vtkContingencyStatistics* cs = vtkContingencyStatistics::New(); // First verify that absence of input does not cause trouble cout << "## Verifying that absence of input does not cause trouble... "; cs->Update(); cout << "done.\n"; // Prepare first test with data cs->SetInput( vtkStatisticsAlgorithm::INPUT_DATA, datasetTable ); vtkTable* outputData = cs->GetOutput( vtkStatisticsAlgorithm::OUTPUT_DATA ); datasetTable->Delete(); // Select Column Pair of Interest ( Learn Option ) // 1.1: a valid pair cs->AddColumnPair( "Port", "Protocol" ); // 1.2: the same valid pair, just reversed -- should thus be ignored cs->AddColumnPair( "Protocol", "Port" ); // 2: another valid pair cs->AddColumnPair( "Source", "Port" ); // 3: an invalid pair cs->AddColumnPair( "Source", "Dummy" ); // Test Learn, Derive, Assess, and Test options cs->SetLearnOption( true ); cs->SetDeriveOption( true ); cs->SetAssessOption( true ); cs->SetTestOption( true ); cs->Update(); vtkMultiBlockDataSet* outputModelDS = vtkMultiBlockDataSet::SafeDownCast( cs->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) ); vtkTable* outputSummary = vtkTable::SafeDownCast( outputModelDS->GetBlock( 0 ) ); vtkTable* outputContingency = vtkTable::SafeDownCast( outputModelDS->GetBlock( 1 ) ); vtkTable* outputTest = vtkTable::SafeDownCast( cs->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_TEST ) ); int testIntValue = 0; double testDoubleValue = 0; vtkIdType n = outputContingency->GetValueByName( 0, "Cardinality" ).ToInt(); cout << "## Calculated the following information entropies (grand total: " << n << "):\n"; testIntValue = outputSummary->GetNumberOfColumns(); if ( testIntValue != nEntropies + 2 ) { vtkGenericWarningMacro("Reported an incorrect number of columns in the summary table: " << testIntValue << " != " << nEntropies + 2 << "."); testStatus = 1; } else { // For each row in the summary table, fetch variable names and information entropies for ( vtkIdType r = 0; r < outputSummary->GetNumberOfRows(); ++ r ) { // Variable names cout << " (X,Y) = (" << outputSummary->GetValue( r, 0 ).ToString() << ", " << outputSummary->GetValue( r, 1 ).ToString() << ")"; // Information entropies for ( vtkIdType c = 0; c < nEntropies; ++ c ) { H[c] = outputSummary->GetValue( r, iEntropies[c] ).ToDouble(); cout << ", " << outputSummary->GetColumnName( iEntropies[c] ) << "=" << H[c]; } cout << "\n"; // Make sure that H(X,Y) > H(Y|X)+ H(X|Y) testDoubleValue = H[1] + H[2]; // H(Y|X)+ H(X|Y) if ( testDoubleValue > H[0] ) { vtkGenericWarningMacro("Reported inconsistent information entropies: H(X,Y) = " << H[0] << " < " << testDoubleValue << " = H(Y|X)+ H(X|Y)."); testStatus = 1; } } } cout << " where H(X,Y) = - Sum_{x,y} p(x,y) log p(x,y) and H(X|Y) = - Sum_{x,y} p(x,y) log p(x|y).\n"; cout << "\n"; cout << "## Calculated the following joint and conditional probabilities and mutual informations:\n"; testIntValue = 0; // Skip first row which contains data set cardinality vtkIdType key; for ( vtkIdType r = 1; r < outputContingency->GetNumberOfRows(); ++ r ) { key = outputContingency->GetValue( r, 0 ).ToInt(); cout << " (" << outputSummary->GetValue( key, 0 ).ToString() << "," << outputSummary->GetValue( key, 1 ).ToString() << ") = (" << outputContingency->GetValue( r, 1 ).ToString() << "," << outputContingency->GetValue( r, 2 ).ToString() << ")"; for ( vtkIdType c = 3; c < outputContingency->GetNumberOfColumns(); ++ c ) { cout << ", " << outputContingency->GetColumnName( c ) << "=" << outputContingency->GetValue( r, c ).ToDouble(); } cout << "\n"; // Update total cardinality testIntValue += outputContingency->GetValueByName( r, "Cardinality" ).ToInt(); } if ( testIntValue != nVals * nMetricPairs ) { vtkGenericWarningMacro("Reported an incorrect total cardinality: " << testIntValue << " != " << nVals * nMetricPairs << "."); testStatus = 1; } cout << "\n"; cout << "## Calculated the following marginal probabilities:\n"; testIntValue = 0; for ( unsigned int b = 2; b < outputModelDS->GetNumberOfBlocks(); ++ b ) { outputContingency = vtkTable::SafeDownCast( outputModelDS->GetBlock( b ) ); for ( vtkIdType r = 0; r < outputContingency->GetNumberOfRows(); ++ r ) { cout << " " << outputContingency->GetColumnName( 0 ) << " = " << outputContingency->GetValue( r, 0 ).ToString() << ", " << outputContingency->GetColumnName( 1 ) << "=" << outputContingency->GetValue( r, 1 ).ToDouble() << ", " << outputContingency->GetColumnName( 2 ) << "=" << outputContingency->GetValue( r, 2 ).ToDouble() << "\n"; } cout << "\n"; // Update total cardinality testIntValue += 0;//outputContingency->GetValueByName( r, "Cardinality" ).ToInt(); } // Now inspect results of the Assess option by looking for outliers key = 0; vtkStdString varX = outputSummary->GetValue( key, 0 ).ToString(); vtkStdString varY = outputSummary->GetValue( key, 1 ).ToString(); // List of columns used for outlier detection vtkStdString outlierColumn[] = { "P", "Px|y", "PMI" }; // Corresponding threshold (low) values double threshold[] = { .2, .2, .0 }; // Corresponding known number of outliers int nOutliers[] = { 4, 4, 1 }; int nOutlierTypes = 3; for ( int i = 0; i < nOutlierTypes; ++ i ) { vtkStdString colName = outlierColumn[i] + "(" + varX + "," + varY + ")"; cout << "## Found the following outliers such that " << colName << " < " << threshold[i] << ":\n"; double val; testIntValue = 0; for ( vtkIdType r = 0; r < outputData->GetNumberOfRows(); ++ r ) { val = outputData->GetValueByName( r, colName ).ToDouble(); if ( val >= threshold[i] ) { continue; } ++ testIntValue; cout << " " << outlierColumn[i] << "(" << outputData->GetValueByName( r, varX ).ToString() << "," << outputData->GetValueByName( r, varY ).ToString() << ") = " << val << "\n"; } if ( testIntValue != nOutliers[i] ) { vtkGenericWarningMacro("Reported an incorrect number of outliers: " << testIntValue << " != " << nOutliers[i] << "."); testStatus = 1; } cout << "\n"; } // Last, check some results of the Test option cout << "## Chi square statistics:\n"; // Reference values double testValues[] = { // (Port,Protocol) 10., // number of degrees of freedom 36.896, // Chi square statistic 22.35, // Chi square statistic with Yates correction #ifdef VTK_USE_GNU_R .00005899, // p-valued of Chi square statistic .01341754, // p-value of Chi square statistic with Yates correction #endif // VTK_USE_GNU_R // (Port,Source) 10., // number of degrees of freedom 17.353, // Chi square statistic 7.279, // Chi square statistic with Yates correction #ifdef VTK_USE_GNU_R .06690889, // p-valued of Chi square statistic .69886917 // p-value of Chi square statistic with Yates correction #endif // VTK_USE_GNU_R }; #ifdef VTK_USE_GNU_R double alpha = .05; vtkIdType nv = 5; #else // VTK_USE_GNU_R vtkIdType nv = 3; #endif // VTK_USE_GNU_R // Loop over Test table for ( vtkIdType r = 0; r < outputTest->GetNumberOfRows(); ++ r ) { cout << " (" << outputSummary->GetValue( r, 0 ).ToString() << "," << outputSummary->GetValue( r, 1 ).ToString() << ")"; for ( vtkIdType c = 0; c < nv; ++ c ) { double x = outputTest->GetValue( r, c ).ToDouble(); cout << ", " << outputTest->GetColumnName( c ) << "=" << x; // Verify calculated results if ( fabs ( x - testValues[r * nv + c] ) > 1.e-4 * x ) { vtkGenericWarningMacro("Incorrect " << outputTest->GetColumnName( c ) << ": " << x << " != " << testValues[r * nv + c]); testStatus = 1; } } #ifdef VTK_USE_GNU_R // Check if null hypothesis is rejected at specified significance level double p = outputTest->GetValueByName( r, "P Yates" ).ToDouble(); // Must verify that p value is valid (it is set to -1 if R has failed) if ( p > -1 && p < alpha ) { cout << ", Null hypothesis (independence) rejected at " << alpha << " significance level"; } #endif // VTK_USE_GNU_R cout << "\n"; } // Clean up delete [] H; cs->Delete(); return testStatus; }