001 package numal.lowprecision.test; 002 003 import numal.lowprecision.AnalyticProblems; 004 import numal.lowprecision.Function; 005 006 // http://www.crcpress.com/shopping_cart/products/product_detail.asp?sku=C4304 007 008 import java.text.DecimalFormat; 009 import static java.lang.Math.*; 010 011 /** Test the least-square fit method {@link AnalyticProblems#marquardt}. 012 * Should return the following results : 013 {@code 014 Parameters: 015 5.23219E2 -1.56833E2 -1.99784E-1 016 017 OUT: 018 7.22065E7 019 1.15716E2 020 1.74713E-3 021 1.65459E2 022 2.30000E1 023 2.20000E1 024 0.00000E0 025 026 Last residual vector: 027 -29.6 86.6 -47.3 -26.2 -22.9 39.5 028 } 029 */ 030 public class TestMarquardt extends Object implements Function { 031 032 /** Ordinate array corresponding to absica y (zero position not significant) */ 033 float[] x; 034 035 /** Array of values corresponding to x which are to be fitted with f(x,par) (zero position not significant) */ 036 float[] y; 037 038 public TestMarquardt(float[] x, float[] y) { 039 this.x = x; 040 this.y = y; 041 } 042 043 public static void main(String args[]) { 044 float[] parameter = {0f, 580f, -180f, -0.160f}; // 3 parameters 045 float[] x = {0f, -5f, -3f, -1f, 1f, 3f, 5f}; 046 float[] y = {0f, 127f, 151f, 379f, 421f, 460f, 426f}; 047 float rv[] = new float[x.length]; // residual vector f(x[i],par) - y[i] 048 final float jjinv[][] = new float[parameter.length][parameter.length]; 049 Function function = new TestMarquardt(x, y); 050 float[] out = new float[8]; 051 052 AnalyticProblems.marquardt(x.length - 1, parameter.length - 1, parameter, rv, jjinv, function, out); 053 054 DecimalFormat fiveDigit = new DecimalFormat("0.00000E0"); 055 DecimalFormat oneDigit = new DecimalFormat("0.0"); 056 System.out.println("===============Test Marquardt================"); 057 System.out.println("Parameters:\n " + 058 fiveDigit.format(parameter[1]) + " " + 059 fiveDigit.format(parameter[2]) + " " + 060 fiveDigit.format(parameter[3]) + "\n\nOUT:\n " + 061 fiveDigit.format(out[7]) + "\n " + 062 fiveDigit.format(out[2]) + "\n " + 063 fiveDigit.format(out[6]) + "\n " + 064 fiveDigit.format(out[3]) + "\n " + 065 fiveDigit.format(out[4]) + "\n " + 066 fiveDigit.format(out[5]) + "\n " + 067 fiveDigit.format(out[1]) + "\n\nLast residual vector:\n " + 068 oneDigit.format(rv[1]) + " " + 069 oneDigit.format(rv[2]) + " " + 070 oneDigit.format(rv[3]) + " " + 071 oneDigit.format(rv[4]) + " " + 072 oneDigit.format(rv[5]) + " " + 073 oneDigit.format(rv[6])); 074 } 075 076 077 078 /** Get the value of this function for a given abcissa and parameters. 079 * 080 * @param x the abcissa at which this function is to be evaluated 081 * @param parameter an array of function parameters 082 * @return the value of this function for the given x and parameters 083 */ 084 public static float f(float x, float parameter[]) { 085 return parameter[1] + parameter[2] * (float) exp(parameter[3] * x); 086 } 087 088 // implementation of the Function interface : 089 090 public float[] getAnalyticValues(int m, float parameter[]) { 091 float[] fit = new float[m]; 092 for (int i = 1; i <= m; i++) { 093 fit[i] = f(x[i], parameter); 094 }; 095 return fit; 096 } 097 098 public boolean computeResidualVector(int m, int n, float parameter[], float rv[]) { 099 for (int i = 1; i <= m; i++) { 100 if (!inRange(x[i], parameter)) return false; 101 rv[i] = f(x[i], parameter) - y[i]; 102 } 103 return true; 104 } 105 106 boolean inRange(float x, float parameter[]) { 107 if (parameter[3] * x > 680f) return false; 108 return true; 109 } 110 111 public void computeJacobian(int m, int n, float parameter[], float rv[], float jac[][]) { 112 for (int i = 1; i <= m; i++) { 113 jac[i][1] = 1f; 114 float ex = (float) exp(parameter[3] * x[i]); 115 jac[i][2] = ex; 116 jac[i][3] = x[i] * parameter[2] * ex; 117 } 118 } 119 120 public int getInvocations() {return 75;}; 121 public float getMachinePrecision() {return 1e-6f;}; 122 public float getRelativeTolerance() {return 1e-4f;}; 123 public float getAbsoluteTolerance() {return 1e-1f;}; 124 public float getStartingValue() {return 1e-2f;}; 125 }