KickJava   Java API By Example, From Geeks To Geeks.

Java > Open Source Codes > JSci > maths > KarhunenLoeve


1 package JSci.maths;
2
3 import JSci.maths.matrices.DoubleSquareMatrix;
4 import JSci.maths.vectors.AbstractDoubleVector;
5 import JSci.maths.vectors.DoubleVector;
6
7 /**
8 * This class implements the Karhunen-Loeve expansion.
9 * @author Daniel Lemire
10 */
11 public final class KarhunenLoeve {
12     double[][] data;
13
14     /**
15      * @param v an array where [i][j] is the jth component of the ith vector.
16     */
17     public KarhunenLoeve(double[][] v) {
18         setData(v);
19     }
20         public double[][] getProductMatrix() {
21                 return(getProductMatrix(data));
22         }
23         private static double[][] vectorToSquare (double[] v) {
24                 double[][] ans=new double[v.length][v.length];
25                 for(int k=0;k<v.length;k++) {
26                         for(int l=0;l<v.length;l++) {
27                                 ans[l][k]=v[k]*v[l];
28                         }
29                 }
30                 return(ans);
31         }
32         private static void add(double[][] a, double c, double[][] b) {
33                 for(int k=0;k<a.length;k++) {
34                         for(int l=0;l<a[k].length;l++) {
35                                 a[k][l]+=b[k][l]*c;
36                         }
37                 }
38         }
39     /**
40      * @param v an array where [i][j] is the jth component of the ith vector.
41      */
42     public static double[][] getProductMatrix (double[][] v) {
43         double[][] ans=new double[v[0].length][v[0].length];
44                 for(int k=0;k<v.length;k++) {
45                         add(ans,1.0/v.length,vectorToSquare(v[k]));
46                 }
47         return(ans);
48     }
49     public static double[][] getProductMatrix (double[] v) {
50         return(vectorToSquare(v));
51     }
52     /**
53     * Careful: doesn't generate a copy.
54     */
55     public double[][] getData() {
56         return(data);
57     }
58     /**
59     * Careful: doesn't generate a copy.
60     * @param v an array where [i][j] is the jth component of the ith vector.
61     */
62     public void setData(double[][] v) {
63         data=v;
64     }
65     /**
66     * Returns the eigenvectors ordered by the norm of the eigenvalues
67         * (from max to min).
68         * @exception MaximumIterationsExceededException if it can't compute
69         * the eigenvectors within the limited number of iterations allowed.
70     */
71     public AbstractDoubleVector[] getEigenvectors() throws MaximumIterationsExceededException {
72                 double[][] test=getProductMatrix(data);
73         DoubleSquareMatrix alpha=new DoubleSquareMatrix(test);
74                 DoubleVector[] beta=new DoubleVector[data[0].length];
75                 double[] eigen=LinearMath.eigenSolveSymmetric(alpha,beta);
76                 tri(eigen,beta);
77         return beta;
78     }
79     /**
80     * Returns the eigenvectors ordered by the norm of the eigenvalues
81         * (from max to min).
82     * @return an array where [i][j] is the jth component of the ith eigenvector.
83         * @exception MaximumIterationsExceededException if it can't compute
84         * the eigenvectors within the limited number of iterations allowed.
85     */
86     public double[][] getEigenvectorArrays() throws MaximumIterationsExceededException {
87                 AbstractDoubleVector[] beta=getEigenvectors();
88                 double[][] ans=new double[beta.length][beta[0].dimension()];
89                 for(int k=0;k<beta.length;k++) {
90                         for(int l=0;l<beta[k].dimension();l++) {
91                                 ans[k][l]=beta[k].getComponent(l);
92                         }
93                 }
94                 return(ans);
95     }
96     /** sort (french) */
97     private static void tri(double[] v, AbstractDoubleVector[] mat) {
98         double temp;
99         AbstractDoubleVector arraytemp;
100         boolean doitTrier=true;
101         while(doitTrier) {
102             doitTrier=false;
103             for(int k=0;k<v.length-1;k++) {
104                 if(v[k]<v[k+1]) {
105                     temp=v[k+1];
106                                         v[k+1]=v[k];
107                                         v[k]=temp;
108                                         doitTrier=true;
109                     arraytemp=mat[k+1];
110                     mat[k+1]=mat[k];
111                     mat[k]=arraytemp;
112                 }
113             }
114         }
115     }
116 }
117
Popular Tags