com.actelion.research.util.NumericalRecipes Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of openchemlib Show documentation
Show all versions of openchemlib Show documentation
Open Source Chemistry Library
/*
* Copyright (c) 1997 - 2016
* Actelion Pharmaceuticals Ltd.
* Gewerbestrasse 16
* CH-4123 Allschwil, Switzerland
*
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
* 3. Neither the name of the the copyright holder nor the
* names of its contributors may be used to endorse or promote products
* derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/*
* @(#)NumericalRecipes.java
*
* contains routines from 'Numerical Recipes in C' ported from C to Java
*
* @original authors W.H.Press, B.P.Fannery, S.A.Teukolsky, W.T.Vettering
*/
package com.actelion.research.util;
public class NumericalRecipes {
private double _fit_siga,_fit_sigb,_fit_a,_fit_b,_fit_chi2,_fit_q; // by reference params of fit()
private double _gser_gamser, _gser_gln; // by reference params of gser()
private double _gcf_gammcf, _gcf_gln; // by reference params of gcf()
private double ochisq; // used by mrqmin
private double[] atry,beta,da;
private double[][] oneda;
private int mfit;
public static void svbksb(double[][] u, double[] w, double[][] v, int m, int n, double[] b, double[] x) {
// Solves A·X = B for a vector X, where A is specified by the arrays u[1..m][1..n], w[1..n],
// v[1..n][1..n] as returned by svdcmp. m and n are the dimensions of a, and will be equal for
// square matrices. b[1..m] is the input right-hand side. x[1..n] is the output solution vector.
// No input quantities are destroyed, so the routine may be called sequentially with different b's.
double[] tmp = new double[n];
for (int j=0; j=0; i--) { // Accumulation of right-hand transformations.
if (i < n-1) {
if (g != 0.0) {
for (int j=l; j=0; i--) { // Accumulation of left-hand transformations.
l=i+1;
g=w[i];
for (int j=l; j=0; k--) { // Diagonalization of the bidiagonal form: Loop over
// singular values, and over allowed iterations.
for (int its=1;its<=30;its++) {
boolean flag = true;
int nm=0;
for (l=k; l>=0; l--) { // Test for splitting.
nm=l-1; // Note that rv1[1] is always zero.
if ((double)(Math.abs(rv1[l])+anorm) == anorm) {
flag = false;
break;
}
if ((double)(Math.abs(w[nm])+anorm) == anorm)
break;
}
if (flag) {
c=0.0; // Cancellation of rv1[l], if l > 1.
s=1.0;
for (int i=l; i2)
_fit_q = gammq(0.5*(ndata-2),0.5*(_fit_chi2)); // Equation (15.2.12).
}
}
public void mrqmin(double x[], double y[], double sig[], int ndata, double a[], int ia[],
int ma, double[][] covar, double[][] alpha, FittingFunction funcs) throws Exception {
// Levenberg-Marquardt method, attempting to reduce the value x^2 of a fit between a set of data
// points x[1..ndata], y[1..ndata] with individual standard deviations sig[1..ndata],
// and a nonlinear function dependent on ma coeficients a[1..ma]. The input array ia[1..ma]
// indicates by nonzero entries those components of a that should be fitted for, and by zero
// entries those components that should be held fixed at their input values. The program returns
// current best-fit values for the parameters a[1..ma], and x^2 = chisq. The arrays
// covar[1..ma][1..ma], alpha[1..ma][1..ma] are used as working space during most
// iterations. Supply a routine funcs(x,a,yfit,dyda,ma) that evaluates the fitting function
// yfit, and its derivatives dyda[1..ma] with respect to the fitting parameters a at x. On
// the first call provide an initial guess for the parameters a, and set alamda<0 for initialization
// (which then sets alamda=.001). If a step succeeds chisq becomes smaller and alamda decreases
// by a factor of 10. If a step fails alamda grows by a factor of 10. You must call this
// routine repeatedly until convergence is achieved. Then, make one final call with alamda=0, so
// that covar[1..ma][1..ma] returns the covariance matrix, and alpha the curvature matrix.
// (Parameters held fixed will return zero covariances.)
if (funcs.alamda < 0.0) { /* Initialization. */
atry = new double[ma];
beta = new double[ma];
da = new double[ma];
mfit = 0;
for (int j=0;j ITMAX)
nrerror("a too large, ITMAX too small in gcf");
_gcf_gammcf = Math.exp(-x+a*Math.log(x)-(_gcf_gln))*h; // Put factors in front.
}
private double gammq(double a, double x) throws Exception {
// Returns the incomplete gamma function Q(a, x) \uFFFDß 1 \uFFFD| P(a, x).
if (x < 0.0 || a <= 0.0)
nrerror("Invalid arguments in routine gammq");
if (x < (a+1.0)) { // Use the series representation
gser(a, x);
return 1.0 - _gser_gamser; // and take its complement.
}
else { // Use the continued fraction representation.
gcf(a, x);
return _gcf_gammcf;
}
}
private double gammln(double xx) {
// Returns the value ln[Ã(xx)] for xx > 0.
// Internal arithmetic will be done in double precision,
// a nicety that you can omit if .ve-.gure
// accuracy is good enough.
final double[] cof = { 76.18009172947146,
-86.50532032941677,
24.01409824083091,
-1.231739572450155,
0.1208650973866179e-2,
-0.5395239384953e-5 };
double x = xx;
double y = xx;
double tmp = x + 5.5;
tmp -= (x+0.5) * Math.log(tmp);
double ser = 1.000000000190015;
for (int j=0; j<=5; j++)
ser += cof[j] / ++y;
return -tmp + Math.log(2.5066282746310005*ser/x);
}
private static void covsrt(double covar[][], int ma, int ia[], int mfit) {
// Expand in storage the covariance matrix covar, so as to take into account parameters
// that are being held fixed. (For the latter, return zero covariances.)
for (int i=mfit; i=0;j--) {
if (ia[j] != 0) {
for (int i=0;i= big) {
big=Math.abs(a[j][k]);
irow=j;
icol=k;
}
}
}
++(ipiv[icol]);
/* We now have the pivot element, so we interchange rows, if needed, to put the pivot
element on the diagonal. The columns are not physically interchanged, only relabeled:
indxc[i], the column of the ith pivot element, is the ith column that is reduced, while
indxr[i] is the row in which that pivot element was originally located. If indxr[i] =
indxc[i] there is an implied column interchange. With this form of bookkeeping, the
solution b’s will end up in the correct order, and the inverse matrix will be scrambled
by columns. */
if (irow != icol) {
for (int l=0;l=0;l--) {
if (indxr[l] != indxc[l])
for (int k=0;k absb) ? absa*Math.sqrt(1.0+SQR(absb/absa))
: ((absb == 0.0) ? 0.0 : absb*Math.sqrt(1.0+SQR(absa/absb)));
}
private static double SQR(double a) {
return (a == 0.0) ? 0.0 : a*a;
}
private static double SIGN(double a, double b) {
return (b >= 0.0) ? Math.abs(a) : -Math.abs(a);
}
private static void SWAP(double a, double b) {
double temp = a;
a = b;
b = temp;
}
private static void nrerror(String err) throws Exception {
throw new Exception("Numerical Recipes run-time error: "+err);
}
}