TestRandomPMomentStatisticsMPI.cxx 28.8 KB
Newer Older
1
2
/*=========================================================================

3
4
Program:   Visualization Toolkit
Module:    TestRandomMomentStatisticsMPI.cxx
5

6
7
8
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9

10
11
12
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.
13
14
15

=========================================================================*/
/*
Philippe Pébay's avatar
Philippe Pébay committed
16
 * Copyright 2011 Sandia Corporation.
17
18
19
20
21
22
23
 * 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
Philippe Pébay's avatar
Philippe Pébay committed
24
// Thanks to Philippe Pebay, David Thompson and Janine Bennett from Sandia National Laboratories
25
26
27
28
29
30
// for implementing this test.

#include <mpi.h>

#include "vtkDescriptiveStatistics.h"
#include "vtkPDescriptiveStatistics.h"
31
32
#include "vtkCorrelativeStatistics.h"
#include "vtkPCorrelativeStatistics.h"
33
34
#include "vtkMultiCorrelativeStatistics.h"
#include "vtkPMultiCorrelativeStatistics.h"
Janine Bennett's avatar
   
Janine Bennett committed
35
#include "vtkPPCAStatistics.h"
36
37
38
39

#include "vtkDoubleArray.h"
#include "vtkMath.h"
#include "vtkMPIController.h"
40
#include "vtkMultiBlockDataSet.h"
41
42
#include "vtkStdString.h"
#include "vtkTable.h"
43
#include "vtkTimerLog.h"
44
#include "vtkVariantArray.h"
45

46
47
#include "vtksys/CommandLineArguments.hxx"

48
49
50
struct RandomSampleStatisticsArgs
{
  int nVals;
51
  bool skipDescriptive;
52
53
54
55
  bool skipPDescriptive;
  bool skipPCorrelative;
  bool skipPMultiCorrelative;
  bool skipPPCA;
56
  int* retVal;
57
  int ioRank;
58
59
};

60
// This will be called by all processes
61
void RandomSampleStatistics( vtkMultiProcessController* controller, void* arg )
62
{
63
64
  // Get test parameters
  RandomSampleStatisticsArgs* args = reinterpret_cast<RandomSampleStatisticsArgs*>( arg );
65
  *(args->retVal) = 0;
66

67
68
69
  // Get MPI communicator
  vtkMPICommunicator* com = vtkMPICommunicator::SafeDownCast( controller->GetCommunicator() );

70
  // Get local rank
71
  int myRank = com->GetLocalProcessId();
72
73

  // Seed random number generator
74
  vtkMath::RandomSeed( static_cast<int>( vtkTimerLog::GetUniversalTime() ) * ( myRank + 1 ) );
75

76
  // Generate an input table that contains samples of mutually independent random variables
77
78
79
  int nUniform = 2;
  int nNormal  = 2;
  int nVariables = nUniform + nNormal;
80

81
82
  vtkTable* inputData = vtkTable::New();
  vtkDoubleArray* doubleArray[4];
Philippe Pébay's avatar
Philippe Pébay committed
83
  vtkStdString columnNames[] = { "Standard Uniform 0",
84
85
86
                                 "Standard Uniform 1",
                                 "Standard Normal 0",
                                 "Standard Normal 1" };
Philippe Pébay's avatar
Philippe Pébay committed
87

88
89
  // Standard uniform samples
  for ( int c = 0; c < nUniform; ++ c )
90
91
92
93
94
95
    {
    doubleArray[c] = vtkDoubleArray::New();
    doubleArray[c]->SetNumberOfComponents( 1 );
    doubleArray[c]->SetName( columnNames[c] );

    double x;
96
    for ( int r = 0; r < args->nVals; ++ r )
97
98
99
      {
      x = vtkMath::Random();
      doubleArray[c]->InsertNextValue( x );
100
      }
Philippe Pébay's avatar
Philippe Pébay committed
101

102
103
104
105
106
107
108
109
110
111
112
113
    inputData->AddColumn( doubleArray[c] );
    doubleArray[c]->Delete();
    }

  // Standard normal samples
  for ( int c = nUniform; c < nVariables; ++ c )
    {
    doubleArray[c] = vtkDoubleArray::New();
    doubleArray[c]->SetNumberOfComponents( 1 );
    doubleArray[c]->SetName( columnNames[c] );

    double x;
114
    for ( int r = 0; r < args->nVals; ++ r )
115
116
117
      {
      x = vtkMath::Gaussian();
      doubleArray[c]->InsertNextValue( x );
118
      }
Philippe Pébay's avatar
Philippe Pébay committed
119

120
121
122
123
    inputData->AddColumn( doubleArray[c] );
    doubleArray[c]->Delete();
    }

124
  // Create timer to be used by all tests
125
  vtkTimerLog *timer=vtkTimerLog::New();
126

127
  // ************************** Serial descriptive Statistics **************************
Janine Bennett's avatar
   
Janine Bennett committed
128

129
  // Skip serial descriptive statistics if requested
130
  if ( ! args->skipDescriptive )
Janine Bennett's avatar
   
Janine Bennett committed
131
    {
132
133
134
    // Synchronize and start clock
    com->Barrier();
    timer->StartTimer();
Janine Bennett's avatar
   
Janine Bennett committed
135

136
137
138
    // For verification, instantiate a serial descriptive statistics engine and set its ports
    vtkDescriptiveStatistics* ds = vtkDescriptiveStatistics::New();
    ds->SetInput( vtkStatisticsAlgorithm::INPUT_DATA, inputData );
Janine Bennett's avatar
   
Janine Bennett committed
139

140
141
142
143
144
    // Select all columns
    for ( int c = 0; c < nVariables; ++ c )
      {
      ds->AddColumn( columnNames[c] );
      }
145

146
    // Test (serially) with Learn operation only (this is only to verify parallel statistics)
147
    ds->SetLearnOption( true );
148
    ds->SetDeriveOption( false );
149
150
151
    ds->SetAssessOption( false );
    ds->SetTestOption( false );
    ds->Update();
152

153
154
155
    // Get output data and meta tables
    vtkMultiBlockDataSet* outputMetaDS = vtkMultiBlockDataSet::SafeDownCast( ds->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );
    vtkTable* outputPrimary = vtkTable::SafeDownCast( outputMetaDS->GetBlock( 0 ) );
156

157
158
159
160
    // Collect (local) cardinalities, extrema, and means
    int nRows = outputPrimary->GetNumberOfRows();
    int np = com->GetNumberOfProcesses();
    int n2Rows = 2 * nRows;
161

162
163
    double* extrema_l = new double[n2Rows];
    double* extrema_g = new double[n2Rows];
164

165
166
    double* cardsAndMeans_l = new double[n2Rows];
    double* cardsAndMeans_g = new double[n2Rows];
Philippe Pébay's avatar
Philippe Pébay committed
167

168
169
170
171
172
173
    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();
174

175
176
177
178
      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();
      }
179

180
181
    // Reduce all extremal values, and gather all cardinalities and means, on process calcProc
    int calcProc = np - 1;
182

183
184
185
186
187
    com->Reduce( extrema_l,
                 extrema_g,
                 n2Rows,
                 vtkCommunicator::MIN_OP,
                 calcProc );
188

189
190
191
192
193
194
195
196
    com->Reduce( cardsAndMeans_l,
                 cardsAndMeans_g,
                 n2Rows,
                 vtkCommunicator::SUM_OP,
                 calcProc );

    // Have process calcProc calculate global cardinality and mean, and send all results to I/O process
    if ( myRank == calcProc )
197
      {
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
      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;
        }
215
216
      }

217
218
    // Have I/O process receive results from process calcProc
    if ( myRank == args->ioRank )
219
      {
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
      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;
        }

      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;
        }
237
238
      }

239
240
241
242
243
    // Synchronize and stop clock
    com->Barrier();
    timer->StopTimer();

    if ( com->GetLocalProcessId() == args->ioRank )
244
      {
245
      cout << "\n## Completed serial calculations of descriptive statistics:\n"
246
247
248
249
250
251
252
           << "   With partial aggregation calculated on process "
           << calcProc
           << "\n"
           << "   Wall time: "
           << timer->GetElapsedTime()
           << " sec.\n";

253
      cout << "   Calculated the following primary statistics:\n";
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
      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";
        }
278
279
      }

280
281
282
    // Clean up
    delete [] cardsAndMeans_l;
    delete [] cardsAndMeans_g;
283

284
285
    delete [] extrema_l;
    delete [] extrema_g;
286

287
    ds->Delete();
288
289
290
291
292
    } // if ( ! args->skipDescriptive )
  else if ( com->GetLocalProcessId() == args->ioRank )
    {
    cout << "\n## Skipped serial calculations of descriptive statistics.\n";
    }
Philippe Pébay's avatar
Philippe Pébay committed
293

294
295
296
297
298
  // ************************** Parallel Descriptive Statistics **************************

  // Skip parallel descriptive statistics if requested
  if ( ! args->skipPDescriptive )
    {
299
    // Now on to the actual parallel descriptive engine
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
    // "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 };

319

320
321
322
    // Synchronize and start clock
    com->Barrier();
    timer->StartTimer();
323

324
325
326
    // Instantiate a parallel descriptive statistics engine and set its input data
    vtkPDescriptiveStatistics* pds = vtkPDescriptiveStatistics::New();
    pds->SetInput( vtkStatisticsAlgorithm::INPUT_DATA, inputData );
327

328
329
330
331
332
    // Select all columns
    for ( int c = 0; c < nVariables; ++ c )
      {
      pds->AddColumn( columnNames[c] );
      }
Janine Bennett's avatar
   
Janine Bennett committed
333

334
    // Test (in parallel) with Learn, Derive, and Assess operations turned on
335
336
337
338
339
340
    pds->SetLearnOption( true );
    pds->SetDeriveOption( true );
    pds->SetAssessOption( true );
    pds->SetTestOption( false );
    pds->SignedDeviationsOff(); // Use unsigned deviations
    pds->Update();
Philippe Pébay's avatar
Philippe Pébay committed
341

342
343
344
    // Synchronize and stop clock
    com->Barrier();
    timer->StopTimer();
345

346
    // Get output data and meta tables
347
348
349
    vtkMultiBlockDataSet* outputMetaDS = vtkMultiBlockDataSet::SafeDownCast( pds->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );
    vtkTable* outputPrimary = vtkTable::SafeDownCast( outputMetaDS->GetBlock( 0 ) );
    vtkTable* outputDerived = vtkTable::SafeDownCast( outputMetaDS->GetBlock( 1 ) );
350
    vtkTable* outputData = pds->GetOutput( vtkStatisticsAlgorithm::OUTPUT_DATA );
351

352
    if ( com->GetLocalProcessId() == args->ioRank )
353
      {
354
355
356
357
358
359
360
361
362
363
      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 )
364
        {
365
366
367
368
369
370
371
372
373
        cout << "   ";
        for ( int i = 0; i < outputPrimary->GetNumberOfColumns(); ++ i )
          {
          cout << outputPrimary->GetColumnName( i )
               << "="
               << outputPrimary->GetValue( r, i ).ToString()
               << "  ";
          }
        cout << "\n";
374
        }
Philippe Pébay's avatar
Philippe Pébay committed
375

376
377
      cout << "   Calculated the following derived statistics:\n";
      for ( vtkIdType r = 0; r < outputDerived->GetNumberOfRows(); ++ r )
378
        {
379
380
381
382
383
384
385
386
387
        cout << "   ";
        for ( int i = 0; i < outputDerived->GetNumberOfColumns(); ++ i )
          {
          cout << outputDerived->GetColumnName( i )
               << "="
               << outputDerived->GetValue( r, i ).ToString()
               << "  ";
          }
        cout << "\n";
388
389
        }
      }
Philippe Pébay's avatar
Philippe Pébay committed
390

391
392
393
394
395
    // 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";
      }
396

397
398
399
400
401
    vtkDoubleArray* relDev[2];
    relDev[0] = vtkDoubleArray::SafeDownCast(
                                             outputData->GetColumnByName( "d(Standard Normal 0)" ) );
    relDev[1] = vtkDoubleArray::SafeDownCast(
                                             outputData->GetColumnByName( "d(Standard Normal 1)" ) );
402

403
    if ( !relDev[0] || ! relDev[1] )
404
      {
405
406
407
408
      cout << "*** Error: "
           << "Empty output column(s) on process "
           << com->GetLocalProcessId()
           << ".\n";
409
410
      }

411
412
    // For each normal variable, count deviations of more than 1, ..., numRuleVal standard deviations from the mean
    for ( int c = 0; c < nNormal; ++ c )
413
      {
414
415
416
417
418
419
      // Allocate and initialize counters
      int* outsideStdv_l = new int[numRuleVal];
      for ( int i = 0; i < numRuleVal; ++ i )
        {
        outsideStdv_l[i] = 0;
        }
420

421
422
423
424
      // Count outliers
      double dev;
      int n = outputData->GetNumberOfRows();
      for ( vtkIdType r = 0; r < n; ++ r )
425
        {
426
        dev = relDev[c]->GetValue( r );
427

428
        if ( dev >= 1. )
429
          {
430
          ++ outsideStdv_l[0];
431

432
          if ( dev >= 2. )
433
            {
434
            ++ outsideStdv_l[1];
435

436
            if ( dev >= 3. )
437
              {
438
              ++ outsideStdv_l[2];
439

440
              if ( dev >= 4. )
441
                {
442
443
444
445
446
447
448
449
450
451
452
                ++ outsideStdv_l[3];

                if ( dev >= 5. )
                  {
                  ++ outsideStdv_l[4];
                  } // if ( dev >= 5. )
                } // if ( dev >= 4. )
              } // if ( dev >= 3. )
            } // if ( dev >= 2. )
          } // if ( dev >= 1. )
        }
453

454
455
456
457
458
459
      // Sum all local counters
      int* outsideStdv_g = new int[numRuleVal];
      com->AllReduce( outsideStdv_l,
                      outsideStdv_g,
                      numRuleVal,
                      vtkCommunicator::SUM_OP );
460

461
462
      // Print out percentages of sample points within 1, ..., numRuleVal standard deviations from the mean.
      if ( com->GetLocalProcessId() == args->ioRank )
463
        {
464
465
466
467
468
469
        cout << "   "
             << outputData->GetColumnName( nUniform + c )
             << ":\n";
        for ( int i = 0; i < numRuleVal; ++ i )
          {
          double testVal = ( 1. - outsideStdv_g[i] / static_cast<double>( outputPrimary->GetValueByName( 0, "Cardinality" ).ToInt() ) ) * 100.;
470

471
472
473
474
475
          cout << "      "
               << testVal
               << "% within "
               << i + 1
               << " standard deviation(s) from the mean.\n";
476

477
478
479
480
481
482
          // Test some statistics
          if ( fabs ( testVal - sigmaRuleVal[i] ) > sigmaRuleTol[i] )
            {
            vtkGenericWarningMacro("Incorrect value.");
            *(args->retVal) = 1;
            }
483
          }
484
        }
485

486
487
488
489
      // Clean up
      delete [] outsideStdv_l;
      delete [] outsideStdv_g;
      } // for ( int c = 0; c < nNormal; ++ c )
490

491
492
    // Clean up
    pds->Delete();
493

494
495
496
497
498
    } // if ( ! args->skipPDescriptive )
  else if ( com->GetLocalProcessId() == args->ioRank )
    {
    cout << "\n## Skipped calculation of parallel descriptive statistics.\n";
    }
499

500
  // ************************** Parallel Correlative Statistics **************************
501

502
503
  // Skip parallel correlative statistics if requested
  if ( ! args->skipPCorrelative )
504
505
506
507
508
509
510
511
512
513
514
515
516
    {
    // Synchronize and start clock
    com->Barrier();
    timer->StartTimer();

    // Instantiate a parallel correlative statistics engine and set its input
    vtkPCorrelativeStatistics* pcs = vtkPCorrelativeStatistics::New();
    pcs->SetInput( vtkStatisticsAlgorithm::INPUT_DATA, inputData );

    // Select column pairs (uniform vs. uniform, normal vs. normal)
    pcs->AddColumnPair( columnNames[0], columnNames[1] );
    pcs->AddColumnPair( columnNames[2], columnNames[3] );

517
    // Test (in parallel) with Learn, Derive operations turned on
518
519
520
521
522
523
524
525
526
527
    pcs->SetLearnOption( true );
    pcs->SetDeriveOption( true );
    pcs->SetAssessOption( false );
    pcs->SetTestOption( false );
    pcs->Update();

    // Get output data and meta tables
    vtkMultiBlockDataSet* outputMetaDS = vtkMultiBlockDataSet::SafeDownCast( pcs->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );
    vtkTable* outputPrimary = vtkTable::SafeDownCast( outputMetaDS->GetBlock( 0 ) );
    vtkTable* outputDerived = vtkTable::SafeDownCast( outputMetaDS->GetBlock( 1 ) );
528

529
    // Synchronize and stop clock
530
531
    com->Barrier();
    timer->StopTimer();
532

533
    if ( com->GetLocalProcessId() == args->ioRank )
534
      {
535
536
537
538
539
540
541
542
543
544
      cout << "\n## Completed parallel calculation of correlative 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 )
545
        {
546
547
548
549
550
551
552
553
554
        cout << "   ";
        for ( int i = 0; i < outputPrimary->GetNumberOfColumns(); ++ i )
          {
          cout << outputPrimary->GetColumnName( i )
               << "="
               << outputPrimary->GetValue( r, i ).ToString()
               << "  ";
          }
        cout << "\n";
555
        }
Philippe Pébay's avatar
Philippe Pébay committed
556

557
558
      cout << "   Calculated the following derived statistics:\n";
      for ( vtkIdType r = 0; r < outputDerived->GetNumberOfRows(); ++ r )
559
        {
560
561
562
563
564
565
566
567
568
        cout << "   ";
        for ( int i = 0; i < outputDerived->GetNumberOfColumns(); ++ i )
          {
          cout << outputDerived->GetColumnName( i )
               << "="
               << outputDerived->GetValue( r, i ).ToString()
               << "  ";
          }
        cout << "\n";
569
570
        }
      }
Philippe Pébay's avatar
Philippe Pébay committed
571

572
573
    // Clean up
    pcs->Delete();
574
575
576
577
578
    } // if ( ! args->skipPCorrelative )
  else if ( com->GetLocalProcessId() == args->ioRank )
    {
    cout << "\n## Skipped calculation of parallel correlative statistics.\n";
    }
579

580
  // ************************** Parallel Multi-Correlative Statistics **************************
581

582
583
  // Skip parallel multi-correlative statistics if requested
  if ( ! args->skipPMultiCorrelative )
584
    {
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
    // Synchronize and start clock
    com->Barrier();
    timer->StartTimer();

    // Instantiate a parallel correlative statistics engine and set its ports
    vtkPMultiCorrelativeStatistics* pmcs = vtkPMultiCorrelativeStatistics::New();
    pmcs->SetInput( 0, inputData );

    // Select column pairs (uniform vs. uniform, normal vs. normal)
    pmcs->SetColumnStatus( columnNames[0], true );
    pmcs->SetColumnStatus( columnNames[1], true );
    pmcs->RequestSelectedColumns();

    pmcs->ResetAllColumnStates();
    pmcs->SetColumnStatus( columnNames[2], true );
    pmcs->SetColumnStatus( columnNames[3], true );
    pmcs->RequestSelectedColumns();

    pmcs->ResetAllColumnStates();
    pmcs->SetColumnStatus( columnNames[0], true );
    pmcs->SetColumnStatus( columnNames[1], true );
    pmcs->SetColumnStatus( columnNames[2], true );
    pmcs->SetColumnStatus( columnNames[3], true );
    pmcs->RequestSelectedColumns();

610
    // Test (in parallel) with Learn, Derive, and Assess operations turned on
611
612
613
    pmcs->SetLearnOption( true );
    pmcs->SetDeriveOption( true );
    pmcs->SetAssessOption( true );
614
    pmcs->SetTestOption( false );
615
616
617
618
619
620
621
622
623
624
    pmcs->Update();

    // Get output meta tables
    vtkMultiBlockDataSet* outputMetaDS = vtkMultiBlockDataSet::SafeDownCast( pmcs->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );

    // Synchronize and stop clock
    com->Barrier();
    timer->StopTimer();

    if ( com->GetLocalProcessId() == args->ioRank )
625
      {
626
627
628
629
630
631
632
633
634
635
636
637
638
639
      cout << "\n## Completed parallel calculation of multi-correlative statistics (with assessment):\n"
           << "   Total sample size: "
           << vtkTable::SafeDownCast( outputMetaDS->GetBlock( 0 ) )->GetValueByName( 0, "Entries").ToInt()
           << " \n"
           << "   Wall time: "
           << timer->GetElapsedTime()
           << " sec.\n";

      vtkTable* outputMeta;
      for ( unsigned int b = 1; b < outputMetaDS->GetNumberOfBlocks(); ++ b )
        {
        outputMeta = vtkTable::SafeDownCast( outputMetaDS->GetBlock( b ) );
        outputMeta->Dump();
        }
640
641
      }

642
643
    // Clean up
    pmcs->Delete();
644
645
646
647
648
    } // if ( ! args->skipPMultiCorrelative )
  else if ( com->GetLocalProcessId() == args->ioRank )
    {
    cout << "\n## Skipped calculation of parallel multi-correlative statistics.\n";
    }
Janine Bennett's avatar
   
Janine Bennett committed
649

650
  // ************************** Parallel PCA Statistics **************************
Janine Bennett's avatar
   
Janine Bennett committed
651

652
653
  // Skip parallel PCA statistics if requested
  if ( ! args->skipPPCA )
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
    {
    // Synchronize and start clock
    com->Barrier();
    timer->StartTimer();

    // Instantiate a parallel pca statistics engine and set its ports
    vtkPPCAStatistics* pcas = vtkPPCAStatistics::New();
    pcas->SetInput( vtkStatisticsAlgorithm::INPUT_DATA, inputData );

    // Select column pairs (uniform vs. uniform, normal vs. normal)
    pcas->SetColumnStatus( columnNames[0], true );
    pcas->SetColumnStatus( columnNames[1], true );
    pcas->RequestSelectedColumns();

    pcas->ResetAllColumnStates();
    pcas->SetColumnStatus( columnNames[2], true );
    pcas->SetColumnStatus( columnNames[3], true );
    pcas->RequestSelectedColumns();

    pcas->ResetAllColumnStates();
    pcas->SetColumnStatus( columnNames[0], true );
    pcas->SetColumnStatus( columnNames[1], true );
    pcas->SetColumnStatus( columnNames[2], true );
    pcas->SetColumnStatus( columnNames[3], true );
    pcas->RequestSelectedColumns();

680
    // Test (in parallel) with all operations turned on (because Test for parallel PCA allows to detect problems)
681
682
683
684
685
686
687
688
    pcas->SetLearnOption( true );
    pcas->SetDeriveOption( true );
    pcas->SetAssessOption( true );
    pcas->SetTestOption( true );
    pcas->Update();

    // Get output meta tables
    vtkMultiBlockDataSet* outputMetaDS = vtkMultiBlockDataSet::SafeDownCast( pcas->GetOutputDataObject( vtkStatisticsAlgorithm::OUTPUT_MODEL ) );
Philippe Pébay's avatar
Philippe Pébay committed
689

Janine Bennett's avatar
   
Janine Bennett committed
690
    // Synchronize and stop clock
691
692
    com->Barrier();
    timer->StopTimer();
Janine Bennett's avatar
   
Janine Bennett committed
693

694
    if ( com->GetLocalProcessId() == args->ioRank )
Janine Bennett's avatar
   
Janine Bennett committed
695
      {
696
697
698
699
700
701
702
703
704
705
706
707
708
709
      cout << "\n## Completed parallel calculation of pca statistics (with assessment):\n"
           << "   Total sample size: "
           << vtkTable::SafeDownCast( outputMetaDS->GetBlock( 0 ) )->GetValueByName( 0, "Entries").ToInt()
           << " \n"
           << "   Wall time: "
           << timer->GetElapsedTime()
           << " sec.\n";

      vtkTable* outputMeta;
      for ( unsigned int b = 1; b < outputMetaDS->GetNumberOfBlocks(); ++ b )
        {
        outputMeta = vtkTable::SafeDownCast( outputMetaDS->GetBlock( b ) );
        outputMeta->Dump();
        }
Janine Bennett's avatar
   
Janine Bennett committed
710
      }
711
712
713

    // Clean up
    pcas->Delete();
714
715
716
717
718
    } // if ( ! args->skipPPCA )
  else if ( com->GetLocalProcessId() == args->ioRank )
    {
    cout << "\n## Skipped calculation of parallel PCA statistics.\n";
    }
Janine Bennett's avatar
   
Janine Bennett committed
719
720

  // Clean up
721
  inputData->Delete();
722
  timer->Delete();
723
724
725
726
727
}

//----------------------------------------------------------------------------
int main( int argc, char** argv )
{
Philippe Pébay's avatar
Philippe Pébay committed
728
  // **************************** MPI Initialization ***************************
729
730
731
  vtkMPIController* controller = vtkMPIController::New();
  controller->Initialize( &argc, &argv );

732
733
  // If an MPI controller was not created, terminate in error.
  if ( ! controller->IsA( "vtkMPIController" ) )
734
    {
735
736
737
    vtkGenericWarningMacro("Failed to initialize a MPI controller.");
    controller->Delete();
    return 1;
Philippe Pébay's avatar
Philippe Pébay committed
738
    }
739

Philippe Pebay's avatar
Philippe Pebay committed
740
741
  vtkMPICommunicator* com = vtkMPICommunicator::SafeDownCast( controller->GetCommunicator() );

Philippe Pébay's avatar
Philippe Pébay committed
742
  // ************************** Find an I/O node ********************************
743
744
745
746
  int* ioPtr;
  int ioRank;
  int flag;

Philippe Pébay's avatar
Philippe Pébay committed
747
  MPI_Attr_get( MPI_COMM_WORLD,
748
749
750
751
752
753
                MPI_IO,
                &ioPtr,
                &flag );

  if ( ( ! flag ) || ( *ioPtr == MPI_PROC_NULL ) )
    {
Philippe Pebay's avatar
Philippe Pebay committed
754
    // Getting MPI attributes did not return any I/O node found.
755
756
757
758
759
760
761
    ioRank = MPI_PROC_NULL;
    vtkGenericWarningMacro("No MPI I/O nodes found.");

    // As no I/O node was found, we need an unambiguous way to report the problem.
    // This is the only case when a testValue of -1 will be returned
    controller->Finalize();
    controller->Delete();
Philippe Pébay's avatar
Philippe Pébay committed
762

763
764
    return -1;
    }
Philippe Pébay's avatar
Philippe Pébay committed
765
  else
766
    {
Philippe Pebay's avatar
Philippe Pebay committed
767
    if ( *ioPtr == MPI_ANY_SOURCE )
768
769
770
771
772
773
774
      {
      // Anyone can do the I/O trick--just pick node 0.
      ioRank = 0;
      }
    else
      {
      // Only some nodes can do I/O. Make sure everyone agrees on the choice (min).
Philippe Pebay's avatar
Philippe Pebay committed
775
776
777
778
      com->AllReduce( ioPtr,
                      &ioRank,
                      1,
                      vtkCommunicator::MIN_OP );
779
780
781
      }
    }

782
  if ( com->GetLocalProcessId() == ioRank )
783
    {
784
    cout << "\n# Process "
785
         << ioRank
786
         << " will be the I/O node.\n";
787
    }
Philippe Pébay's avatar
Philippe Pébay committed
788

789
790
  // Check how many processes have been made available
  int numProcs = controller->GetNumberOfProcesses();
791
  if ( controller->GetLocalProcessId() == ioRank )
792
    {
793
    cout << "\n# Running test with "
794
795
796
797
         << numProcs
         << " processes...\n";
    }

798
  // **************************** Parse command line ***************************
799
  // Set default argument values
800
  int nVals = 100000;
801
  bool skipDescriptive = false;
802
803
804
805
  bool skipPDescriptive = false;
  bool skipPCorrelative = false;
  bool skipPMultiCorrelative = false;
  bool skipPPCA = false;
806

807
808
809
810
811
  // Initialize command line argument parser
  vtksys::CommandLineArguments clArgs;
  clArgs.Initialize( argc, argv );
  clArgs.StoreUnusedArguments( false );

812
  // Parse per-process cardinality of each pseudo-random sample
813
814
815
816
  clArgs.AddArgument("--n-per-proc",
                     vtksys::CommandLineArguments::SPACE_ARGUMENT,
                     &nVals, "Per-process cardinality of each pseudo-random sample");

817
818
819
820
821
822
823
  // Parse whether serial descriptive statistics should be skipped (for faster testing)
  clArgs.AddArgument("--skip-Descriptive",
                     vtksys::CommandLineArguments::NO_ARGUMENT,
                     &skipDescriptive, "Skip serial descriptive statistics");

  // Parse whether parallel descriptive statistics should be skipped (for faster testing)
  clArgs.AddArgument("--skip-PDescriptive",
824
                     vtksys::CommandLineArguments::NO_ARGUMENT,
825
                     &skipPDescriptive, "Skip parallel descriptive statistics");
826

827
828
  // Parse whether parallel correlative statistics should be skipped (for faster testing)
  clArgs.AddArgument("--skip-PCorrelative",
829
                     vtksys::CommandLineArguments::NO_ARGUMENT,
830
                     &skipPCorrelative, "Skip parallel correlative statistics");
831

832
833
  // Parse whether parallel multi-correlative statistics should be skipped (for faster testing)
  clArgs.AddArgument("--skip-PMultiCorrelative",
834
                     vtksys::CommandLineArguments::NO_ARGUMENT,
835
                     &skipPMultiCorrelative, "Skip parallel multi-correlative statistics");
836

837
838
  // Parse whether parallel PCA statistics should be skipped (for faster testing)
  clArgs.AddArgument("--skip-PPCA",
839
                     vtksys::CommandLineArguments::NO_ARGUMENT,
840
                     &skipPPCA, "Skip parallel PCA statistics");
841

842
843
844
845
846
847
848
  // If incorrect arguments were provided, terminate in error.
  if ( ! clArgs.Parse() )
    {
    vtkGenericWarningMacro("Incorrect input data arguments were provided.");
    return 1;
    }

Philippe Pébay's avatar
Philippe Pébay committed
849
  // ************************** Initialize test *********************************
850
  // Parameters for regression test.
851
  int testValue = 0;
852
  RandomSampleStatisticsArgs args;
853
  args.nVals = nVals;
854
  args.skipDescriptive = skipDescriptive;
855
856
857
858
  args.skipPDescriptive = skipPDescriptive;
  args.skipPCorrelative = skipPCorrelative;
  args.skipPMultiCorrelative = skipPMultiCorrelative;
  args.skipPPCA = skipPPCA;
859
  args.retVal = &testValue;
860
  args.ioRank = ioRank;
861

862
  // Execute the function named "process" on both processes
863
  controller->SetSingleMethod( RandomSampleStatistics, &args );
864
  controller->SingleMethodExecute();
865

866
  // Clean up and exit
867
  if ( com->GetLocalProcessId() == ioRank )
868
869
870
871
    {
    cout << "\n# Test completed.\n\n";
    }

872
873
  controller->Finalize();
  controller->Delete();
Philippe Pébay's avatar
Philippe Pébay committed
874

875
876
  return testValue;
}