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    }