/*
* Open Source Physics software is free software as described near the bottom of this code file.
*
* For additional information and documentation on Open Source Physics please see:
* <http://www.opensourcephysics.org/>
*/
/*
@Author J E Hasbun 2007.
Applies the Hessian method to find the parameters that minimize a
function of those parameters also but uses LUPDecomposition's solve
method instead of the inverse to get the new guesses
@Copyright (c) 2007
This software is to support the Open Source Physics library
http://www.opensourcephysics.org under the terms of the GNU General Public
License (GPL) as published by the Free Software Foundation.
*/
package org.opensourcephysics.numerics;
/**
* Class description
*
*/
public class HessianMinimize {
int Iterations;
double[][] H;
double[] xp;
double[] xm;
double[] xpp;
double[] xpm;
double[] xmp;
double[] xmm;
private double rmsd_tmp, rmsd;
private double[] xtmp;
/* Inputs
Veq - the function of m parameters whose minimum is sought
x - the array containing the guess to the solutions
max - the maximum iteration number
tol - the tolerance level
*/
public double minimize(MultiVarFunction Veq, double[] x, int max, double tol) {
int m = x.length;
double[] xxn = new double[m];
double[] D = new double[m];
double[] dx = new double[m];
xtmp = new double[m];
System.arraycopy(x, 0, xtmp, 0, m);
rmsd_tmp = Veq.evaluate(x);
rmsd = 0;
crudeGuess(Veq, x); //obtain a crude guess
check_rmsd(Veq, xtmp, x, m); //check if x is better, else keep old one
for(int i = 0; i<m; i++) {
dx[i] = (Math.abs(x[i])+1.0)/1e5; //step sizes for the finite differences
}
double err, relerr;
err = 9999.;
relerr = 9999.;
//Use the Hessian method for an equation of several variables
//start with a good guess.
Iterations = 0;
while((err>tol*1.e-6)&&(relerr>tol*1.e-6)&&(Iterations<max)) {
Iterations++;
LUPDecomposition lu = new LUPDecomposition(getHessian(Veq, x, D, dx));
// use the LUPDecomposition's solve method
xxn = lu.solve(D); //the corrections
for(int i = 0; i<m; i++) {
xxn[i] = xxn[i]+x[i]; //new guesses
}
err = (x[0]-xxn[0])*(x[0]-xxn[0]);
relerr = x[0]*x[0];
x[0] = xxn[0];
for(int i = 1; i<m; i++) {
err = err+(x[i]-xxn[i])*(x[i]-xxn[i]);
relerr = relerr+x[i]*x[i];
x[i] = xxn[i]; //copy to go back
//dx[i]=0.5*dx[i]; //could decrease the variation at each iteration
}
err = Math.sqrt(err); //the error
relerr = err/(relerr+tol);
}
check_rmsd(Veq, xtmp, x, m); //check if x is better, else keep old one
return err;
}
private void allocateArrays(int m) {
H = new double[m][m];
xp = new double[m];
xm = new double[m];
xpp = new double[m];
xpm = new double[m];
xmp = new double[m];
xmm = new double[m];
}
void crudeGuess(MultiVarFunction Veq, double[] x) {
/*
@Author J E Hasbun 2007.
This is a crude method to obtain better starting guesses for the
multiple funtion variable minimization problem. It uses a Naive
Newton Raphson step for each of the variables.
Ref: Computational Physics by P. L. Devries (J. Wiley, 1993)
@Copyright (c) 2007
This software is to support the Open Source Physics library
http://www.opensourcephysics.org under the terms of the GNU General Public
License (GPL) as published by the Free Software Foundation.
*/
double sp, s0, sm;
int m = x.length; //array size
int Nc = 5; //cycles to make
double f = 0.35; //moderates the step size to next crude guess
int n = 0; //cycle counter
double[] xp = new double[m];
double[] xm = new double[m];
double[] dx = new double[m];
for(int i = 0; i<m; i++) {
dx[i] = (Math.abs(x[i])+1.0)/1.e3; //step sizes for the finite derivatives
}
//Cycle through each parameter Nc times
while(n<Nc) {
n++;
for(int i = 0; i<m; i++) {
//The SUM will be evaluated with the parameters
//xp, a, and xm:
for(int k = 0; k<m; k++) {
if(k==i) {
xp[i] = x[i]+dx[i];
xm[i] = x[i]-dx[i];
} else {
xp[k] = x[k];
xm[k] = x[k];
}
}
sp = Veq.evaluate(xp); // Evaluate the sum.
s0 = Veq.evaluate(x);
sm = Veq.evaluate(xm);
//make the crude Newton-Raphson step next
x[i] = x[i]-f*0.5*dx[i]*(sp-sm)/(sp-2.0*s0+sm);
//As we move towards a minimum, we should decrease
//step size used in calculating the derivative.
dx[i] = 0.5*dx[i];
}
}
}
void check_rmsd(MultiVarFunction Veq, double[] xtmp, double[] xx, int mx) {
//checks whether xtmp or xx is better, and keep the better one
if(java.lang.Double.isNaN(ArrayLib.sum(xx))) {
rmsd = rmsd_tmp;
System.arraycopy(xtmp, 0, xx, 0, mx);
} else {
rmsd = Veq.evaluate(xx);
if(rmsd<=rmsd_tmp) {
rmsd_tmp = rmsd;
System.arraycopy(xx, 0, xtmp, 0, mx);
} else {
rmsd = rmsd_tmp;
System.arraycopy(xtmp, 0, xx, 0, mx);
}
}
}
public int getIterations() {
return Iterations;
}
public double[][] getHessian(MultiVarFunction Veq, double[] x, double[] D, double[] dx) {
/*
@Author J E Hasbun 2007.
Finds the Hessian of a function of several variables
The method is similar to the Newton method but for the derivative of a
general function. The idea is that to "minimize" a multivariable function
Veq({X}), we need the guess {Xnew}={Xold}+inverse{H}*{D}, so here we find
H and D for Veq. Ref: Computational Physics by P. L. Devries (J. Wiley, 1993)
Input:
Veq - the function of m parameters whose minimum is sought
x - the parameters
dx - the size of the variations used in the derivatives
Uses:
H - the square matrix containing the calculated Hessian as a return
the Hessian is a square matrix of 2nd order partial derivatives of
of Veq with respect to the variables
D - the 1st order negative partial derivative of Veq
builds the Hessian H - used in a multidimensional newton method
the arrays contain:
x - the variable parameters of function Veq
x - x[1], ..., x[i],..., x[m]
xp - x[1], ..., x[i]+dx[i],..., x[m]
xm - x[1], ..., x[i]-dx[i],..., x[m]
xpp - x[1],..,x[i]+dx[i],..., x[j]+dx[j],..., x[m]
xpm - x[1],..,x[i]+dx[i],..., x[j]-dx[j],..., x[m]
xmp - x[1],..,x[i]-dx[i],..., x[j]+dx[j],..., x[m]
xmm - x[1],..,x[i]-dx[i],..., x[j]-dx[j],..., x[m]
*/
//The Hessian H is calculated by the finite difference method
int m = x.length;
if((xp==null)||(xp.length!=m)) {
allocateArrays(m);
}
// Compute the Hessian:
for(int i = 0; i<m; i++) {
for(int j = i; j<m; j++) {
if(i==j) {
for(int k = 0; k<m; k++) { //reset the x's
xp[k] = x[k];
xm[k] = x[k];
}
xp[i] = x[i]+dx[i]; //change the ith one
xm[i] = x[i]-dx[i];
H[i][i] = (Veq.evaluate(xp)-2.0*Veq.evaluate(x)+Veq.evaluate(xm))/(dx[i]*dx[i]);
} else {
for(int k = 0; k<m; k++) { //reset the x's
xpp[k] = x[k];
xpm[k] = x[k];
xmp[k] = x[k];
xmm[k] = x[k];
}
xpp[i] = x[i]+dx[i]; //change the ith, jth ones
xpp[j] = x[j]+dx[j];
xpm[i] = x[i]+dx[i];
xpm[j] = x[j]-dx[j];
xmp[i] = x[i]-dx[i];
xmp[j] = x[j]+dx[j];
xmm[i] = x[i]-dx[i];
xmm[j] = x[j]-dx[j];
H[i][j] = ((Veq.evaluate(xpp)-Veq.evaluate(xpm))/(2.0*dx[j])-(Veq.evaluate(xmp)-Veq.evaluate(xmm))/(2.0*dx[j]))/(2.0*dx[i]);
H[j][i] = H[i][j];
}
}
}
// note the D function is the negative of the partial derivative
for(int i = 0; i<m; i++) {
for(int k = 0; k<m; k++) { //reset the x's
xp[k] = x[k];
xm[k] = x[k];
}
xp[i] = x[i]+dx[i]; //change the ith one
xm[i] = x[i]-dx[i];
D[i] = -(Veq.evaluate(xp)-Veq.evaluate(xm))/(2.0*dx[i]);
}
return H;
}
}
/*
* Open Source Physics software is free software; you can redistribute
* it and/or modify it under the terms of the GNU General Public License (GPL) as
* published by the Free Software Foundation; either version 2 of the License,
* or(at your option) any later version.
* Code that uses any portion of the code in the org.opensourcephysics package
* or any subpackage (subdirectory) of this package must must also be be released
* under the GNU GPL license.
*
* This software is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston MA 02111-1307 USA
* or view the license online at http://www.gnu.org/copyleft/gpl.html
*
* Copyright (c) 2007 The Open Source Physics project
* http://www.opensourcephysics.org
*/