Mercurial > hg > segmenter-vamp-plugin

annotate armadillo-2.4.4/include/armadillo_bits/op_princomp_meat.hpp @ 5:79b343f3e4b8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
In thi version the problem of letters assigned to each segment has been solved.
author maxzanoni76 <max.zanoni@eecs.qmul.ac.uk>
date Wed, 11 Apr 2012 13:48:13 +0100
parents 8b6102e2a9b0
children
rev   line source
max@0 1 // Copyright (C) 2010-2011 NICTA (www.nicta.com.au)
max@0 2 // Copyright (C) 2010-2011 Conrad Sanderson
max@0 3 // Copyright (C) 2010 Dimitrios Bouzas
max@0 4 // Copyright (C) 2011 Stanislav Funiak
max@0 5 //
max@0 6 // This file is part of the Armadillo C++ library.
max@0 7 // It is provided without any warranty of fitness
max@0 8 // for any purpose. You can redistribute this file
max@0 9 // and/or modify it under the terms of the GNU
max@0 10 // Lesser General Public License (LGPL) as published
max@0 11 // by the Free Software Foundation, either version 3
max@0 12 // of the License or (at your option) any later version.
max@0 13 // (see http://www.opensource.org/licenses for more info)
max@0 14
max@0 15
max@0 16 //! \addtogroup op_princomp
max@0 17 //! @{
max@0 18
max@0 19
max@0 20
max@0 21 //! \brief
max@0 22 //! principal component analysis -- 4 arguments version
max@0 23 //! computation is done via singular value decomposition
max@0 24 //! coeff_out -> principal component coefficients
max@0 25 //! score_out -> projected samples
max@0 26 //! latent_out -> eigenvalues of principal vectors
max@0 27 //! tsquared_out -> Hotelling's T^2 statistic
max@0 28 template<typename eT>
max@0 29 inline
max@0 30 bool
max@0 31 op_princomp::direct_princomp
max@0 32 (
max@0 33 Mat<eT>& coeff_out,
max@0 34 Mat<eT>& score_out,
max@0 35 Col<eT>& latent_out,
max@0 36 Col<eT>& tsquared_out,
max@0 37 const Mat<eT>& in
max@0 38 )
max@0 39 {
max@0 40 arma_extra_debug_sigprint();
max@0 41
max@0 42 const uword n_rows = in.n_rows;
max@0 43 const uword n_cols = in.n_cols;
max@0 44
max@0 45 if(n_rows > 1) // more than one sample
max@0 46 {
max@0 47 // subtract the mean - use score_out as temporary matrix
max@0 48 score_out = in - repmat(mean(in), n_rows, 1);
max@0 49
max@0 50 // singular value decomposition
max@0 51 Mat<eT> U;
max@0 52 Col<eT> s;
max@0 53
max@0 54 const bool svd_ok = svd(U,s,coeff_out,score_out);
max@0 55
max@0 56 if(svd_ok == false)
max@0 57 {
max@0 58 return false;
max@0 59 }
max@0 60
max@0 61
max@0 62 //U.reset(); // TODO: do we need this ? U will get automatically deleted anyway
max@0 63
max@0 64 // normalize the eigenvalues
max@0 65 s /= std::sqrt( double(n_rows - 1) );
max@0 66
max@0 67 // project the samples to the principals
max@0 68 score_out *= coeff_out;
max@0 69
max@0 70 if(n_rows <= n_cols) // number of samples is less than their dimensionality
max@0 71 {
max@0 72 score_out.cols(n_rows-1,n_cols-1).zeros();
max@0 73
max@0 74 //Col<eT> s_tmp = zeros< Col<eT> >(n_cols);
max@0 75 Col<eT> s_tmp(n_cols);
max@0 76 s_tmp.zeros();
max@0 77
max@0 78 s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
max@0 79 s = s_tmp;
max@0 80
max@0 81 // compute the Hotelling's T-squared
max@0 82 s_tmp.rows(0,n_rows-2) = eT(1) / s_tmp.rows(0,n_rows-2);
max@0 83
max@0 84 const Mat<eT> S = score_out * diagmat(Col<eT>(s_tmp));
max@0 85 tsquared_out = sum(S%S,1);
max@0 86 }
max@0 87 else
max@0 88 {
max@0 89 // compute the Hotelling's T-squared
max@0 90 const Mat<eT> S = score_out * diagmat(Col<eT>( eT(1) / s));
max@0 91 tsquared_out = sum(S%S,1);
max@0 92 }
max@0 93
max@0 94 // compute the eigenvalues of the principal vectors
max@0 95 latent_out = s%s;
max@0 96 }
max@0 97 else // 0 or 1 samples
max@0 98 {
max@0 99 coeff_out.eye(n_cols, n_cols);
max@0 100
max@0 101 score_out.copy_size(in);
max@0 102 score_out.zeros();
max@0 103
max@0 104 latent_out.set_size(n_cols);
max@0 105 latent_out.zeros();
max@0 106
max@0 107 tsquared_out.set_size(n_rows);
max@0 108 tsquared_out.zeros();
max@0 109 }
max@0 110
max@0 111 return true;
max@0 112 }
max@0 113
max@0 114
max@0 115
max@0 116 //! \brief
max@0 117 //! principal component analysis -- 3 arguments version
max@0 118 //! computation is done via singular value decomposition
max@0 119 //! coeff_out -> principal component coefficients
max@0 120 //! score_out -> projected samples
max@0 121 //! latent_out -> eigenvalues of principal vectors
max@0 122 template<typename eT>
max@0 123 inline
max@0 124 bool
max@0 125 op_princomp::direct_princomp
max@0 126 (
max@0 127 Mat<eT>& coeff_out,
max@0 128 Mat<eT>& score_out,
max@0 129 Col<eT>& latent_out,
max@0 130 const Mat<eT>& in
max@0 131 )
max@0 132 {
max@0 133 arma_extra_debug_sigprint();
max@0 134
max@0 135 const uword n_rows = in.n_rows;
max@0 136 const uword n_cols = in.n_cols;
max@0 137
max@0 138 if(n_rows > 1) // more than one sample
max@0 139 {
max@0 140 // subtract the mean - use score_out as temporary matrix
max@0 141 score_out = in - repmat(mean(in), n_rows, 1);
max@0 142
max@0 143 // singular value decomposition
max@0 144 Mat<eT> U;
max@0 145 Col<eT> s;
max@0 146
max@0 147 const bool svd_ok = svd(U,s,coeff_out,score_out);
max@0 148
max@0 149 if(svd_ok == false)
max@0 150 {
max@0 151 return false;
max@0 152 }
max@0 153
max@0 154
max@0 155 // U.reset();
max@0 156
max@0 157 // normalize the eigenvalues
max@0 158 s /= std::sqrt( double(n_rows - 1) );
max@0 159
max@0 160 // project the samples to the principals
max@0 161 score_out *= coeff_out;
max@0 162
max@0 163 if(n_rows <= n_cols) // number of samples is less than their dimensionality
max@0 164 {
max@0 165 score_out.cols(n_rows-1,n_cols-1).zeros();
max@0 166
max@0 167 Col<eT> s_tmp = zeros< Col<eT> >(n_cols);
max@0 168 s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
max@0 169 s = s_tmp;
max@0 170 }
max@0 171
max@0 172 // compute the eigenvalues of the principal vectors
max@0 173 latent_out = s%s;
max@0 174
max@0 175 }
max@0 176 else // 0 or 1 samples
max@0 177 {
max@0 178 coeff_out.eye(n_cols, n_cols);
max@0 179
max@0 180 score_out.copy_size(in);
max@0 181 score_out.zeros();
max@0 182
max@0 183 latent_out.set_size(n_cols);
max@0 184 latent_out.zeros();
max@0 185 }
max@0 186
max@0 187 return true;
max@0 188 }
max@0 189
max@0 190
max@0 191
max@0 192 //! \brief
max@0 193 //! principal component analysis -- 2 arguments version
max@0 194 //! computation is done via singular value decomposition
max@0 195 //! coeff_out -> principal component coefficients
max@0 196 //! score_out -> projected samples
max@0 197 template<typename eT>
max@0 198 inline
max@0 199 bool
max@0 200 op_princomp::direct_princomp
max@0 201 (
max@0 202 Mat<eT>& coeff_out,
max@0 203 Mat<eT>& score_out,
max@0 204 const Mat<eT>& in
max@0 205 )
max@0 206 {
max@0 207 arma_extra_debug_sigprint();
max@0 208
max@0 209 const uword n_rows = in.n_rows;
max@0 210 const uword n_cols = in.n_cols;
max@0 211
max@0 212 if(n_rows > 1) // more than one sample
max@0 213 {
max@0 214 // subtract the mean - use score_out as temporary matrix
max@0 215 score_out = in - repmat(mean(in), n_rows, 1);
max@0 216
max@0 217 // singular value decomposition
max@0 218 Mat<eT> U;
max@0 219 Col<eT> s;
max@0 220
max@0 221 const bool svd_ok = svd(U,s,coeff_out,score_out);
max@0 222
max@0 223 if(svd_ok == false)
max@0 224 {
max@0 225 return false;
max@0 226 }
max@0 227
max@0 228 // U.reset();
max@0 229
max@0 230 // normalize the eigenvalues
max@0 231 s /= std::sqrt( double(n_rows - 1) );
max@0 232
max@0 233 // project the samples to the principals
max@0 234 score_out *= coeff_out;
max@0 235
max@0 236 if(n_rows <= n_cols) // number of samples is less than their dimensionality
max@0 237 {
max@0 238 score_out.cols(n_rows-1,n_cols-1).zeros();
max@0 239
max@0 240 Col<eT> s_tmp = zeros< Col<eT> >(n_cols);
max@0 241 s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
max@0 242 s = s_tmp;
max@0 243 }
max@0 244 }
max@0 245 else // 0 or 1 samples
max@0 246 {
max@0 247 coeff_out.eye(n_cols, n_cols);
max@0 248 score_out.copy_size(in);
max@0 249 score_out.zeros();
max@0 250 }
max@0 251
max@0 252 return true;
max@0 253 }
max@0 254
max@0 255
max@0 256
max@0 257 //! \brief
max@0 258 //! principal component analysis -- 1 argument version
max@0 259 //! computation is done via singular value decomposition
max@0 260 //! coeff_out -> principal component coefficients
max@0 261 template<typename eT>
max@0 262 inline
max@0 263 bool
max@0 264 op_princomp::direct_princomp
max@0 265 (
max@0 266 Mat<eT>& coeff_out,
max@0 267 const Mat<eT>& in
max@0 268 )
max@0 269 {
max@0 270 arma_extra_debug_sigprint();
max@0 271
max@0 272 if(in.n_elem != 0)
max@0 273 {
max@0 274 // singular value decomposition
max@0 275 Mat<eT> U;
max@0 276 Col<eT> s;
max@0 277
max@0 278 const Mat<eT> tmp = in - repmat(mean(in), in.n_rows, 1);
max@0 279
max@0 280 const bool svd_ok = svd(U,s,coeff_out, tmp);
max@0 281
max@0 282 if(svd_ok == false)
max@0 283 {
max@0 284 return false;
max@0 285 }
max@0 286 }
max@0 287 else
max@0 288 {
max@0 289 coeff_out.eye(in.n_cols, in.n_cols);
max@0 290 }
max@0 291
max@0 292 return true;
max@0 293 }
max@0 294
max@0 295
max@0 296
max@0 297 //! \brief
max@0 298 //! principal component analysis -- 4 arguments complex version
max@0 299 //! computation is done via singular value decomposition
max@0 300 //! coeff_out -> principal component coefficients
max@0 301 //! score_out -> projected samples
max@0 302 //! latent_out -> eigenvalues of principal vectors
max@0 303 //! tsquared_out -> Hotelling's T^2 statistic
max@0 304 template<typename T>
max@0 305 inline
max@0 306 bool
max@0 307 op_princomp::direct_princomp
max@0 308 (
max@0 309 Mat< std::complex<T> >& coeff_out,
max@0 310 Mat< std::complex<T> >& score_out,
max@0 311 Col<T>& latent_out,
max@0 312 Col< std::complex<T> >& tsquared_out,
max@0 313 const Mat< std::complex<T> >& in
max@0 314 )
max@0 315 {
max@0 316 arma_extra_debug_sigprint();
max@0 317
max@0 318 typedef std::complex<T> eT;
max@0 319
max@0 320 const uword n_rows = in.n_rows;
max@0 321 const uword n_cols = in.n_cols;
max@0 322
max@0 323 if(n_rows > 1) // more than one sample
max@0 324 {
max@0 325 // subtract the mean - use score_out as temporary matrix
max@0 326 score_out = in - repmat(mean(in), n_rows, 1);
max@0 327
max@0 328 // singular value decomposition
max@0 329 Mat<eT> U;
max@0 330 Col<T> s;
max@0 331
max@0 332 const bool svd_ok = svd(U,s,coeff_out,score_out);
max@0 333
max@0 334 if(svd_ok == false)
max@0 335 {
max@0 336 return false;
max@0 337 }
max@0 338
max@0 339
max@0 340 //U.reset();
max@0 341
max@0 342 // normalize the eigenvalues
max@0 343 s /= std::sqrt( double(n_rows - 1) );
max@0 344
max@0 345 // project the samples to the principals
max@0 346 score_out *= coeff_out;
max@0 347
max@0 348 if(n_rows <= n_cols) // number of samples is less than their dimensionality
max@0 349 {
max@0 350 score_out.cols(n_rows-1,n_cols-1).zeros();
max@0 351
max@0 352 Col<T> s_tmp = zeros< Col<T> >(n_cols);
max@0 353 s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
max@0 354 s = s_tmp;
max@0 355
max@0 356 // compute the Hotelling's T-squared
max@0 357 s_tmp.rows(0,n_rows-2) = 1.0 / s_tmp.rows(0,n_rows-2);
max@0 358 const Mat<eT> S = score_out * diagmat(Col<T>(s_tmp));
max@0 359 tsquared_out = sum(S%S,1);
max@0 360 }
max@0 361 else
max@0 362 {
max@0 363 // compute the Hotelling's T-squared
max@0 364 const Mat<eT> S = score_out * diagmat(Col<T>(T(1) / s));
max@0 365 tsquared_out = sum(S%S,1);
max@0 366 }
max@0 367
max@0 368 // compute the eigenvalues of the principal vectors
max@0 369 latent_out = s%s;
max@0 370
max@0 371 }
max@0 372 else // 0 or 1 samples
max@0 373 {
max@0 374 coeff_out.eye(n_cols, n_cols);
max@0 375
max@0 376 score_out.copy_size(in);
max@0 377 score_out.zeros();
max@0 378
max@0 379 latent_out.set_size(n_cols);
max@0 380 latent_out.zeros();
max@0 381
max@0 382 tsquared_out.set_size(n_rows);
max@0 383 tsquared_out.zeros();
max@0 384 }
max@0 385
max@0 386 return true;
max@0 387 }
max@0 388
max@0 389
max@0 390
max@0 391 //! \brief
max@0 392 //! principal component analysis -- 3 arguments complex version
max@0 393 //! computation is done via singular value decomposition
max@0 394 //! coeff_out -> principal component coefficients
max@0 395 //! score_out -> projected samples
max@0 396 //! latent_out -> eigenvalues of principal vectors
max@0 397 template<typename T>
max@0 398 inline
max@0 399 bool
max@0 400 op_princomp::direct_princomp
max@0 401 (
max@0 402 Mat< std::complex<T> >& coeff_out,
max@0 403 Mat< std::complex<T> >& score_out,
max@0 404 Col<T>& latent_out,
max@0 405 const Mat< std::complex<T> >& in
max@0 406 )
max@0 407 {
max@0 408 arma_extra_debug_sigprint();
max@0 409
max@0 410 typedef std::complex<T> eT;
max@0 411
max@0 412 const uword n_rows = in.n_rows;
max@0 413 const uword n_cols = in.n_cols;
max@0 414
max@0 415 if(n_rows > 1) // more than one sample
max@0 416 {
max@0 417 // subtract the mean - use score_out as temporary matrix
max@0 418 score_out = in - repmat(mean(in), n_rows, 1);
max@0 419
max@0 420 // singular value decomposition
max@0 421 Mat<eT> U;
max@0 422 Col< T> s;
max@0 423
max@0 424 const bool svd_ok = svd(U,s,coeff_out,score_out);
max@0 425
max@0 426 if(svd_ok == false)
max@0 427 {
max@0 428 return false;
max@0 429 }
max@0 430
max@0 431
max@0 432 // U.reset();
max@0 433
max@0 434 // normalize the eigenvalues
max@0 435 s /= std::sqrt( double(n_rows - 1) );
max@0 436
max@0 437 // project the samples to the principals
max@0 438 score_out *= coeff_out;
max@0 439
max@0 440 if(n_rows <= n_cols) // number of samples is less than their dimensionality
max@0 441 {
max@0 442 score_out.cols(n_rows-1,n_cols-1).zeros();
max@0 443
max@0 444 Col<T> s_tmp = zeros< Col<T> >(n_cols);
max@0 445 s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
max@0 446 s = s_tmp;
max@0 447 }
max@0 448
max@0 449 // compute the eigenvalues of the principal vectors
max@0 450 latent_out = s%s;
max@0 451
max@0 452 }
max@0 453 else // 0 or 1 samples
max@0 454 {
max@0 455 coeff_out.eye(n_cols, n_cols);
max@0 456
max@0 457 score_out.copy_size(in);
max@0 458 score_out.zeros();
max@0 459
max@0 460 latent_out.set_size(n_cols);
max@0 461 latent_out.zeros();
max@0 462 }
max@0 463
max@0 464 return true;
max@0 465 }
max@0 466
max@0 467
max@0 468
max@0 469 //! \brief
max@0 470 //! principal component analysis -- 2 arguments complex version
max@0 471 //! computation is done via singular value decomposition
max@0 472 //! coeff_out -> principal component coefficients
max@0 473 //! score_out -> projected samples
max@0 474 template<typename T>
max@0 475 inline
max@0 476 bool
max@0 477 op_princomp::direct_princomp
max@0 478 (
max@0 479 Mat< std::complex<T> >& coeff_out,
max@0 480 Mat< std::complex<T> >& score_out,
max@0 481 const Mat< std::complex<T> >& in
max@0 482 )
max@0 483 {
max@0 484 arma_extra_debug_sigprint();
max@0 485
max@0 486 typedef std::complex<T> eT;
max@0 487
max@0 488 const uword n_rows = in.n_rows;
max@0 489 const uword n_cols = in.n_cols;
max@0 490
max@0 491 if(n_rows > 1) // more than one sample
max@0 492 {
max@0 493 // subtract the mean - use score_out as temporary matrix
max@0 494 score_out = in - repmat(mean(in), n_rows, 1);
max@0 495
max@0 496 // singular value decomposition
max@0 497 Mat<eT> U;
max@0 498 Col< T> s;
max@0 499
max@0 500 const bool svd_ok = svd(U,s,coeff_out,score_out);
max@0 501
max@0 502 if(svd_ok == false)
max@0 503 {
max@0 504 return false;
max@0 505 }
max@0 506
max@0 507 // U.reset();
max@0 508
max@0 509 // normalize the eigenvalues
max@0 510 s /= std::sqrt( double(n_rows - 1) );
max@0 511
max@0 512 // project the samples to the principals
max@0 513 score_out *= coeff_out;
max@0 514
max@0 515 if(n_rows <= n_cols) // number of samples is less than their dimensionality
max@0 516 {
max@0 517 score_out.cols(n_rows-1,n_cols-1).zeros();
max@0 518 }
max@0 519
max@0 520 }
max@0 521 else // 0 or 1 samples
max@0 522 {
max@0 523 coeff_out.eye(n_cols, n_cols);
max@0 524
max@0 525 score_out.copy_size(in);
max@0 526 score_out.zeros();
max@0 527 }
max@0 528
max@0 529 return true;
max@0 530 }
max@0 531
max@0 532
max@0 533
max@0 534 //! \brief
max@0 535 //! principal component analysis -- 1 argument complex version
max@0 536 //! computation is done via singular value decomposition
max@0 537 //! coeff_out -> principal component coefficients
max@0 538 template<typename T>
max@0 539 inline
max@0 540 bool
max@0 541 op_princomp::direct_princomp
max@0 542 (
max@0 543 Mat< std::complex<T> >& coeff_out,
max@0 544 const Mat< std::complex<T> >& in
max@0 545 )
max@0 546 {
max@0 547 arma_extra_debug_sigprint();
max@0 548
max@0 549 typedef typename std::complex<T> eT;
max@0 550
max@0 551 if(in.n_elem != 0)
max@0 552 {
max@0 553 // singular value decomposition
max@0 554 Mat<eT> U;
max@0 555 Col< T> s;
max@0 556
max@0 557 const Mat<eT> tmp = in - repmat(mean(in), in.n_rows, 1);
max@0 558
max@0 559 const bool svd_ok = svd(U,s,coeff_out, tmp);
max@0 560
max@0 561 if(svd_ok == false)
max@0 562 {
max@0 563 return false;
max@0 564 }
max@0 565 }
max@0 566 else
max@0 567 {
max@0 568 coeff_out.eye(in.n_cols, in.n_cols);
max@0 569 }
max@0 570
max@0 571 return true;
max@0 572 }
max@0 573
max@0 574
max@0 575
max@0 576 template<typename T1>
max@0 577 inline
max@0 578 void
max@0 579 op_princomp::apply
max@0 580 (
max@0 581 Mat<typename T1::elem_type>& out,
max@0 582 const Op<T1,op_princomp>& in
max@0 583 )
max@0 584 {
max@0 585 arma_extra_debug_sigprint();
max@0 586
max@0 587 typedef typename T1::elem_type eT;
max@0 588
max@0 589 const unwrap_check<T1> tmp(in.m, out);
max@0 590 const Mat<eT>& A = tmp.M;
max@0 591
max@0 592 const bool status = op_princomp::direct_princomp(out, A);
max@0 593
max@0 594 if(status == false)
max@0 595 {
max@0 596 out.reset();
max@0 597
max@0 598 arma_bad("princomp(): failed to converge");
max@0 599 }
max@0 600 }
max@0 601
max@0 602
max@0 603
max@0 604 //! @}