26 #include "jama/jama_eig.h" 27 #include "jama/jama_lu.h" 33 PrincipalComponentAnalysis::PrincipalComponentAnalysis(
const int n) {
36 p_statistics.resize(n * n);
37 for(
int i = 0; i < n; i++) {
38 for(
int j = 0; j <= i; j++) {
40 p_statistics[n*j+i] = p_statistics[n*i+j];
44 p_hasTransform =
false;
48 PrincipalComponentAnalysis::PrincipalComponentAnalysis(TNT::Array2D<double> transform) {
49 if(transform.dim1() != transform.dim2()) {
50 std::string m =
"Illegal transform matrix";
54 p_dimensions = transform.dim1();
55 p_transform = transform;
58 p_hasTransform =
true;
64 void PrincipalComponentAnalysis::AddData(
const double *data,
const unsigned int count) {
68 std::string m =
"Cannot add data to a PCA that has a defined transform matrix";
69 throw IException(IException::Programmer, m,
_FILEINFO_);
73 for(
int i = 0; i < p_dimensions; i++) {
74 for(
int j = 0; j <= i; j++) {
75 p_statistics[p_dimensions*i+j]->AddData(&data[count*i], &data[count*j], count);
81 void PrincipalComponentAnalysis::ComputeTransform() {
83 std::string m =
"This PCA already has a computed transform";
84 throw IException(IException::Programmer, m,
_FILEINFO_);
87 TNT::Array2D<double> C(p_dimensions, p_dimensions);
88 for(
int i = 0; i < p_dimensions; i++) {
89 for(
int j = 0; j < p_dimensions; j++) {
90 C[i][j] = (p_statistics[p_dimensions*i+j])->Correlation();
94 JAMA::Eigenvalue<double>
E(C);
95 TNT::Array2D<double> D, V;
99 p_transform = TNT::Array2D<double>(V.dim1(), V.dim2());
105 for(
int i = 0; i < V.dim1(); i++) {
106 for(
int j = 0; j < V.dim2(); j++) {
107 p_transform[i][j] = V[i][V.dim2()-j-1];
113 p_hasTransform =
true;
117 void PrincipalComponentAnalysis::ComputeInverse() {
118 TNT::Array2D<double> id(p_transform.dim1(), p_transform.dim2(), 0.0);
119 for(
int i = 0; i < p_transform.dim1(); i++)
id[i][i] = 1;
121 JAMA::LU<double> lu(p_transform);
122 if(lu.det() == 0.0) {
123 std::string m =
"Cannot take the inverse of the transform matrix";
124 throw IException(IException::Programmer, m,
_FILEINFO_);
127 p_inverse = lu.solve(
id);
131 TNT::Array2D<double> PrincipalComponentAnalysis::Transform(TNT::Array2D<double> data) {
132 if(data.dim1() != 1 || data.dim2() != p_dimensions) {
133 QString m =
"Transform input must be of dimension 1 x " + QString::number(p_dimensions);
134 throw IException(IException::Programmer, m,
_FILEINFO_);
138 return TNT::matmult(data, p_transform);
142 TNT::Array2D<double> PrincipalComponentAnalysis::Inverse(TNT::Array2D<double> data) {
143 if(data.dim1() != 1 || data.dim2() != p_dimensions) {
144 QString m =
"Transform input must be of dimension 1 x " + QString::number(p_dimensions);
145 throw IException(IException::Programmer, m,
_FILEINFO_);
149 return TNT::matmult(data, p_inverse);
Namespace for the standard library.
Container of multivariate statistics.
#define _FILEINFO_
Macro for the filename and line number.
const double E
Sets some basic constants for use in ISIS programming.
Namespace for ISIS/Bullet specific routines.