001    package util.wavelet;
002    
003    import numal.lowprecision.AnalyticProblems;
004    import numal.lowprecision.Function;
005    
006    import java.text.DecimalFormat;
007    import static java.lang.Math.*;
008    
009    /** 
010     * A function representing a gaussian to be used in {@link AnalyticProblems#marquardt}.
011     *
012     * @author John Talbot
013     */
014    public class GaussianFit implements Function {
015    
016      /** Ordinate array corresponding to absica y (zero position not significant) */
017      float[] x; 
018    
019      /** Array of values corresponding to x which are to be fitted with f(x, par) (zero position not significant) */
020      float[] y;
021    
022      public GaussianFit(float[] x, float[] y) {
023        this.x = x;
024        this.y = y;
025      }
026    
027      static final int HEIGHT = 1;
028      static final int CENTER = 2;
029      static final int SIGMA  = 3;
030    
031      static final float C0 = 3.5801f;
032      static final float C1 = 0.0001f;
033      
034      /** 
035       * Get the wavelength as a function of SDSS spectra pixel. The function is 10^(C0 + pixel*C1).
036       * @param pixel an integer representing the position in the SDSS spectra
037       * @return wavelength
038       */
039      public static float getX(int pixel) {
040          return (float) pow(10d, (double) (C0 + pixel * C1));
041      }
042    
043      public static void main(String args[]) {
044    
045        int m = 3848;
046        float[] x = new float[m + 1];
047        float[] y = new float[m + 1];
048     
049        float lambda = (float)(getX(1) + (getX(m) - getX(1)) / 2d);
050        float[] originalParameters = {0f, 110f, lambda + 20, 140f};  // height, center, sigma
051    
052        for (int i = 1; i <= m; i++) {
053            x[i] = getX(i);
054            y[i] = f(x[i], originalParameters) + WaveletTransform.getGaussianDeviate();
055        }
056        Function function = new GaussianFit(x, y);
057    
058        // these parameters are guesses
059        float[] parameter = {0f, 100f, lambda, 150f};  // height, center, sigma
060        float rv[] = new float[x.length];  // residual vector f(x[i],par) - y[i]
061        final float jjinv[][] = new float[parameter.length][parameter.length];
062        float[] out = new float[8];
063    
064        AnalyticProblems.marquardt(x.length - 1, parameter.length - 1, parameter, rv, jjinv, function, out);
065    
066    /*
067        float[] yfit = function.getAnalyticValues(m, parameter);
068        float[] ya= new float[m + 1];  
069        for (int i = 1; i <= m; i++) {
070            ya[i] = f(x[i], originalParameters);
071        }
072        for (int i = 1; i <= m; i++) System.out.println(x[i] + " " + y[i] + " " + ya[i] + " " + yfit[i]);
073    */
074    
075        DecimalFormat fiveDigit = new DecimalFormat("00.00000E0");
076        DecimalFormat oneDigit = new DecimalFormat("00.0");
077        System.out.println("===============Gaussian fit================");
078        System.out.println("Parameters:\n  " +
079          "Height= " + fiveDigit.format(parameter[HEIGHT]) + "   " +
080          "(" + fiveDigit.format(originalParameters[HEIGHT]) + ") " +
081          "Center= " + fiveDigit.format(parameter[CENTER]) + "   " +
082          "(" + fiveDigit.format(originalParameters[CENTER]) + ") " +
083          "Sigma= "  + fiveDigit.format(parameter[SIGMA]) + 
084          "(" + fiveDigit.format(originalParameters[SIGMA]) + ") " +
085          "\n\nOUT:\n    " +
086          fiveDigit.format(out[7]) + "\n    " +
087          fiveDigit.format(out[2]) + "\n    " +
088          fiveDigit.format(out[6]) + "\n    " +
089          fiveDigit.format(out[3]) + "\n    " +
090          fiveDigit.format(out[4]) + "\n    " +
091          fiveDigit.format(out[5]) + "\n    " +
092          fiveDigit.format(out[1]) + "\n\nLast residual vector:\n  " +
093          oneDigit.format(rv[1]) + "  " + 
094          oneDigit.format(rv[2]) + "  " + 
095          oneDigit.format(rv[3]) + "  " + 
096          oneDigit.format(rv[4]) + "  " + 
097          oneDigit.format(rv[5]) + "  " + 
098          oneDigit.format(rv[6]));
099        }
100    
101      
102      /** 
103       * Get the value of this function for a given abcissa and parameters.
104       *
105       * @param    x          abcissa at which this function is to be evaluated
106       * @param    parameter  array of function parameters
107       * @return   the value of this function for the given x and parameters
108       */
109      public static float f(float x, float parameter[]) {
110          float x0 = x - parameter[CENTER];
111          return parameter[HEIGHT] * (float) exp(-x0 * x0 / (2f * parameter[SIGMA] * parameter[SIGMA]));    
112      }
113    
114      // implementation of the Function interface :
115      public float[] getAnalyticValues(int m, float parameter[]) {
116          float[] fit = new float[m + 1];
117          for (int i = 1; i <= m; i++) {
118              fit[i] = f(x[i], parameter);
119          };
120          return fit;
121      }
122    
123      boolean inRange(float x, float parameter[]) {
124    
125          // TODO: constrain the parameters 
126          //  5  < parameter[SIGMA] < 500   (or within 50% of maximum wavelet scale width)
127          //  0.5 * peak in wavelet space   < parameter[HEIGHT] < 1.5 * peak in wavelet space
128          //  0.8 *  peak position in wavelet space < parameter[CENTER] < peak position in wavelet space
129    
130          if (parameter[SIGMA] == 0f || parameter[HEIGHT] == 0f) return false;
131          else return true;
132      }
133    
134      // BUG: this is probably wrong !
135      public boolean computeResidualVector(int m, int n, float parameter[], float rv[]) {
136          for (int i = 1; i <= m; i++) {
137              if (!inRange(x[i], parameter)) return false;
138              rv[i] = f(x[i], parameter) - y[i];
139          }
140          return true;
141      }
142    
143      // BUG: this is probably wrong (but it seems right, checked twice) !
144      public void computeJacobian(int m, int n, float parameter[], float rv[], float jac[][]) {
145          for (int i = 1; i <= m; i++) {
146              jac[i][HEIGHT] = f(x[i], parameter) / parameter[HEIGHT];
147              jac[i][CENTER] = f(x[i], parameter) * (x[i] - parameter[HEIGHT]) / (parameter[SIGMA] * parameter[SIGMA]);
148              jac[i][SIGMA]  = jac[i][CENTER] * (x[i] - parameter[HEIGHT]) / parameter[SIGMA];
149          }
150      }
151    
152      public int getInvocations() {return 75;};
153      public float getMachinePrecision()  {return 1e-6f;};
154      public float getRelativeTolerance() {return 1e-4f;};
155      public float getAbsoluteTolerance() {return 1e-1f;};
156      public float getStartingValue()     {return 1e-2f;};
157    }