001 package numal.lowprecision;
002
003 import java.util.*;
004 import static java.lang.Math.*;
005 import static numal.lowprecision.Basic.*;
006 import static numal.lowprecision.LinearAlgebra.qrisngvaldec;
007
008 public class AnalyticProblems extends Object {
009
010 /**
011 * The standard deviation of each parameter in the fit is sqrt(v[j][j] / m)
012 * and the correlation between parameter i and parameter j is
013 * v[i][j]/ sqrt(v[i][i]*v[j][j]) where i=1,...,n and j=i+1,...,n.
014 * // http://www.crcpress.com/shopping_cart/products/product_detail.asp?sku=C4304
015 * <ul>
016 * <lh>Output values</lh>
017 * <li>out[1] information about the termination of the process
018 * <li>0 : normal termination</li>
019 * <li>1 : number of invocations of {@link Function#computeResidualVector} exceeded {@link Function#getInvocations()}</li>
020 * <li>2 : {@link Function#computeResidualVector} returned false</li>
021 * <li>3 : {@link Function#computeResidualVector} returned false when called with the initial estimates of parameter[1:n];
022 * the iteration process was not started and so v[1:n][1:n] can't be used</li>
023 * <li>4 : precision asked for can't be attained; this precision is possibly chosen too high, relative to the
024 * precision in which the residual vector is calculated (see {@link Function#getRelativeTolerance()}</li>
025 * </li>
026 * <li>out[2] the Euclidean norm of the residual vector calculated with values of the unknowns delivered</li>
027 * <li>out[3] the Euclidean norm of the residual vector calculated with the initial values of the unknown variables</li>
028 * <li>out[4] the number of calls of {@link Function#computeResidualVector} necessary to attain the calculated results
029 * </li>
030 * <li>out[5] the total number of iterations performed; note that in each iteration one call to {@link Function#computeJacobian} had to be made</li>
031 * <li>out[6] the improvement vector in the last iteration</li>
032 * <li>out[7] the condition number of J^T J, i.e., the ratio of its largest to smallest eigenvalues</li>
033 * </ul>
034 *
035 * Dependencies: {@link Basic#mulcol},
036 * {@link Basic#dupvec},
037 * {@link Basic#vecvec},
038 * {@link Basic#matvec},
039 * {@link Basic#tamvec},
040 * {@link Basic#mattam},
041 * {@link LinearAlgebra#qrisngvaldec}.
042 *
043 * @param m sample points used in the fitting
044 * @param n number of parameters, must be less than or equal to m
045 * @param parameter <b>input</b>: initial approximation to the set of parameters <b>output</b>: parameters producing a least square fit
046 * @param rv 1xm vector <b>output</b> : residual vector obtained with current value of unknowns
047 * @param v nxn matrix <b>output</b> : inverse of matrix J^T J where J denotes the transpose of the matrix of partial derivatives dg[1..m]/dp[1..n]
048 * @param function residual vector of a given function, jacobian and options such as precision and maximum number of iterations
049 * @param out an array of 7 output values (see description above)
050 */
051 public static void marquardt(int m, int n, float parameter[], float rv[], float v[][], Function function, float out[]) {
052 int invocation = 1;
053 int errorCode = 0;
054 int iteration = 1;
055 if (!function.computeResidualVector(m, n, parameter, rv)) {
056 errorCode = 3;
057 out[4] = invocation;
058 out[5] = iteration - 1;
059 out[1] = errorCode;
060 return;
061 }
062 final float em[] = { function.getMachinePrecision(), 0f, function.getMachinePrecision(), 0f, 10 * n, 0f, function.getMachinePrecision(), 0f}; // 8 entries
063 final float mu = 0.01f;
064 final float w = 0.5f;
065 final float vv = 10f;
066 final float ww = (function.getStartingValue() < 1.0e-7) ? 1e-8f : 1e-1f * function.getStartingValue();
067 final float reltolres = function.getRelativeTolerance();
068 final float abstolres = function.getAbsoluteTolerance() * function.getAbsoluteTolerance();
069 final float pw = (float) (-log(ww * function.getMachinePrecision()) / 2.3d);
070 float fparpres;
071 float val[] = new float[n + 1];
072 float b[] = new float[n + 1];
073 float bb[] = new float[n + 1];
074 float parpres[] = new float[n + 1];
075 float jac[][] = new float[m + 1][n + 1];
076 float p = 0f;
077 float res = 0f;
078 boolean emergency = false;
079 float lambda = 0f;
080 float lambdamin;
081 float fpar = vecvec(1, m, 0, rv, rv);
082 out[3] = (float) sqrt(fpar);
083
084 do {
085 function.computeJacobian(m, n, parameter, rv, jac);
086 qrisngvaldec(jac, m, n, val, v, em);
087 if (iteration == 1)
088 lambda = function.getStartingValue() * vecvec(1, n, 0, val, val);
089 else
090 if (p == 0f) lambda *= w;
091 for (int i = 1; i <= n; i++) b[i] = val[i] * tamvec(1, m, i, jac, rv);
092 while (true) {
093 for (int i = 1; i <= n; i++) bb[i] = b[i] / (val[i] * val[i] + lambda);
094 for (int i = 1; i <= n; i++) parpres[i] = parameter[i] - matvec(1, n, i, v, bb);
095 invocation++;
096 if (invocation >= function.getInvocations())
097 errorCode = 1;
098 else
099 if (!function.computeResidualVector(m, n, parpres, rv)) errorCode = 2;
100 if (errorCode != 0) {
101 emergency = true;
102 break;
103 }
104 fparpres = vecvec(1, m, 0, rv, rv);
105 res = fpar - fparpres;
106 if (res < mu * vecvec(1, n, 0, b, bb)) {
107 p += 1f;
108 lambda *= vv;
109 if (p == 1f) {
110 lambdamin = ww * vecvec(1, n, 0, val, val);
111 if (lambda < lambdamin) lambda = lambdamin;
112 }
113 if (p >= pw) {
114 errorCode = 4;
115 emergency = true;
116 break;
117 }
118 } else {
119 dupvec(1, n, 0, parameter, parpres);
120 fpar = fparpres;
121 break;
122 }
123 System.out.println(iteration + " " + parameter[1] + " " + parameter[2] + " " + parameter[3]);
124 }
125 if (emergency) break;
126 iteration++;
127
128 } while (fpar > abstolres && res > reltolres * fpar + abstolres);
129
130 for (int i = 1; i <= n; i++)
131 mulcol(1, n, i, i, jac, v, 1f / (val[i] + function.getMachinePrecision()));
132 for (int i = 1; i <= n; i++)
133 for (int j = 1; j <= i; j++)
134 v[i][j] = v[j][i] = mattam(1, n, i, j, jac, jac);
135
136 lambda = lambdamin = val[1];
137 for (int i = 2; i <= n; i++)
138 if (val[i] > lambda) lambda = val[i];
139 else if (val[i] < lambdamin) lambdamin = val[i];
140
141 float temp = lambda / (lambdamin + function.getMachinePrecision());
142 out[7] = temp * temp;
143 out[2] = (float) sqrt(fpar);
144 out[6] = (float) sqrt(res + fpar) - out[2];
145 out[4] = invocation;
146 out[5] = iteration - 1;
147 out[1] = errorCode;
148 }
149
150
151 }