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}