File: vtkContingencyStatisticsGnuR.cxx

package info (click to toggle)
vtk7 7.1.1%2Bdfsg1-12
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 125,776 kB
  • sloc: cpp: 1,539,582; ansic: 106,521; python: 78,038; tcl: 47,013; xml: 8,142; yacc: 5,040; java: 4,439; perl: 3,132; lex: 1,926; sh: 1,500; makefile: 122; objc: 83
file content (110 lines) | stat: -rw-r--r-- 3,760 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
/*=========================================================================

Program:   Visualization Toolkit
Module:    vtkContingencyStatisticsGnuR.cxx

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.

=========================================================================*/
/*-------------------------------------------------------------------------
  Copyright 2011 Sandia Corporation.
  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
  the U.S. Government retains certain rights in this software.
  -------------------------------------------------------------------------*/

#include "vtkContingencyStatisticsGnuR.h"
#include "vtkRInterface.h"
#include "vtkToolkits.h"
#include "vtkStatisticsAlgorithmPrivate.h"
#include "vtkDoubleArray.h"
#include "vtkIdTypeArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkObjectFactory.h"
#include "vtkStringArray.h"
#include "vtkStdString.h"
#include "vtkTable.h"
#include "vtkVariantArray.h"

#include <map>
#include <vector>

#include <sstream>


vtkStandardNewMacro(vtkContingencyStatisticsGnuR);

vtkContingencyStatisticsGnuR::vtkContingencyStatisticsGnuR()
{

}

// ----------------------------------------------------------------------
vtkContingencyStatisticsGnuR::~vtkContingencyStatisticsGnuR()
{

}

// ----------------------------------------------------------------------
void vtkContingencyStatisticsGnuR::PrintSelf( ostream &os, vtkIndent indent )
{
  this->Superclass::PrintSelf( os, indent );
}

// ----------------------------------------------------------------------
void vtkContingencyStatisticsGnuR::CalculatePValues( vtkTable* outTab)
{
  vtkIdTypeArray* dimCol = vtkArrayDownCast<vtkIdTypeArray>( outTab->GetColumn(0) );
  vtkDoubleArray* chi2Col = vtkArrayDownCast<vtkDoubleArray>( outTab->GetColumn(1));
  vtkDoubleArray* chi2yCol = vtkArrayDownCast<vtkDoubleArray>( outTab->GetColumn(2));

  // Prepare VTK - R interface
  vtkRInterface* ri = vtkRInterface::New();

  // Use the calculated DOFs and Chi square statistics as inputs to the Chi square function
  ri->AssignVTKDataArrayToRVariable( dimCol, "d" );
  ri->AssignVTKDataArrayToRVariable( chi2Col, "chi2" );
  ri->AssignVTKDataArrayToRVariable( chi2yCol, "chi2y" );

  // Now prepare R script and calculate the p-values (in a single R script evaluation for efficiency)
  std::ostringstream rs;
  rs << "p<-c();"
     << "py<-c();"
     << "for(i in 1:"
     << dimCol->GetNumberOfTuples()
     << "){"
     << "p<-c(p,1-pchisq(chi2[i],d[i]));"
     << "py<-c(py,1-pchisq(chi2y[i],d[i]))"
     << "}";
  ri->EvalRscript( rs.str().c_str() );

  // Retrieve the p-values
  vtkDoubleArray* testChi2Col = vtkArrayDownCast<vtkDoubleArray>( ri->AssignRVariableToVTKDataArray( "p" ) );
  vtkDoubleArray* testChi2yCol = vtkArrayDownCast<vtkDoubleArray>( ri->AssignRVariableToVTKDataArray( "py" ) );
  if ( ! testChi2Col || ! testChi2yCol
       || testChi2Col->GetNumberOfTuples() != dimCol->GetNumberOfTuples()
       || testChi2yCol->GetNumberOfTuples() != dimCol->GetNumberOfTuples() )
  {
    vtkWarningMacro( "Something went wrong with the R calculations. Reported p-values will be invalid." );
    this->Superclass::CalculatePValues( outTab );
  }
  else
  {
    outTab->AddColumn( testChi2Col );
    outTab->AddColumn( testChi2yCol );
  }

  testChi2Col->SetName( "P" );
  testChi2yCol->SetName( "P Yates" );

  // Clean up
  ri->Delete();
}