00001 #include <stdio.h>
00002 #include <math.h>
00003 #include <AR/matrix.h>
00004
00005 #define VZERO 1e-16
00006 #define EPS 1e-6
00007 #define MAX_ITER 100
00008 #define xmalloc(V,T,S) { if( ((V) = (T *)malloc( sizeof(T) * (S) ))\
00009 == NULL ) {printf("malloc error\n"); exit(1);} }
00010
00011 static int EX( ARMat *input, ARVec *mean );
00012 static int CENTER( ARMat *inout, ARVec *mean );
00013 static int PCA( ARMat *input, ARMat *output, ARVec *ev );
00014 static int x_by_xt( ARMat *input, ARMat *output );
00015 static int xt_by_x( ARMat *input, ARMat *output );
00016 static int EV_create( ARMat *input, ARMat *u, ARMat *output, ARVec *ev );
00017 static int QRM( ARMat *u, ARVec *ev );
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 int arMatrixPCA( ARMat *input, ARMat *evec, ARVec *ev, ARVec *mean )
00050 {
00051 ARMat *work;
00052 double srow, sum;
00053 int row, clm;
00054 int check, rval;
00055 int i;
00056
00057 row = input->row;
00058 clm = input->clm;
00059 check = (row < clm)? row: clm;
00060 if( row < 2 || clm < 2 ) return(-1);
00061 if( evec->clm != input->clm || evec->row != check ) return(-1);
00062 if( ev->clm != check ) return(-1);
00063 if( mean->clm != input->clm ) return(-1);
00064
00065 work = arMatrixAllocDup( input );
00066 if( work == NULL ) return -1;
00067
00068 srow = sqrt((double)row);
00069 if( EX( work, mean ) < 0 ) {
00070 arMatrixFree( work );
00071 return(-1);
00072 }
00073 if( CENTER( work, mean ) < 0 ) {
00074 arMatrixFree( work );
00075 return(-1);
00076 }
00077 for(i=0; i<row*clm; i++) work->m[i] /= srow;
00078
00079 rval = PCA( work, evec, ev );
00080 arMatrixFree( work );
00081
00082 sum = 0.0;
00083 for( i = 0; i < ev->clm; i++ ) sum += ev->v[i];
00084 for( i = 0; i < ev->clm; i++ ) ev->v[i] /= sum;
00085
00086 return( rval );
00087 }
00088
00089 int arMatrixPCA2( ARMat *input, ARMat *evec, ARVec *ev )
00090 {
00091 ARMat *work;
00092
00093 double sum;
00094 int row, clm;
00095 int check, rval;
00096 int i;
00097
00098 row = input->row;
00099 clm = input->clm;
00100 check = (row < clm)? row: clm;
00101 if( row < 2 || clm < 2 ) return(-1);
00102 if( evec->clm != input->clm || evec->row != check ) return(-1);
00103 if( ev->clm != check ) return(-1);
00104
00105 work = arMatrixAllocDup( input );
00106 if( work == NULL ) return -1;
00107
00108
00109
00110
00111
00112
00113 rval = PCA( work, evec, ev );
00114 arMatrixFree( work );
00115
00116 sum = 0.0;
00117 for( i = 0; i < ev->clm; i++ ) sum += ev->v[i];
00118 for( i = 0; i < ev->clm; i++ ) ev->v[i] /= sum;
00119
00120 return( rval );
00121 }
00122
00123 static int PCA( ARMat *input, ARMat *output, ARVec *ev )
00124 {
00125 ARMat *u;
00126 double *m1, *m2;
00127 int row, clm, min;
00128 int i, j;
00129
00130 row = input->row;
00131 clm = input->clm;
00132 min = (clm < row)? clm: row;
00133 if( row < 2 || clm < 2 ) return(-1);
00134 if( output->clm != input->clm ) return(-1);
00135 if( output->row != min ) return(-1);
00136 if( ev->clm != min ) return(-1);
00137
00138 u = arMatrixAlloc( min, min );
00139 if( u->row != min || u->clm != min ) return(-1);
00140 if( row < clm ) {
00141 if( x_by_xt( input, u ) < 0 ) { arMatrixFree(u); return(-1); }
00142 }
00143 else {
00144 if( xt_by_x( input, u ) < 0 ) { arMatrixFree(u); return(-1); }
00145 }
00146
00147 if( QRM( u, ev ) < 0 ) { arMatrixFree(u); return(-1); }
00148
00149 if( row < clm ) {
00150 if( EV_create( input, u, output, ev ) < 0 ) {
00151 arMatrixFree(u);
00152 return(-1);
00153 }
00154 }
00155 else{
00156 m1 = u->m;
00157 m2 = output->m;
00158 for( i = 0; i < min; i++){
00159 if( ev->v[i] < VZERO ) break;
00160 for( j = 0; j < min; j++ ) *(m2++) = *(m1++);
00161 }
00162 for( ; i < min; i++){
00163 ev->v[i] = 0.0;
00164 for( j = 0; j < min; j++ ) *(m2++) = 0.0;
00165 }
00166 }
00167
00168 arMatrixFree(u);
00169
00170 return( 0 );
00171 }
00172
00173 static int EX( ARMat *input, ARVec *mean )
00174 {
00175 double *m, *v;
00176 int row, clm;
00177 int i, j;
00178
00179 row = input->row;
00180 clm = input->clm;
00181 if( row <= 0 || clm <= 0 ) return(-1);
00182 if( mean->clm != clm ) return(-1);
00183
00184 for( i = 0; i < clm; i++ ) mean->v[i] = 0.0;
00185
00186 m = input->m;
00187 for( i = 0; i < row; i++ ) {
00188 v = mean->v;
00189 for( j = 0; j < clm; j++ ) *(v++) += *(m++);
00190 }
00191
00192 for( i = 0; i < clm; i++ ) mean->v[i] /= row;
00193
00194 return(0);
00195 }
00196
00197 static int CENTER( ARMat *inout, ARVec *mean )
00198 {
00199 double *m, *v;
00200 int row, clm;
00201 int i, j;
00202
00203 row = inout->row;
00204 clm = inout->clm;
00205 if( mean->clm != clm ) return(-1);
00206
00207 m = inout->m;
00208 for( i = 0; i < row; i++ ) {
00209 v = mean->v;
00210 for( j = 0; j < clm; j++ ) *(m++) -= *(v++);
00211 }
00212
00213 return(0);
00214 }
00215
00216 static int x_by_xt( ARMat *input, ARMat *output )
00217 {
00218 double *in1, *in2, *out;
00219 int row, clm;
00220 int i, j, k;
00221
00222 row = input->row;
00223 clm = input->clm;
00224 if( output->row != row || output->clm != row ) return(-1);
00225
00226 out = output->m;
00227 for( i = 0; i < row; i++ ) {
00228 for( j = 0; j < row; j++ ) {
00229 if( j < i ) {
00230 *out = output->m[j*row+i];
00231 }
00232 else {
00233 in1 = &(input->m[clm*i]);
00234 in2 = &(input->m[clm*j]);
00235 *out = 0.0;
00236 for( k = 0; k < clm; k++ ) {
00237 *out += *(in1++) * *(in2++);
00238 }
00239 }
00240 out++;
00241 }
00242 }
00243
00244 return(0);
00245 }
00246
00247 static int xt_by_x( ARMat *input, ARMat *output )
00248 {
00249 double *in1, *in2, *out;
00250 int row, clm;
00251 int i, j, k;
00252
00253 row = input->row;
00254 clm = input->clm;
00255 if( output->row != clm || output->clm != clm ) return(-1);
00256
00257 out = output->m;
00258 for( i = 0; i < clm; i++ ) {
00259 for( j = 0; j < clm; j++ ) {
00260 if( j < i ) {
00261 *out = output->m[j*clm+i];
00262 }
00263 else {
00264 in1 = &(input->m[i]);
00265 in2 = &(input->m[j]);
00266 *out = 0.0;
00267 for( k = 0; k < row; k++ ) {
00268 *out += *in1 * *in2;
00269 in1 += clm;
00270 in2 += clm;
00271 }
00272 }
00273 out++;
00274 }
00275 }
00276
00277 return(0);
00278 }
00279
00280 static int EV_create( ARMat *input, ARMat *u, ARMat *output, ARVec *ev )
00281 {
00282 double *m, *m1, *m2;
00283 double sum, work;
00284 int row, clm;
00285 int i, j, k;
00286
00287 row = input->row;
00288 clm = input->clm;
00289 if( row <= 0 || clm <= 0 ) return(-1);
00290 if( u->row != row || u->clm != row ) return(-1);
00291 if( output->row != row || output->clm != clm ) return(-1);
00292 if( ev->clm != row ) return(-1);
00293
00294 m = output->m;
00295 for( i = 0; i < row; i++ ) {
00296 if( ev->v[i] < VZERO ) break;
00297 work = 1 / sqrt(fabs(ev->v[i]));
00298 for( j = 0; j < clm; j++ ) {
00299 sum = 0.0;
00300 m1 = &(u->m[i*row]);
00301 m2 = &(input->m[j]);
00302 for( k = 0; k < row; k++ ) {
00303 sum += *m1 * *m2;
00304 m1++;
00305 m2 += clm;
00306 }
00307 *(m++) = sum * work;
00308 }
00309 }
00310 for( ; i < row; i++ ) {
00311 ev->v[i] = 0.0;
00312 for( j = 0; j < clm; j++ ) *(m++) = 0.0;
00313 }
00314
00315 return(0);
00316 }
00317
00318 static int QRM( ARMat *a, ARVec *dv )
00319 {
00320 ARVec *ev, ev1;
00321 double w, t, s, x, y, c;
00322 double *v1, *v2;
00323 int dim, iter;
00324 int i, j, k, h;
00325
00326 dim = a->row;
00327 if( dim != a->clm || dim < 2 ) return(-1);
00328 if( dv->clm != dim ) return(-1);
00329
00330 ev = arVecAlloc( dim );
00331 if( ev == NULL ) return(-1);
00332
00333 ev1.clm = dim-1;
00334 ev1.v = &(ev->v[1]);
00335 if( arVecTridiagonalize( a, dv, &ev1 ) < 0 ) {
00336 arVecFree( ev );
00337 return(-1);
00338 }
00339
00340 ev->v[0] = 0.0;
00341 for( h = dim-1; h > 0; h-- ) {
00342 j = h;
00343 while(j>0 && fabs(ev->v[j]) > EPS*(fabs(dv->v[j-1])+fabs(dv->v[j]))) j--;
00344 if( j == h ) continue;
00345
00346 iter = 0;
00347 do{
00348 iter++;
00349 if( iter > MAX_ITER ) break;
00350
00351 w = (dv->v[h-1] - dv->v[h]) / 2;
00352 t = ev->v[h] * ev->v[h];
00353 s = sqrt(w*w+t);
00354 if( w < 0 ) s = -s;
00355 x = dv->v[j] - dv->v[h] + t/(w+s);
00356 y = ev->v[j+1];
00357 for( k = j; k < h; k++ ) {
00358 if( fabs(x) >= fabs(y) ) {
00359 if( fabs(x) > VZERO ) {
00360 t = -y / x;
00361 c = 1 / sqrt(t*t+1);
00362 s = t * c;
00363 }
00364 else{
00365 c = 1.0;
00366 s = 0.0;
00367 }
00368 }
00369 else{
00370 t = -x / y;
00371 s = 1.0 / sqrt(t*t+1);
00372 c = t * s;
00373 }
00374 w = dv->v[k] - dv->v[k+1];
00375 t = (w * s + 2 * c * ev->v[k+1]) * s;
00376 dv->v[k] -= t;
00377 dv->v[k+1] += t;
00378 if( k > j) ev->v[k] = c * ev->v[k] - s * y;
00379 ev->v[k+1] += s * (c * w - 2 * s * ev->v[k+1]);
00380
00381 for( i = 0; i < dim; i++ ) {
00382 x = a->m[k*dim+i];
00383 y = a->m[(k+1)*dim+i];
00384 a->m[k*dim+i] = c * x - s * y;
00385 a->m[(k+1)*dim+i] = s * x + c * y;
00386 }
00387 if( k < h-1 ) {
00388 x = ev->v[k+1];
00389 y = -s * ev->v[k+2];
00390 ev->v[k+2] *= c;
00391 }
00392 }
00393 } while(fabs(ev->v[h]) > EPS*(fabs(dv->v[h-1])+fabs(dv->v[h])));
00394 }
00395
00396 for( k = 0; k < dim-1; k++ ) {
00397 h = k;
00398 t = dv->v[h];
00399 for( i = k+1; i < dim; i++ ) {
00400 if( dv->v[i] > t ) {
00401 h = i;
00402 t = dv->v[h];
00403 }
00404 }
00405 dv->v[h] = dv->v[k];
00406 dv->v[k] = t;
00407 v1 = &(a->m[h*dim]);
00408 v2 = &(a->m[k*dim]);
00409 for( i = 0; i < dim; i++ ) {
00410 w = *v1;
00411 *v1 = *v2;
00412 *v2 = w;
00413 v1++;
00414 v2++;
00415 }
00416 }
00417
00418 arVecFree( ev );
00419 return(0);
00420 }
