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 }