MatchingPursuit


1   package JSci.maths.wavelet;
2   
3   import JSci.maths.*;
4   /*************************************
5   * A Wavelet (and more) matching pursuit class
6   * Uses adaptative Morse coding for better
7   * performance. The MatchingPursuit is used
8   * to obtain a Time-Frequency representation (TF)
9   * through a fast algorithm.
10  * @author Daniel Lemire
11  *************************************/
12  public class MatchingPursuit extends BasisFunctionLibrary  implements Cloneable    {
13    private int[] Record;
14    private double[] Coefs;
15    private double[][][] TF=new double[0][][];
16  
17      public void add(MultiscaleFunction fprimary, MultiscaleFunction fdual) {
18      super.add(fprimary, fdual);
19      double[][][] TFback=TF;
20      TF=new double[getSize()][][];
21      System.arraycopy(TFback,0,TF,0,TFback.length);
22      TF[getSize()-1]=getTF(super.Fprimary[getSize()-1]);
23    }
24  
25    public Object   clone() {
26      MatchingPursuit c=(MatchingPursuit) super.clone();
27      if(this.Record!=null) {
28        c.Record=ArrayMath.copy(this.Record);
29      }
30      if(this.Coefs!=null) {
31        c.Coefs=ArrayMath.copy(this.Coefs);
32      }
33      if(this.TF!=null) {
34        TF=new double[this.TF.length][][];
35        System.arraycopy(TF,0,c.TF,0,TF.length);
36      }
37      return(c);
38    }
39  
40    /*****************************
41    * Allows one to trace back the
42    * matches
43    ******************************/
44    private void addToRecord (int K, double d) {
45      if(Record==null) {
46        Record=new int[1];
47        Record[0]=K;
48        Coefs=new double[1];
49        Coefs[0]=d;
50      } else {
51        int[] back=Record;
52        double[] backednorm=Coefs;
53        Record=new int[Record.length+1];
54        Coefs=new double[Record.length+1];
55        System.arraycopy(back,0,Record,0,Record.length-1);
56        System.arraycopy(backednorm,0,Coefs,0,Record.length-1);
57        //for(int k=0;k<Record.length-1;k++) {
58          //Record[k]=back[k];
59          //Coefs[k]=backednorm[k];
60       // }
61        Record[Record.length-1]=K;
62        Coefs[Record.length-1]=d;
63      }
64    }
65  
66    /******************************
67    * all matches are recorded
68    * so one can trace them back
69    *******************************/
70    public int[] getRecord() {
71      return(Record);
72    }
73  
74    /******************************
75    * Trace back how much of the
76    * norm was taken out at each
77    * match.
78    *******************************/
79    public double[] getRecordedNorms() {
80      double[] ans=new double[Record.length];
81      for(int k=0; k<Record.length;k++) {
82        ans[k]=Math.abs(Coefs[k])*Fprimary[k].norm();
83      }
84      return(ans);
85    }
86  
87    /******************************
88    * Recover the matching coefficients.
89    *******************************/
90    public double[] getCoefs() {
91      return(Coefs);
92    }
93  
94      /************************************
95      * Constructor
96      ************************************/
97      public MatchingPursuit (DiscreteFunction f) {
98          super(f);
99      }
100 
101   public void setData (DiscreteFunction f)   {
102     if(DFunction!=null) {
103       if(f.dimension()!=DFunction.dimension()) {
104         throw new IllegalArgumentException  ("You cannot change the dimension of the data object. Please create a new object.");
105       }
106     }
107         DFunction= (DiscreteFunction) f.clone();
108     Record=null;
109     Coefs=null;
110   }
111 
112   private MatchingPursuit() {
113   }
114 
115   /**************************************************
116   * Check the validity of the
117   * current matching algorithm.
118   * Will return an exception if the diagnostic
119   * fails. This makes sense: the software should
120   * stop if the algorithm isn't safely within
121   * the given tolerance.
122   * @exception IllegalArgumentException if
123   *   the matching fails
124   * @exception MaximumIterationsExceededException if it can't
125   *   match one of the elements of the dictionnary
126   * @exception IllegalArgumentException if tol is
127   *   is negative
128   ***************************************************/
129   public void diagnostic(double tol) throws MaximumIterationsExceededException {
130     if(tol<0) {
131       throw new IllegalArgumentException  ("A tolerance cannot be negative :"+tol+" < 0");
132     }
133     MatchingPursuit testmatch=(MatchingPursuit) this.clone();
134     double ener1, ener2;
135     for(int k=0;k<getSize();k++) {
136       testmatch.DFunction=(DiscreteFunction) super.Fprimary[k].clone();
137       ener1=testmatch.DFunction.norm();
138       testmatch.match();
139       ener2=testmatch.DFunction.norm();
140       if(ener2>ener1*tol) {
141         throw new IllegalArgumentException  ("Fail to match dictionnary element number "+k+" got a 'de facto' tolerance of "+(ener2/ener1));
142       }
143     }
144   }
145 
146     /************************************
147     * Convert the an element to a
148     * time-frequency representation
149     * using Fourier as a basis multiply
150     * everything by a coefficient b.
151     *************************************/
152     private double[][] getTF(int pos, double b) {
153     double[][] ans=ArrayMath.copy(TF[pos]);
154     for(int k=0;k<ans.length;k++) {
155       for(int l=0;l<ans[k].length;l++) {
156         ans[k][l]*=b;
157       }
158     }
159     return(ans);
160     }
161 
162     /************************************
163     * Convert the an element to a
164     * time-frequency representation
165     * using Fourier as a basis.
166     *************************************/
167     private double[][] getTF(DiscreteFunction a) {
168         int freqlength=(int) Math.ceil(DFunction.dimension()/2.0);
169         double freq;
170         double[][] TF=new double[freqlength][DFunction.dimension()];
171     double[] time=a.evaluate(0);
172         for(int k=0;k<freqlength;k++) {
173             freq=norm(DiscreteHilbertSpace.integrate(a,new Cosine(DFunction.dimension(),k)),DiscreteHilbertSpace.integrate(a,new Sine(DFunction.dimension(),k)));
174             for(int l=0;l<DFunction.dimension();l++) {
175                 TF[k][l]=time[l]*freq;
176             }
177         }
178         return(TF);
179     }
180 
181 
182     /************************************
183     * Does the matching 1 time and return
184     * the TF representation. The TF representation
185   * may contain negative values and so,
186   * taking the absolute value of the
187   * result will often be useful.
188     *************************************/
189     public double[][] match() {
190         double coef=getWeigth(0);
191     int pos=0;
192     double min=getResidue(0);
193         double resitemp,weigthtemp;
194     double[] soltemp;
195         for(int k=1;k<getSize();k++) {
196       weigthtemp=DiscreteHilbertSpace.integrate(DFunction,Fdual[k]);
197       soltemp=DiscreteHilbertSpace.add(DFunction,-weigthtemp,Fprimary[k]);
198             resitemp=ArrayMath.norm(soltemp);
199             if(resitemp<min) {
200                 min=resitemp;
201                 coef=weigthtemp;
202                 pos=k;
203             }
204         }
205     addToRecord(pos,coef);
206     DFunction=new DiscreteFunction(DiscreteHilbertSpace.add(DFunction,-coef,Fprimary[pos]));
207         return(getTF(pos,coef));
208     }
209 
210   /******************************************
211   * Repeatly match until it remains less
212   * than tol * 100 percent of the original
213   * L2 norm; no matter what, at least one
214   * match will be done. It will return
215   * the TF representaiton. The TF representation
216   * may contain negative values and so,
217   * taking the absolute value of the
218   * result will often be useful.
219   * @param tol percentile of energy
220   * @exception IllegalArgumentException
221   *   if tol is not within the interval [0,1]
222   * @exception MaximumIterationsExceededException
223   *   if the number of required match exceeds 5 times
224   *   the size of the dictionnary (it should be a more
225   *   than confortable margin unless the problem is
226   *   ill-posed, change the dictionnary if it doesn't work)
227   *******************************************/
228   public double[][] matchAll (double tol) throws MaximumIterationsExceededException {
229 
230     if((tol>1)||(tol<0)) {
231       throw new IllegalArgumentException  ("The percentile should be between 0 and 1: "+tol);
232     }
233     int TimesAround=0;
234     double ener0=DFunction.norm();
235     double ts=ener0*tol;
236     int maxTimesAround=5*getSize();
237     double[][] ans=new double[0][0];
238     while(DFunction.norm()>ts) {
239       TimesAround++;
240       if(TimesAround>maxTimesAround) {
241         throw new MaximumIterationsExceededException("Impossible to match to the desired precision ("+tol+") with this dictionnary of size "+getSize()+" after "+TimesAround+" iterations. You might want to expand the dictionnary.");
242       }
243       if(ans.length!=0) {
244         ans=add(ans,match());
245       } else {
246         ans=match();
247       }
248     }
249     return(ans);
250   }
251 
252 
253 
254   /******************************************
255   * Force the system to select the given
256   * element as the best match. It will return
257   * the TF representaiton. The TF representation
258   * may contain negative values and so,
259   * taking the absolute value of the
260   * result will often be useful.
261   *******************************************/
262   public double[][] forcedMatch(int pos) {
263     double coef=DiscreteHilbertSpace.integrate(DFunction,Fdual[pos]);
264     addToRecord(pos,coef);
265     DFunction=new DiscreteFunction(DiscreteHilbertSpace.add(DFunction,-coef,Fprimary[pos]));
266         return(getTF(pos,coef));
267   }
268 
269     /*********************************
270     * Does the matching j times and
271     * return the TF representation.
272   * The TF representation
273   * may contain negative values and so,
274   * taking the absolute value of the
275   * result will often be useful.
276   * @param j number of iterations
277   * @exception IllegalArgumentException if j is not positive
278     **********************************/
279     public double[][] match(int j) {
280     if(j<=0) {
281       throw new IllegalArgumentException  ("Can't do "+j+" matching... Must be a positive number.");
282     }
283         double[][] ans=match();
284         for(int k=1;k<j;k++) {
285             ans=add(ans,match());
286         }
287         return(ans);
288     }
289 
290 }
291 
292
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags