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 }