00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 template<typename eT>
00031 inline
00032 void
00033 op_princomp::direct_princomp
00034 (
00035 Mat<eT>& coeff_out,
00036 Mat<eT>& score_out,
00037 Col<eT>& latent_out,
00038 Col<eT>& tsquared_out,
00039 const Mat<eT>& in
00040 )
00041 {
00042 arma_extra_debug_sigprint();
00043
00044 const u32 n_rows = in.n_rows;
00045 const u32 n_cols = in.n_cols;
00046
00047 if(n_rows > 1)
00048 {
00049
00050 score_out = in - repmat(mean(in), n_rows, 1);
00051
00052
00053 Mat<eT> U;
00054 Col<eT> s;
00055
00056 const bool svd_ok = svd(U,s,coeff_out,score_out);
00057
00058 if(svd_ok == false)
00059 {
00060 arma_print("princomp(): singular value decomposition failed");
00061
00062 coeff_out.reset();
00063 score_out.reset();
00064 latent_out.reset();
00065 tsquared_out.reset();
00066
00067 return;
00068 }
00069
00070
00071
00072
00073
00074 s /= std::sqrt(n_rows - 1);
00075
00076
00077 score_out *= coeff_out;
00078
00079 if(n_rows <= n_cols)
00080 {
00081 score_out.cols(n_rows-1,n_cols-1).zeros();
00082
00083
00084 Col<eT> s_tmp(n_cols);
00085 s_tmp.zeros();
00086
00087 s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
00088 s = s_tmp;
00089
00090
00091 s_tmp.rows(0,n_rows-2) = eT(1) / s_tmp.rows(0,n_rows-2);
00092
00093 const Mat<eT> S = score_out * diagmat(Col<eT>(s_tmp));
00094 tsquared_out = sum(S%S,1);
00095 }
00096 else
00097 {
00098
00099 const Mat<eT> S = score_out * diagmat(Col<eT>( eT(1) / s));
00100 tsquared_out = sum(S%S,1);
00101 }
00102
00103
00104 latent_out = s%s;
00105 }
00106 else
00107 {
00108 if(n_rows == 1)
00109 {
00110 coeff_out = eye< Mat<eT> >(n_cols, n_cols);
00111
00112 score_out.copy_size(in);
00113 score_out.zeros();
00114
00115 latent_out.set_size(n_cols);
00116 latent_out.zeros();
00117
00118 tsquared_out.set_size(1);
00119 tsquared_out.zeros();
00120 }
00121 else
00122 {
00123 coeff_out.reset();
00124 score_out.reset();
00125 latent_out.reset();
00126 tsquared_out.reset();
00127 }
00128 }
00129
00130 }
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 template<typename eT>
00141 inline
00142 void
00143 op_princomp::direct_princomp
00144 (
00145 Mat<eT>& coeff_out,
00146 Mat<eT>& score_out,
00147 Col<eT>& latent_out,
00148 const Mat<eT>& in
00149 )
00150 {
00151 arma_extra_debug_sigprint();
00152
00153 const u32 n_rows = in.n_rows;
00154 const u32 n_cols = in.n_cols;
00155
00156 if(n_rows > 1)
00157 {
00158
00159 score_out = in - repmat(mean(in), n_rows, 1);
00160
00161
00162 Mat<eT> U;
00163 Col<eT> s;
00164
00165 const bool svd_ok = svd(U,s,coeff_out,score_out);
00166
00167 if(svd_ok == false)
00168 {
00169 arma_print("princomp(): singular value decomposition failed");
00170
00171 coeff_out.reset();
00172 score_out.reset();
00173 latent_out.reset();
00174
00175 return;
00176 }
00177
00178
00179
00180
00181
00182 s /= std::sqrt(n_rows - 1);
00183
00184
00185 score_out *= coeff_out;
00186
00187 if(n_rows <= n_cols)
00188 {
00189 score_out.cols(n_rows-1,n_cols-1).zeros();
00190
00191 Col<eT> s_tmp = zeros< Col<eT> >(n_cols);
00192 s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
00193 s = s_tmp;
00194 }
00195
00196
00197 latent_out = s%s;
00198
00199 }
00200 else
00201 {
00202 if(n_rows == 1)
00203 {
00204 coeff_out = eye< Mat<eT> >(n_cols, n_cols);
00205
00206 score_out.copy_size(in);
00207 score_out.zeros();
00208
00209 latent_out.set_size(n_cols);
00210 latent_out.zeros();
00211 }
00212 else
00213 {
00214 coeff_out.reset();
00215 score_out.reset();
00216 latent_out.reset();
00217 }
00218 }
00219
00220 }
00221
00222
00223
00224
00225
00226
00227
00228
00229 template<typename eT>
00230 inline
00231 void
00232 op_princomp::direct_princomp
00233 (
00234 Mat<eT>& coeff_out,
00235 Mat<eT>& score_out,
00236 const Mat<eT>& in
00237 )
00238 {
00239 arma_extra_debug_sigprint();
00240
00241 const u32 n_rows = in.n_rows;
00242 const u32 n_cols = in.n_cols;
00243
00244 if(n_rows > 1)
00245 {
00246
00247 score_out = in - repmat(mean(in), n_rows, 1);
00248
00249
00250 Mat<eT> U;
00251 Col<eT> s;
00252
00253 const bool svd_ok = svd(U,s,coeff_out,score_out);
00254
00255 if(svd_ok == false)
00256 {
00257 arma_print("princomp(): singular value decomposition failed");
00258
00259 coeff_out.reset();
00260 score_out.reset();
00261
00262 return;
00263 }
00264
00265
00266
00267
00268 s /= std::sqrt(n_rows - 1);
00269
00270
00271 score_out *= coeff_out;
00272
00273 if(n_rows <= n_cols)
00274 {
00275 score_out.cols(n_rows-1,n_cols-1).zeros();
00276
00277 Col<eT> s_tmp = zeros< Col<eT> >(n_cols);
00278 s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
00279 s = s_tmp;
00280 }
00281 }
00282 else
00283 {
00284 if(n_rows == 1)
00285 {
00286 coeff_out = eye< Mat<eT> >(n_cols, n_cols);
00287 score_out.copy_size(in);
00288 score_out.zeros();
00289 }
00290 else
00291 {
00292 coeff_out.reset();
00293 score_out.reset();
00294 }
00295 }
00296 }
00297
00298
00299
00300
00301
00302
00303
00304 template<typename eT>
00305 inline
00306 void
00307 op_princomp::direct_princomp
00308 (
00309 Mat<eT>& coeff_out,
00310 const Mat<eT>& in
00311 )
00312 {
00313 arma_extra_debug_sigprint();
00314
00315 if(in.n_elem != 0)
00316 {
00317
00318 Mat<eT> U;
00319 Col<eT> s;
00320
00321 const Mat<eT> tmp = in - repmat(mean(in), in.n_rows, 1);
00322
00323 const bool svd_ok = svd(U,s,coeff_out, tmp);
00324
00325 if(svd_ok == false)
00326 {
00327 arma_print("princomp(): singular value decomposition failed");
00328
00329 coeff_out.reset();
00330 }
00331 }
00332 else
00333 {
00334 coeff_out.reset();
00335 }
00336 }
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347 template<typename T>
00348 inline
00349 void
00350 op_princomp::direct_princomp
00351 (
00352 Mat< std::complex<T> >& coeff_out,
00353 Mat< std::complex<T> >& score_out,
00354 Col<T>& latent_out,
00355 Col< std::complex<T> >& tsquared_out,
00356 const Mat< std::complex<T> >& in
00357 )
00358 {
00359 arma_extra_debug_sigprint();
00360
00361 typedef std::complex<T> eT;
00362
00363 const u32 n_rows = in.n_rows;
00364 const u32 n_cols = in.n_cols;
00365
00366 if(n_rows > 1)
00367 {
00368
00369 score_out = in - repmat(mean(in), n_rows, 1);
00370
00371
00372 Mat<eT> U;
00373 Col<T> s;
00374
00375 const bool svd_ok = svd(U,s,coeff_out,score_out);
00376
00377 if(svd_ok == false)
00378 {
00379 arma_print("princomp(): singular value decomposition failed");
00380
00381 coeff_out.reset();
00382 score_out.reset();
00383 latent_out.reset();
00384 tsquared_out.reset();
00385
00386 return;
00387 }
00388
00389
00390
00391
00392
00393 s /= std::sqrt(n_rows - 1);
00394
00395
00396 score_out *= coeff_out;
00397
00398 if(n_rows <= n_cols)
00399 {
00400 score_out.cols(n_rows-1,n_cols-1).zeros();
00401
00402 Col<T> s_tmp = zeros< Col<T> >(n_cols);
00403 s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
00404 s = s_tmp;
00405
00406
00407 s_tmp.rows(0,n_rows-2) = 1.0 / s_tmp.rows(0,n_rows-2);
00408 const Mat<eT> S = score_out * diagmat(Col<T>(s_tmp));
00409 tsquared_out = sum(S%S,1);
00410 }
00411 else
00412 {
00413
00414 const Mat<eT> S = score_out * diagmat(Col<T>(T(1) / s));
00415 tsquared_out = sum(S%S,1);
00416 }
00417
00418
00419 latent_out = s%s;
00420
00421 }
00422 else
00423 {
00424 if(n_rows == 1)
00425 {
00426 coeff_out = eye< Mat<eT> >(n_cols, n_cols);
00427
00428 score_out.copy_size(in);
00429 score_out.zeros();
00430
00431 latent_out.set_size(n_cols);
00432 latent_out.zeros();
00433
00434 tsquared_out.set_size(1);
00435 tsquared_out.zeros();
00436 }
00437 else
00438 {
00439 coeff_out.reset();
00440 score_out.reset();
00441 latent_out.reset();
00442 tsquared_out.reset();
00443 }
00444 }
00445 }
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455 template<typename T>
00456 inline
00457 void
00458 op_princomp::direct_princomp
00459 (
00460 Mat< std::complex<T> >& coeff_out,
00461 Mat< std::complex<T> >& score_out,
00462 Col<T>& latent_out,
00463 const Mat< std::complex<T> >& in
00464 )
00465 {
00466 arma_extra_debug_sigprint();
00467
00468 typedef std::complex<T> eT;
00469
00470 const u32 n_rows = in.n_rows;
00471 const u32 n_cols = in.n_cols;
00472
00473 if(n_rows > 1)
00474 {
00475
00476 score_out = in - repmat(mean(in), n_rows, 1);
00477
00478
00479 Mat<eT> U;
00480 Col< T> s;
00481
00482 const bool svd_ok = svd(U,s,coeff_out,score_out);
00483
00484 if(svd_ok == false)
00485 {
00486 arma_print("princomp(): singular value decomposition failed");
00487
00488 coeff_out.reset();
00489 score_out.reset();
00490 latent_out.reset();
00491
00492 return;
00493 }
00494
00495
00496
00497
00498
00499 s /= std::sqrt(n_rows - 1);
00500
00501
00502 score_out *= coeff_out;
00503
00504 if(n_rows <= n_cols)
00505 {
00506 score_out.cols(n_rows-1,n_cols-1).zeros();
00507
00508 Col<T> s_tmp = zeros< Col<T> >(n_cols);
00509 s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
00510 s = s_tmp;
00511 }
00512
00513
00514 latent_out = s%s;
00515
00516 }
00517 else
00518 {
00519 if(n_rows == 1)
00520 {
00521 coeff_out = eye< Mat<eT> >(n_cols, n_cols);
00522 score_out.copy_size(in);
00523 score_out.zeros();
00524 latent_out.set_size(n_cols);
00525 latent_out.zeros();
00526 }
00527 else
00528 {
00529 coeff_out.reset();
00530 score_out.reset();
00531 latent_out.reset();
00532 }
00533 }
00534 }
00535
00536
00537
00538
00539
00540
00541
00542
00543 template<typename T>
00544 inline
00545 void
00546 op_princomp::direct_princomp
00547 (
00548 Mat< std::complex<T> >& coeff_out,
00549 Mat< std::complex<T> >& score_out,
00550 const Mat< std::complex<T> >& in
00551 )
00552 {
00553 arma_extra_debug_sigprint();
00554
00555 typedef std::complex<T> eT;
00556
00557 const u32 n_rows = in.n_rows;
00558 const u32 n_cols = in.n_cols;
00559
00560 if(n_rows > 1)
00561 {
00562
00563 score_out = in - repmat(mean(in), n_rows, 1);
00564
00565
00566 Mat<eT> U;
00567 Col< T> s;
00568
00569 const bool svd_ok = svd(U,s,coeff_out,score_out);
00570
00571 if(svd_ok == false)
00572 {
00573 arma_print("princomp(): singular value decomposition failed");
00574
00575 coeff_out.reset();
00576 score_out.reset();
00577
00578 return;
00579 }
00580
00581
00582
00583
00584 s /= std::sqrt(n_rows - 1);
00585
00586
00587 score_out *= coeff_out;
00588
00589 if(n_rows <= n_cols)
00590 {
00591 score_out.cols(n_rows-1,n_cols-1).zeros();
00592 }
00593
00594 }
00595 else
00596 {
00597 if(n_rows == 1)
00598 {
00599 coeff_out = eye< Mat<eT> >(n_cols, n_cols);
00600
00601 score_out.copy_size(in);
00602 score_out.zeros();
00603 }
00604 else
00605 {
00606 coeff_out.reset();
00607 score_out.reset();
00608 }
00609 }
00610 }
00611
00612
00613
00614
00615
00616
00617
00618 template<typename T>
00619 inline
00620 void
00621 op_princomp::direct_princomp
00622 (
00623 Mat< std::complex<T> >& coeff_out,
00624 const Mat< std::complex<T> >& in
00625 )
00626 {
00627 arma_extra_debug_sigprint();
00628
00629 typedef typename std::complex<T> eT;
00630
00631 if(in.n_elem != 0)
00632 {
00633
00634 Mat<eT> U;
00635 Col< T> s;
00636
00637 const Mat<eT> tmp = in - repmat(mean(in), in.n_rows, 1);
00638
00639 const bool svd_ok = svd(U,s,coeff_out, tmp);
00640
00641 if(svd_ok == false)
00642 {
00643 arma_print("princomp(): singular value decomposition failed");
00644
00645 coeff_out.reset();
00646 }
00647 }
00648 else
00649 {
00650 coeff_out.reset();
00651 }
00652 }
00653
00654
00655
00656 template<typename T1>
00657 inline
00658 void
00659 op_princomp::apply
00660 (
00661 Mat<typename T1::elem_type>& out,
00662 const Op<T1,op_princomp>& in
00663 )
00664 {
00665 arma_extra_debug_sigprint();
00666
00667 typedef typename T1::elem_type eT;
00668
00669 const unwrap_check<T1> tmp(in.m, out);
00670 const Mat<eT>& A = tmp.M;
00671
00672 op_princomp::direct_princomp(out, A);
00673 }
00674
00675
00676
00677