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 }