1 // .SECTION Thanks
2 // This test was implemented by Philippe Pebay, Kitware SAS 2012
3
4 #include "vtkAutoCorrelativeStatistics.h"
5 #include "vtkDataObjectCollection.h"
6 #include "vtkDoubleArray.h"
7 #include "vtkIdTypeArray.h"
8 #include "vtkInformation.h"
9 #include "vtkMath.h"
10 #include "vtkMultiBlockDataSet.h"
11 #include "vtkStringArray.h"
12 #include "vtkTable.h"
13 #include "vtkTimerLog.h"
14 #include "vtkVariantArray.h"
15
16 //=============================================================================
TestAutoCorrelativeStatistics(int,char * [])17 int TestAutoCorrelativeStatistics( int, char *[] )
18 {
19 int testStatus = 0;
20
21 // ************** Test with 2 columns of input data **************
22
23 // Input data
24 double mingledData[] =
25 {
26 46,
27 45,
28 47,
29 49,
30 46,
31 47,
32 46,
33 46,
34 47,
35 46,
36 47,
37 49,
38 49,
39 49,
40 47,
41 45,
42 50,
43 50,
44 46,
45 46,
46 51,
47 50,
48 48,
49 48,
50 52,
51 54,
52 48,
53 47,
54 52,
55 52,
56 49,
57 49,
58 53,
59 54,
60 50,
61 50,
62 53,
63 54,
64 50,
65 52,
66 53,
67 53,
68 50,
69 51,
70 54,
71 54,
72 49,
73 49,
74 52,
75 52,
76 50,
77 51,
78 52,
79 52,
80 49,
81 47,
82 48,
83 48,
84 48,
85 50,
86 46,
87 48,
88 47,
89 47,
90 };
91
92 // Test with entire data set
93 int nVals1 = 32;
94
95 vtkDoubleArray* dataset1Arr = vtkDoubleArray::New();
96 dataset1Arr->SetNumberOfComponents( 1 );
97 dataset1Arr->SetName( "Metric 0" );
98
99 vtkDoubleArray* dataset2Arr = vtkDoubleArray::New();
100 dataset2Arr->SetNumberOfComponents( 1 );
101 dataset2Arr->SetName( "Metric 1" );
102
103 for ( int i = 0; i < nVals1; ++ i )
104 {
105 int ti = i << 1;
106 dataset1Arr->InsertNextValue( mingledData[ti] );
107 dataset2Arr->InsertNextValue( mingledData[ti + 1] );
108 }
109
110 // Create input data table
111 vtkTable* datasetTable1 = vtkTable::New();
112 datasetTable1->AddColumn( dataset1Arr );
113 dataset1Arr->Delete();
114 datasetTable1->AddColumn( dataset2Arr );
115 dataset2Arr->Delete();
116
117 // Create input parameter table for the stationary case
118 vtkIdTypeArray* timeLags = vtkIdTypeArray::New();
119 timeLags->SetName( "Time Lags" );
120 timeLags->SetNumberOfTuples( 1 );
121 timeLags->SetValue( 0, 0 );
122 vtkTable* paramTable = vtkTable::New();
123 paramTable->AddColumn( timeLags );
124 timeLags->Delete();
125
126 // Columns of interest
127 int nMetrics1 = 2;
128 vtkStdString columns1[] =
129 {
130 "Metric 1",
131 "Metric 0"
132 };
133
134 // Reference values
135 // Means for metrics 0 and 1 respectively
136 double meansXs1[] = { 49.21875 , 49.5 };
137
138 // Standard deviations for metrics 0 and 1, respectively
139 double varsXs1[] = { 5.9828629, 7.548397 };
140
141 // Set autocorrelative statistics algorithm and its input data port
142 vtkAutoCorrelativeStatistics* as1 = vtkAutoCorrelativeStatistics::New();
143
144 // First verify that absence of input does not cause trouble
145 cout << "\n## Verifying that absence of input does not cause trouble... ";
146 as1->Update();
147 cout << "done.\n";
148
149 // Prepare first test with data
150 as1->SetInputData( vtkStatisticsAlgorithm::INPUT_DATA, datasetTable1 );
151 datasetTable1->Delete();
152
153 // Select columns of interest
154 for ( int i = 0; i < nMetrics1; ++ i )
155 {
156 as1->AddColumn( columns1[i] );
157 }
158
159 // Set spatial cardinality
160 as1->SetSliceCardinality( nVals1 );
161
162 // Set parameters for autocorrelation of whole data set with respect to itself
163 as1->SetInputData( vtkStatisticsAlgorithm::LEARN_PARAMETERS, paramTable );
164
165 // Test Learn and Derive options
166 as1->SetLearnOption( true );
167 as1->SetDeriveOption( true );
168 as1->SetAssessOption( false );
169 as1->SetTestOption( false );
170 as1->Update();
171
172 // Get output model tables
173 vtkMultiBlockDataSet* outputModelAS1 = vtkMultiBlockDataSet::SafeDownCast( as1->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );
174
175 cout << "\n## Calculated the following statistics for first data set:\n";
176 for ( unsigned b = 0; b < outputModelAS1->GetNumberOfBlocks(); ++ b )
177 {
178 vtkStdString varName = outputModelAS1->GetMetaData( b )->Get( vtkCompositeDataSet::NAME() );
179
180 vtkTable* modelTab = vtkTable::SafeDownCast( outputModelAS1->GetBlock( b ) );
181 if ( varName == "Autocorrelation FFT" )
182 {
183 if ( modelTab->GetNumberOfRows() )
184 {
185 cout << "\n Autocorrelation FFT:\n";
186 modelTab->Dump();
187 continue;
188 }
189 }
190
191 cout << " Variable="
192 << varName
193 << "\n";
194
195 cout << " ";
196 for ( int i = 0; i < modelTab->GetNumberOfColumns(); ++ i )
197 {
198 cout << modelTab->GetColumnName( i )
199 << "="
200 << modelTab->GetValue( 0, i ).ToString()
201 << " ";
202 }
203
204 // Verify some of the calculated statistics
205 if ( fabs ( modelTab->GetValueByName( 0, "Mean Xs" ).ToDouble() - meansXs1[b] ) > 1.e-6 )
206 {
207 vtkGenericWarningMacro("Incorrect mean for Xs");
208 testStatus = 1;
209 }
210
211 if ( fabs ( modelTab->GetValueByName( 0, "Variance Xs" ).ToDouble() - varsXs1[b] ) > 1.e-5 )
212 {
213 vtkGenericWarningMacro("Incorrect variance for Xs");
214 testStatus = 1;
215 }
216
217 if ( fabs ( modelTab->GetValueByName( 0, "Autocorrelation" ).ToDouble() - 1. ) > 1.e-6 )
218 {
219 vtkGenericWarningMacro("Incorrect autocorrelation");
220 testStatus = 1;
221 }
222
223 cout << "\n";
224 }
225
226 // Test with a slight variation of initial data set (to test model aggregation)
227 int nVals2 = 32;
228
229 vtkDoubleArray* dataset4Arr = vtkDoubleArray::New();
230 dataset4Arr->SetNumberOfComponents( 1 );
231 dataset4Arr->SetName( "Metric 0" );
232
233 vtkDoubleArray* dataset5Arr = vtkDoubleArray::New();
234 dataset5Arr->SetNumberOfComponents( 1 );
235 dataset5Arr->SetName( "Metric 1" );
236
237 for ( int i = 0; i < nVals2; ++ i )
238 {
239 int ti = i << 1;
240 dataset4Arr->InsertNextValue( mingledData[ti] + 1. );
241 dataset5Arr->InsertNextValue( mingledData[ti + 1] );
242 }
243
244 vtkTable* datasetTable2 = vtkTable::New();
245 datasetTable2->AddColumn( dataset4Arr );
246 dataset4Arr->Delete();
247 datasetTable2->AddColumn( dataset5Arr );
248 dataset5Arr->Delete();
249
250 // Set auto-correlative statistics algorithm and its input data port
251 vtkAutoCorrelativeStatistics* as2 = vtkAutoCorrelativeStatistics::New();
252 as2->SetInputData( vtkStatisticsAlgorithm::INPUT_DATA, datasetTable2 );
253
254 // Select columns of interest
255 for ( int i = 0; i < nMetrics1; ++ i )
256 {
257 as2->AddColumn( columns1[i] );
258 }
259
260 // Set spatial cardinality
261 as2->SetSliceCardinality( nVals2 );
262
263 // Set parameters for autocorrelation of whole data set with respect to itself
264 as2->SetInputData( vtkStatisticsAlgorithm::LEARN_PARAMETERS, paramTable );
265
266 // Update with Learn option only
267 as2->SetLearnOption( true );
268 as2->SetDeriveOption( false );
269 as2->SetTestOption( false );
270 as2->SetAssessOption( false );
271 as2->Update();
272
273 // Get output meta tables
274 vtkMultiBlockDataSet* outputModelAS2 = vtkMultiBlockDataSet::SafeDownCast( as2->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );
275
276 cout << "\n## Calculated the following statistics for second data set:\n";
277 for ( unsigned b = 0; b < outputModelAS2->GetNumberOfBlocks(); ++ b )
278 {
279 vtkStdString varName = outputModelAS2->GetMetaData( b )->Get( vtkCompositeDataSet::NAME() );
280
281 vtkTable* modelTab = vtkTable::SafeDownCast( outputModelAS2->GetBlock( b ) );
282 if ( varName == "Autocorrelation FFT" )
283 {
284 if ( modelTab->GetNumberOfRows() )
285 {
286 cout << "\n Autocorrelation FFT:\n";
287 modelTab->Dump();
288 continue;
289 }
290 }
291
292 cout << "\n Variable="
293 << varName
294 << "\n";
295
296 cout << " ";
297 for ( int i = 0; i < modelTab->GetNumberOfColumns(); ++ i )
298 {
299 cout << modelTab->GetColumnName( i )
300 << "="
301 << modelTab->GetValue( 0, i ).ToString()
302 << " ";
303 }
304
305 cout << "\n";
306 }
307
308 // Test model aggregation by adding new data to engine which already has a model
309 as1->SetInputData( vtkStatisticsAlgorithm::INPUT_DATA, datasetTable2 );
310 vtkMultiBlockDataSet* model = vtkMultiBlockDataSet::New();
311 model->ShallowCopy( outputModelAS1 );
312 as1->SetInputData( vtkStatisticsAlgorithm::INPUT_MODEL, model );
313
314 // Clean up
315 model->Delete();
316 datasetTable2->Delete();
317 as2->Delete();
318
319 // Update with Learn and Derive options only
320 as1->SetLearnOption( true );
321 as1->SetDeriveOption( true );
322 as1->SetTestOption( false );
323 as1->SetAssessOption( false );
324 as1->Update();
325
326 // Updated reference values
327 // Means and variances for metrics 0 and 1, respectively
328 double meansXs0[] = { 49.71875 , 49.5 };
329 double varsXs0[] = { 6.1418651 , 7.548397 * 62. / 63. };
330
331 // Get output meta tables
332 outputModelAS1 = vtkMultiBlockDataSet::SafeDownCast( as1->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );
333
334 cout << "\n## Calculated the following statistics for aggregated (first + second) data set:\n";
335 for ( unsigned b = 0; b < outputModelAS1->GetNumberOfBlocks(); ++ b )
336 {
337 vtkStdString varName = outputModelAS1->GetMetaData( b )->Get( vtkCompositeDataSet::NAME() );
338
339 vtkTable* modelTab = vtkTable::SafeDownCast( outputModelAS1->GetBlock( b ) );
340 if ( varName == "Autocorrelation FFT" )
341 {
342 if ( modelTab->GetNumberOfRows() )
343 {
344 cout << "\n Autocorrelation FFT:\n";
345 modelTab->Dump();
346 continue;
347 }
348 }
349
350 cout << "\n Variable="
351 << varName
352 << "\n";
353
354 cout << " ";
355 for ( int i = 0; i < modelTab->GetNumberOfColumns(); ++ i )
356 {
357 cout << modelTab->GetColumnName( i )
358 << "="
359 << modelTab->GetValue( 0, i ).ToString()
360 << " ";
361 }
362
363 // Verify some of the calculated statistics
364 if ( fabs ( modelTab->GetValueByName( 0, "Mean Xs" ).ToDouble() - meansXs0[b] ) > 1.e-6 )
365 {
366 vtkGenericWarningMacro("Incorrect mean for Xs");
367 testStatus = 1;
368 }
369
370 if ( fabs ( modelTab->GetValueByName( 0, "Variance Xs" ).ToDouble() - varsXs0[b] ) > 1.e-5 )
371 {
372 vtkGenericWarningMacro("Incorrect variance for Xs");
373 testStatus = 1;
374 }
375
376 cout << "\n";
377 }
378
379 // Clean up
380 as1->Delete();
381
382 // ************** Test with 2 columns of synthetic data **************
383
384 // Space and time parameters
385 vtkIdType nSteps = 2;
386 vtkIdType cardSlice = 1000;
387 vtkIdType cardTotal = nSteps * cardSlice;
388
389 // Expand parameter table to contain all steps
390 vtkVariantArray* row = vtkVariantArray::New();
391 row->SetNumberOfValues( 1 );
392 for ( vtkIdType p = 1; p < nSteps; ++ p )
393 {
394 row->SetValue( 0, p );
395 paramTable->InsertNextRow( row );
396 }
397 row->Delete();
398
399 vtkDoubleArray* lineArr = vtkDoubleArray::New();
400 lineArr->SetNumberOfComponents( 1 );
401 lineArr->SetName( "Line" );
402
403 vtkDoubleArray* vArr = vtkDoubleArray::New();
404 vArr->SetNumberOfComponents( 1 );
405 vArr->SetName( "V" );
406
407 vtkDoubleArray* circleArr = vtkDoubleArray::New();
408 circleArr->SetNumberOfComponents( 1 );
409 circleArr->SetName( "Circle" );
410
411 // Fill data columns
412 vtkIdType midPoint = cardTotal >> 1;
413 double dAlpha = ( 2.0 * vtkMath::Pi() ) / cardSlice;
414 for ( int i = 0; i < cardTotal; ++ i )
415 {
416 lineArr->InsertNextValue( i );
417 if ( i < midPoint )
418 {
419 vArr->InsertNextValue( cardTotal - i );
420 circleArr->InsertNextValue( cos( i * dAlpha ) );
421 }
422 else
423 {
424 vArr->InsertNextValue( i );
425 circleArr->InsertNextValue( sin( i * dAlpha ) );
426 }
427 }
428
429
430 // Create input data table
431 vtkTable* datasetTable3 = vtkTable::New();
432 datasetTable3->AddColumn( lineArr );
433 lineArr->Delete();
434 datasetTable3->AddColumn( vArr );
435 vArr->Delete();
436 datasetTable3->AddColumn( circleArr );
437 circleArr->Delete();
438
439 // Columns of interest
440 int nMetrics2 = 3;
441 vtkStdString columns2[] =
442 {
443 "Line",
444 "V",
445 "Circle"
446 };
447
448 // Reference values
449 double halfNm1 = .5 * ( cardSlice - 1 );
450
451 // Means of Xs for circle, line, and v-shaped variables respectively
452 double meansXt3[] = { 0., 0.,
453 halfNm1, halfNm1 + cardSlice,
454 cardTotal - halfNm1, cardTotal - halfNm1 - 1. };
455
456 // Autocorrelation values for circle, line, and v-shaped variables respectively
457 double autocorr3[] = { 1., 0.,
458 1., 1.,
459 1., -1. };
460
461 // Prepare autocorrelative statistics algorithm and its input data port
462 vtkAutoCorrelativeStatistics* as3 = vtkAutoCorrelativeStatistics::New();
463 as3->SetInputData( vtkStatisticsAlgorithm::INPUT_DATA, datasetTable3 );
464 datasetTable3->Delete();
465
466 // Select Columns of Interest
467 for ( int i = 0; i < nMetrics2; ++ i )
468 {
469 as3->AddColumn( columns2[i] );
470 }
471
472 // Set spatial cardinality
473 as3->SetSliceCardinality( cardSlice );
474
475 // Set autocorrelation parameters for first slice against slice following midpoint
476 as3->SetInputData( vtkStatisticsAlgorithm::LEARN_PARAMETERS, paramTable );
477
478 // Test Learn, and Derive options
479 as3->SetLearnOption( true );
480 as3->SetDeriveOption( true );
481 as3->SetAssessOption( false );
482 as3->SetTestOption( false );
483 as3->Update();
484
485 // Get output data and meta tables
486 vtkMultiBlockDataSet* outputModelAS3 = vtkMultiBlockDataSet::SafeDownCast( as3->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );
487
488 cout << "\n## Calculated the following statistics for third data set:\n";
489 for ( unsigned b = 0; b < outputModelAS3->GetNumberOfBlocks(); ++ b )
490 {
491 vtkStdString varName = outputModelAS3->GetMetaData( b )->Get( vtkCompositeDataSet::NAME() );
492
493 vtkTable* modelTab = vtkTable::SafeDownCast( outputModelAS3->GetBlock( b ) );
494 if ( varName == "Autocorrelation FFT" )
495 {
496 if ( modelTab->GetNumberOfRows() )
497 {
498 cout << "\n Autocorrelation FFT:\n";
499 modelTab->Dump();
500 continue;
501 }
502 }
503
504 cout << "\n Variable="
505 << varName
506 << "\n";
507
508 for ( int r = 0; r < modelTab->GetNumberOfRows(); ++ r )
509 {
510 cout << " ";
511 for ( int i = 0; i < modelTab->GetNumberOfColumns(); ++ i )
512 {
513 cout << modelTab->GetColumnName( i )
514 << "="
515 << modelTab->GetValue( r, i ).ToString()
516 << " ";
517 }
518
519 // Verify some of the calculated statistics
520 int idx = nSteps * b + r;
521 if ( fabs ( modelTab->GetValueByName( r, "Mean Xt" ).ToDouble() - meansXt3[idx] ) > 1.e-6 )
522 {
523 vtkGenericWarningMacro("Incorrect mean for Xt");
524 testStatus = 1;
525 }
526
527 if ( fabs ( modelTab->GetValueByName( r, "Autocorrelation" ).ToDouble() - autocorr3[idx] ) > 1.e-6 )
528 {
529 vtkGenericWarningMacro("Incorrect autocorrelation "<<autocorr3[idx]);
530 testStatus = 1;
531 }
532
533 cout << "\n";
534 } // i
535 } // r
536
537 // Clean up
538 as3->Delete();
539 paramTable->Delete();
540
541 return testStatus;
542 }
543