001/*
002 *   ModelData  -- A utility class used to model experimental data.
003 *
004 *   Copyright (C) 2002-2025, by Joseph A. Huwaldt. All rights reserved.
005 *   
006 *   This library is free software; you can redistribute it and/or
007 *   modify it under the terms of the GNU Lesser General Public
008 *   License as published by the Free Software Foundation; either
009 *   version 2 of the License, or (at your option) any later version.
010 *   
011 *   This library is distributed in the hope that it will be useful,
012 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
013 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
014 *   Lesser General Public License for more details.
015 *
016 *   You should have received a copy of the GNU Lesser General Public License
017 *   along with this program; if not, write to the Free Software
018 *   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
019 *   Or visit:  http://www.gnu.org/licenses/lgpl.html
020 */
021package jahuwaldt.tools.math;
022
023import java.util.Arrays;
024
025/**
026 * Methods in this class are used to create functions that model experimental data.
027 * 
028 * <p> Modified by: Joseph A. Huwaldt </p>
029 * 
030 * @author Joseph A. Huwaldt    Date:   October 21, 2002
031 * @version February 23, 2025
032 */
033public class ModelData {
034
035    //  Tolerance used in single-value-decomposition to remove singular values.
036    private static final double TOL = 1E-13;
037
038    /**
039     * Method that returns the coefficients of a polynomial of the specified degree that
040     * best models or "fits" the supplied data values in a minimization of chi-squared
041     * sense. This method assumes that each data point is normally distributed and has a
042     * standard deviation of 1.0.
043     *
044     * @param xarr Array of independent parameter data values to fit.
045     * @param yarr Array of dependent data values (associated with each X value) to be modeled.
046     * @param coef An existing array that will be filled in with the coefficients of the
047     *  polynomial, in decreasing power that best fits the sample data. Example: p(x) = A +
048     *  B*x + C*x^2 + D*x^3 + E*x^4 corresponds to coef[] = {A, B, C, D, E}. The number of
049     *  elements in the coef array determines the order of the polynomial created. Number
050     *  of elements be greater than 1.
051     * @return The chi-squared value of the polynomial returned. A value near zero
052     *  indicates a perfect fit.
053     * @throws RootException if there is a problem computing the polynomial.
054     */
055    public static double polynomial(double[] xarr, double[] yarr, double[] coef) throws RootException {
056        if (xarr == null)
057            throw new NullPointerException("xarr == null");
058        double[] sig = new double[xarr.length];
059        Arrays.fill(sig, 1.0);
060        BasisFunction bf = new PolynomialFit();
061        return fit(xarr, yarr, sig, bf, coef);
062    }
063
064    /**
065     * Method that returns the coefficients of a polynomial of the specified degree that
066     * best models or "fits" the supplied data values in a minimization of chi-squared
067     * sense.
068     *
069     * @param xarr Array of independent parameter data values to fit.
070     * @param yarr Array of dependent data values (associated with each X value) to be modeled.
071     * @param sig Array of individual standard deviations for each data point.
072     * @param coef An existing array that will be filled in with the coefficients of the
073     *  polynomial, in decreasing power that best fits the sample data. Example: p(x) = A +
074     *  B*x + C*x^2 + D*x^3 + E*x^4 corresponds to coef[] = {A, B, C, D, E}. The number of
075     *  elements in the coef array determines the order of the polynomial created. Number
076     *  of elements be greater than 1.
077     * @return The chi-squared value of the polynomial returned. A value near zero
078     *  indicates a perfect fit.
079     * @throws RootException if there is a problem computing the polynomial.
080     */
081    public static double polynomial(double[] xarr, double[] yarr, double[] sig, double[] coef) throws RootException {
082        BasisFunction bf = new PolynomialFit();
083        return fit(xarr, yarr, sig, bf, coef);
084    }
085
086    /**
087     * Method that returns the coefficients of an arbitrary basis function that best
088     * models or "fits" the supplied data values in a minimization of chi-squared sense.
089     *
090     * @param xarr Array of independent parameter data values to fit.
091     * @param yarr Array of dependent data values (associated with each X value) to be
092     * modeled.
093     * @param sig Array of individual standard deviations for each data point.
094     * @param func The basis function used to generate the data fit.
095     * @param coef An existing array that will be filled in with the coefficients of the
096     * basis function that best fits the data.
097     * @return The chi-squared value of the function is returned. A value near zero
098     * indicates a perfect fit.
099     * @throws RootException if there is a problem computing the fit.
100     */
101    public static double fit(double[] xarr, double[] yarr, double[] sig, BasisFunction func, double[] coef) throws RootException {
102        if (xarr == null)
103            throw new NullPointerException("xarr == null");
104        if (yarr == null)
105            throw new NullPointerException("yarr == null");
106        if (sig == null)
107            throw new NullPointerException("sig == null");
108        if (func == null)
109            throw new NullPointerException("func == null");
110        if (coef == null)
111            throw new NullPointerException("coef == null");
112
113        int ma = coef.length;
114        if (ma < func.getMinNumCoef())
115            throw new IllegalArgumentException("Number of coefficients must be greater than " + (func.getMinNumCoef() - 1));
116
117        int ndata = xarr.length;
118        if (ndata != yarr.length)
119            throw new IllegalArgumentException("xarr and yarr must be the same length.");
120        if (ndata != sig.length)
121            throw new IllegalArgumentException("sig must be the same length as xarr and yarr.");
122
123        //      Allocate the arrays that we are going to need.
124        double[][] u = new double[ndata][ma];
125        double[][] v = new double[ma][ma];
126        double[] w = new double[ma];
127
128        //      Use Singular Value Decomposition to find the parameters of the basis function.
129        double chisq = svdfit(xarr, yarr, sig, coef, u, v, w, func);
130
131        return chisq;
132    }
133
134    /**
135     * Returns the slope of the line formed by a linear regression through the specified
136     * data arrays.
137     *
138     * @param x The independent array of data points.
139     * @param y The dependent array of data points (must be the same size as the x array.
140     * @return The slope of the line formed by a linear regression through the x and y
141     *  arrays of data points.
142     */
143    public static double linearSlope(double[] x, double[] y) {
144
145        //  Find the mean of the input arrays.
146        int length = x.length;
147        double sumX = 0, sumY = 0;
148        for (int i = 0; i < length; ++i) {
149            sumX += x[i];
150            sumY += y[i];
151        }
152        double xbar = sumX / length;
153        double ybar = sumY / length;
154
155        //  b = sum( (x - xbar)*(y - ybar) ) / sum ( (x - xbar)^2 )
156        double denom = 0;
157        double num = 0;
158        for (int i = 0; i < 4; ++i) {
159            double a = x[i] - xbar;
160            num += a * (y[i] - ybar);
161            denom += a * a;
162        }
163        double b = num / denom;
164
165        return b;
166    }
167
168    /**
169     * Used to test the methods in this class.
170     *
171     * @param args the command-line arguments
172     */
173    public static void main(String[] args) {
174
175        System.out.println();
176        System.out.println("Testing ModelData...");
177
178        try {
179            System.out.println("  Poly Curve Fit #1:  linear");
180            double[] x = {0, 1, 2, 3, 4, 5};
181            double[] y = {1, 3, 5, 7, 9, 11};
182            double[] coef = new double[2];
183            double chisq = polynomial(x, y, coef);
184            System.out.println("    coef = " + (float)coef[0] + ", " + (float)coef[1] + ", chisq = " + (float)chisq);
185            System.out.println("      should be:  p(x) = 1 + 2*x\n");
186
187            System.out.println("  Poly Curve Fit #2:  quadratic");
188            double[] x1 = {3, 4, 5};
189            double[] y1 = {-13.87, -0.09, -13.95};
190            double[] coef1 = new double[3];
191            chisq = polynomial(x1, y1, coef1);
192            System.out.println("    coef = " + (float)coef1[0] + ", " + (float)coef1[1] + ", " + (float)coef1[2] + ", chisq = " + (float)chisq);
193            System.out.println("      should be:  p(x) = -221.05 + 110.52*x - 13.82*x^2");
194
195        } catch (Exception e) {
196            e.printStackTrace();
197        }
198
199        System.out.println("Done.");
200
201    }
202
203    //-----------------------------------------------------------------------------------
204    /**
205     * Given a set of data points x,y with individual standard deviations sig, use chi^2
206     * minimization to determine the coefficients "a" of the supplied basis fitting
207     * function y = sumi(ai*funci(x)).
208     */
209    private static double svdfit(double[] x, double[] y, double[] sig, double[] a, double[][] u, double[][] v, double[] w,
210            BasisFunction f) throws RootException {
211        int ndata = x.length;
212        int ma = a.length;
213
214        double[] beta = new double[ndata];
215        double[] afunc = new double[ma];
216
217        //      Accumulate coefficients of the fitting matrix.
218        for (int i = 0; i < ndata; ++i) {
219            f.func(x[i], afunc);
220            double sigi = 1.0 / sig[i];
221            for (int j = 0; j < ma; ++j) {
222                u[i][j] = afunc[j] * sigi;
223            }
224            beta[i] = y[i] * sigi;
225        }
226
227        //      Singular value decomposition.
228        svdcmp(u, w, v);
229
230        //      Edit the singular values.
231        double wmax = 0;
232        for (int j = 0; j < ma; ++j) {
233            if (w[j] > wmax)
234                wmax = w[j];
235        }
236        double thresh = TOL * wmax;
237        for (int j = 0; j < ma; ++j) {
238            if (w[j] < thresh)
239                w[j] = 0;
240        }
241
242        svdksb(u, w, v, beta, a);
243
244        //      Evaluate chi-square.
245        double chisq = 0;
246        for (int i = 0; i < ndata; ++i) {
247            f.func(x[i], afunc);
248            double sum = 0;
249            for (int j = 0; j < ma; ++j) {
250                sum += a[j] * afunc[j];
251            }
252            double tmp = (y[i] - sum) / sig[i];
253            chisq += tmp * tmp;
254        }
255
256        return chisq;
257    }
258
259    /**
260     * Solves A*X = B for a vector X, where A is specified by the arrays u[][], w[], and
261     * v[][] as returned by svdcmp(). b is the right-hand side. x is the output solution
262     * vector. No input quantities are destroyed, so the routine can be called
263     * sequentially with different b's.
264     */
265    private static void svdksb(double[][] u, double[] w, double[][] v, double b[], double x[]) {
266        int n = u[0].length;
267        int m = u.length;
268
269        double[] tmp = new double[n];
270        for (int j = 0; j < n; ++j) {
271            //  Calculate U^T*B.
272            double s = 0;
273            if (w[j] != 0) {
274                //      Nonzero result only if w[j] is nonzero.
275                for (int i = 0; i < m; ++i) {
276                    s += u[i][j] * b[i];
277                }
278                s /= w[j];
279            }
280            tmp[j] = s;
281        }
282
283        //      Matrix multiply by V to get answer.
284        for (int j = 0; j < n; ++j) {
285            double s = 0;
286            for (int jj = 0; jj < n; ++jj) {
287                s += v[j][jj] * tmp[jj];
288            }
289            x[j] = s;
290        }
291    }
292
293    /**
294     * Given a matrix a[][], this routine computes its singular value decomposition. A =
295     * U*W*V^T. The matrix U replaces "a" on output. The diagonal matrix of singular
296     * values W is output as a vector w[]. The matrix V (not the transpose V^T) is output
297     * as v[][].
298     */
299    private static void svdcmp(double[][] a, double[] w, double[][] v) throws RootException {
300        int n = a[0].length;
301        int m = a.length;
302        int k, l = 0;
303
304        double[] rv1 = new double[n];
305        double g = 0;
306        double scale = 0;
307        double anorm = 0;
308        double s = 0;
309
310        //      Householder reduction to bidiagonal form.
311        for (int i = 0; i < n; ++i) {
312            l = i + 1;
313            rv1[i] = scale * g;
314            s = g = scale = 0;
315
316            if (i < m) {
317                for (k = i; k < m; ++k) {
318                    scale += Math.abs(a[k][i]);
319                }
320                if (scale != 0) {
321                    for (k = i; k < m; ++k) {
322                        a[k][i] /= scale;
323                        s += a[k][i] * a[k][i];
324                    }
325                    double f = a[i][i];
326                    g = -MathTools.sign(Math.sqrt(s), f);
327                    double h = f * g - s;
328                    a[i][i] = f - g;
329                    for (int j = l; j < n; ++j) {
330                        for (s = 0, k = i; k < m; ++k) {
331                            s += a[k][i] * a[k][j];
332                        }
333                        f = s / h;
334                        for (k = i; k < m; ++k) {
335                            a[k][j] += f * a[k][i];
336                        }
337                    }
338                    for (k = i; k < m; ++k) {
339                        a[k][i] *= scale;
340                    }
341                }
342
343            }
344
345            w[i] = scale * g;
346            g = s = scale = 0;
347
348            if (i < m && i + 1 != n) {
349                for (k = l; k < n; ++k) {
350                    scale += Math.abs(a[i][k]);
351                }
352                if (scale != 0) {
353                    for (k = l; k < n; ++k) {
354                        a[i][k] /= scale;
355                        s += a[i][k] * a[i][k];
356                    }
357                    double f = a[i][l];
358                    g = -MathTools.sign(Math.sqrt(s), f);
359                    double h = f * g - s;
360                    a[i][l] = f - g;
361                    for (k = l; k < n; ++k) {
362                        rv1[k] = a[i][k] / h;
363                    }
364                    for (int j = l; j < m; ++j) {
365                        for (s = 0, k = l; k < n; ++k) {
366                            s += a[j][k] * a[i][k];
367                        }
368                        for (k = l; k < n; ++k) {
369                            a[j][k] += s * rv1[k];
370                        }
371                    }
372                    for (k = l; k < n; ++k) {
373                        a[i][k] *= scale;
374                    }
375                }
376
377            }
378
379            anorm = Math.max(anorm, (Math.abs(w[i]) + Math.abs(rv1[i])));
380        }       //      Next i
381
382        //      Accumulation of right-hand transformations.
383        for (int i = n - 1; i >= 0; --i) {
384            if (i + 1 < n) {
385                if (g != 0) {
386                    //  Double division to avoid possible underflow.
387                    for (int j = l; j < n; ++j) {
388                        v[j][i] = (a[i][j] / a[i][l]) / g;
389                    }
390                    for (int j = l; j < n; ++j) {
391                        for (s = 0, k = l; k < n; ++k) {
392                            s += a[i][k] * v[k][j];
393                        }
394                        for (k = l; k < n; ++k) {
395                            v[k][j] += s * v[k][i];
396                        }
397                    }
398                }
399                for (int j = l; j < n; ++j) {
400                    v[i][j] = v[j][i] = 0;
401                }
402            }
403
404            v[i][i] = 1;
405            g = rv1[i];
406            l = i;
407        }       //      Next i
408
409        //      Accumulation of left-hand transformations.
410        for (int i = Math.min(m, n) - 1; i >= 0; --i) {
411            l = i + 1;
412            g = w[i];
413            for (int j = l; j < n; ++j) {
414                a[i][j] = 0;
415            }
416            if (g != 0) {
417                g = 1.0 / g;
418                for (int j = l; j < n; ++j) {
419                    for (s = 0, k = l; k < m; ++k) {
420                        s += a[k][i] * a[k][j];
421                    }
422                    double f = (s / a[i][i]) * g;
423                    for (k = i; k < m; ++k) {
424                        a[k][j] += f * a[k][i];
425                    }
426                }
427                for (int j = i; j < m; ++j) {
428                    a[j][i] *= g;
429                }
430
431            } else
432                for (int j = i; j < m; ++j) {
433                    a[j][i] = 0;
434                }
435
436            ++a[i][i];
437        }       //      Next i
438
439        //      Diagonalization of the bidiagonal form:  Loop over singular values,
440        //      and over allowed iterations.
441        for (k = n - 1; k >= 0; --k) {
442            for (int its = 0; its < 30; ++its) {
443                boolean flag = true;
444                double x, y, z, c, f, h;
445                int nm = 0;
446
447                //      Test for splitting.
448                for (l = k; l >= 0; --l) {
449                    nm = l - 1;                         //      Note that rv1[1] is always zero.
450                    if ((Math.abs(rv1[l]) + anorm) == anorm) {
451                        flag = false;
452                        break;
453                    }
454                    if ((Math.abs(w[nm]) + anorm) == anorm)
455                        break;
456                }
457
458                //      Cancelation of rv1[l], if l > 1.
459                if (flag) {
460                    c = 0;
461                    s = 1;
462                    for (int i = l; i <= k; ++i) {
463                        f = s * rv1[i];
464                        rv1[i] = c * rv1[i];
465                        if ((Math.abs(f) + anorm) == anorm)
466                            break;
467                        g = w[i];
468                        h = pythag(f, g);
469                        w[i] = h;
470                        h = 1 / h;
471                        c = g * h;
472                        s = -f * h;
473                        for (int j = 0; j < m; ++j) {
474                            y = a[j][nm];
475                            z = a[j][i];
476                            a[j][nm] = y * c + z * s;
477                            a[j][i] = z * c - y * s;
478                        }
479                    }
480                }
481
482                z = w[k];
483                if (l == k) {
484                    //  Convergence
485                    //  Singular value is made nonnegative.
486                    if (z < 0) {
487                        w[k] = -z;
488                        for (int j = 0; j < n; ++j) {
489                            v[j][k] = -v[j][k];
490                        }
491                    }
492                    break;
493                }
494
495                if (its == 30)
496                    throw new RootException("No convergence in 30 svdcmp() iterations.");
497
498                //      Shift from bottom 2-by-2 minor.
499                x = w[l];
500                nm = k - 1;
501                y = w[nm];
502                g = rv1[nm];
503                h = rv1[k];
504                f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2 * h * y);
505                g = pythag(f, 1);
506                f = ((x - z) * (x + z) + h * ((y / (f + MathTools.sign(g, f))) - h)) / x;
507                c = s = 1;
508
509                //      Next QR transformation.
510                for (int j = l; j <= nm; ++j) {
511                    int i = j + 1;
512                    g = rv1[i];
513                    y = w[i];
514                    h = s * g;
515                    g = c * g;
516                    z = pythag(f, h);
517                    rv1[j] = z;
518                    c = f / z;
519                    s = h / z;
520                    f = x * c + g * s;
521                    g = g * c - x * s;
522                    h = y * s;
523                    y *= c;
524                    for (int jj = 0; jj < n; ++jj) {
525                        x = v[jj][j];
526                        z = v[jj][i];
527                        v[jj][j] = x * c + z * s;
528                        v[jj][i] = z * c - x * s;
529                    }
530                    z = pythag(f, h);
531                    w[j] = z;
532
533                    //  Rotation can be arbitrary if z == 0.
534                    if (z != 0) {
535                        z = 1 / z;
536                        c = f * z;
537                        s = h * z;
538                    }
539                    f = c * g + s * y;
540                    x = c * y - s * g;
541                    for (int jj = 0; jj < m; ++jj) {
542                        y = a[jj][j];
543                        z = a[jj][i];
544                        a[jj][j] = y * c + z * s;
545                        a[jj][i] = z * c - y * s;
546                    }
547
548                }       //      Next j
549
550                rv1[l] = 0;
551                rv1[k] = f;
552                w[k] = x;
553
554            }   //      Next its
555        }       //      Next k
556
557    }
558
559    /**
560     * Method that calculates sqrt(a^2 + b^2) without destructive underflow or overflow.
561     */
562    private static double pythag(double a, double b) {
563        double absa = Math.abs(a);
564        double absb = Math.abs(b);
565
566        if (absa > absb) {
567            double ratio = absb / absa;
568            return absa * Math.sqrt(1 + ratio * ratio);
569
570        } else {
571            if (absb == 0)
572                return 0.0;
573            else {
574                double ratio = absa / absb;
575                return absb * Math.sqrt(1 + ratio * ratio);
576            }
577        }
578    }
579}