pca.c 2.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102
  1. /*
  2. * principal component analysis (PCA)
  3. * Copyright (c) 2004 Michael Niedermayer <michaelni@gmx.at>
  4. *
  5. * This file is part of FFmpeg.
  6. *
  7. * FFmpeg is free software; you can redistribute it and/or
  8. * modify it under the terms of the GNU Lesser General Public
  9. * License as published by the Free Software Foundation; either
  10. * version 2.1 of the License, or (at your option) any later version.
  11. *
  12. * FFmpeg is distributed in the hope that it will be useful,
  13. * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  15. * Lesser General Public License for more details.
  16. *
  17. * You should have received a copy of the GNU Lesser General Public
  18. * License along with FFmpeg; if not, write to the Free Software
  19. * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  20. */
  21. #include "libavutil/pca.c"
  22. #include "libavutil/lfg.h"
  23. #undef printf
  24. #include <stdio.h>
  25. #include <stdlib.h>
  26. int main(void){
  27. PCA *pca;
  28. int i, j, k;
  29. #define LEN 8
  30. double eigenvector[LEN*LEN];
  31. double eigenvalue[LEN];
  32. AVLFG prng;
  33. av_lfg_init(&prng, 1);
  34. pca= ff_pca_init(LEN);
  35. for(i=0; i<9000000; i++){
  36. double v[2*LEN+100];
  37. double sum=0;
  38. int pos = av_lfg_get(&prng) % LEN;
  39. int v2 = av_lfg_get(&prng) % 101 - 50;
  40. v[0] = av_lfg_get(&prng) % 101 - 50;
  41. for(j=1; j<8; j++){
  42. if(j<=pos) v[j]= v[0];
  43. else v[j]= v2;
  44. sum += v[j];
  45. }
  46. /* for(j=0; j<LEN; j++){
  47. v[j] -= v[pos];
  48. }*/
  49. // sum += av_lfg_get(&prng) % 10;
  50. /* for(j=0; j<LEN; j++){
  51. v[j] -= sum/LEN;
  52. }*/
  53. // lbt1(v+100,v+100,LEN);
  54. ff_pca_add(pca, v);
  55. }
  56. ff_pca(pca, eigenvector, eigenvalue);
  57. for(i=0; i<LEN; i++){
  58. pca->count= 1;
  59. pca->mean[i]= 0;
  60. // (0.5^|x|)^2 = 0.5^2|x| = 0.25^|x|
  61. // pca.covariance[i + i*LEN]= pow(0.5, fabs
  62. for(j=i; j<LEN; j++){
  63. printf("%f ", pca->covariance[i + j*LEN]);
  64. }
  65. printf("\n");
  66. }
  67. for(i=0; i<LEN; i++){
  68. double v[LEN];
  69. double error=0;
  70. memset(v, 0, sizeof(v));
  71. for(j=0; j<LEN; j++){
  72. for(k=0; k<LEN; k++){
  73. v[j] += pca->covariance[FFMIN(k,j) + FFMAX(k,j)*LEN] * eigenvector[i + k*LEN];
  74. }
  75. v[j] /= eigenvalue[i];
  76. error += fabs(v[j] - eigenvector[i + j*LEN]);
  77. }
  78. printf("%f ", error);
  79. }
  80. printf("\n");
  81. for(i=0; i<LEN; i++){
  82. for(j=0; j<LEN; j++){
  83. printf("%9.6f ", eigenvector[i + j*LEN]);
  84. }
  85. printf(" %9.1f %f\n", eigenvalue[i], eigenvalue[i]/eigenvalue[0]);
  86. }
  87. return 0;
  88. }