idwint.h   idwint.h 
skipping to change at line 57 skipping to change at line 57
#include "densesolver.h" #include "densesolver.h"
/************************************************************************* /*************************************************************************
IDW interpolant. IDW interpolant.
*************************************************************************/ *************************************************************************/
struct idwinterpolant struct idwinterpolant
{ {
int n; int n;
int nx; int nx;
int d; int d;
double r;
int nw; int nw;
kdtree tree; kdtree tree;
int modeltype;
ap::real_2d_array q; ap::real_2d_array q;
ap::real_1d_array xbuf; ap::real_1d_array xbuf;
ap::integer_1d_array tbuf; ap::integer_1d_array tbuf;
ap::real_1d_array rbuf; ap::real_1d_array rbuf;
ap::real_2d_array xybuf; ap::real_2d_array xybuf;
int debugsolverfailures; int debugsolverfailures;
double debugworstrcond; double debugworstrcond;
double debugbestrcond; double debugbestrcond;
}; };
skipping to change at line 86 skipping to change at line 88
Result: Result:
IDW interpolant Z(X) IDW interpolant Z(X)
-- ALGLIB -- -- ALGLIB --
Copyright 02.03.2010 by Bochkanov Sergey Copyright 02.03.2010 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
double idwcalc(idwinterpolant& z, const ap::real_1d_array& x); double idwcalc(idwinterpolant& z, const ap::real_1d_array& x);
/************************************************************************* /*************************************************************************
IDW interpolant using modified Shepard method IDW interpolant using modified Shepard method for uniform point
distributions.
INPUT PARAMETERS: INPUT PARAMETERS:
XY - X and Y values, array[0..N-1,0..NX]. XY - X and Y values, array[0..N-1,0..NX].
First NX columns contain X-values, last column contain First NX columns contain X-values, last column contain
Y-values. Y-values.
N - number of nodes, N>0. N - number of nodes, N>0.
NX - space dimension, NX>=1. NX - space dimension, NX>=1.
D - nodal function type, either: D - nodal function type, either:
* 0 constant model. Just for demonstration only, worst * 0 constant model. Just for demonstration only, worst
model ever. model ever.
* 1 linear model, least squares fitting. Simpe model for * 1 linear model, least squares fitting. Simpe model for
datasets too small for quadratic models datasets too small for quadratic models
* 2 quadratic model, least squares fitting. Best model * 2 quadratic model, least squares fitting. Best model
available (if your dataset is large enough). available (if your dataset is large enough).
* -1 "fast" linear model, use with caution!!! It is * -1 "fast" linear model, use with caution!!! It is
significantly faster than linear/quadratic and better significantly faster than linear/quadratic and better
than constant model. But it is less robust (especially than constant model. But it is less robust (especially
in the presence of noise). in the presence of noise).
NQ - number of points used to calculate nodal functions (ignored NQ - number of points used to calculate nodal functions (ignored
for constant models). NQ should be LARGER than 1.5 times for constant models). NQ should be LARGER than:
the number of coefficients in a nodal function: * max(1.5*(1+NX),2^NX+1) for linear model,
* 1.5*(1+NX) for linear model, * max(3/4*(NX+2)*(NX+1),2^NX+1) for quadratic model.
* 3/4*(NX+2)*(NX+1) for quadratic model.
Values less than this threshold will be silently increased. Values less than this threshold will be silently increased.
NW - number of points used to calculate weights and to interpolate. NW - number of points used to calculate weights and to interpolate.
Required: >=2^NX+1, values less than this threshold will be Required: >=2^NX+1, values less than this threshold will be
silently increased. silently increased.
Recommended value: about 2*NQ Recommended value: about 2*NQ
OUTPUT PARAMETERS: OUTPUT PARAMETERS:
Z - IDW interpolant. Z - IDW interpolant.
NOTES: NOTES:
skipping to change at line 131 skipping to change at line 133
models models
* when N is large, NQ and NW must be significantly smaller than N both * when N is large, NQ and NW must be significantly smaller than N both
to obtain optimal performance and to obtain optimal accuracy. In 2 or to obtain optimal performance and to obtain optimal accuracy. In 2 or
3-dimensional tasks NQ=15 and NW=25 are good values to start with. 3-dimensional tasks NQ=15 and NW=25 are good values to start with.
* NQ and NW may be greater than N. In such cases they will be * NQ and NW may be greater than N. In such cases they will be
automatically decreased. automatically decreased.
* this subroutine is always succeeds (as long as correct parameters are * this subroutine is always succeeds (as long as correct parameters are
passed). passed).
* see 'Multivariate Interpolation of Large Sets of Scattered Data' by * see 'Multivariate Interpolation of Large Sets of Scattered Data' by
Robert J. Renka for more information on this algorithm. Robert J. Renka for more information on this algorithm.
* this subroutine assumes that point distribution is uniform at the small
scales. If it isn't - for example, points are concentrated along
"lines", but "lines" distribution is uniform at the larger scale - then
you should use IDWBuildModifiedShepardR()
-- ALGLIB PROJECT -- -- ALGLIB PROJECT --
Copyright 02.03.2010 by Bochkanov Sergey Copyright 02.03.2010 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void idwbuildmodifiedshepard(const ap::real_2d_array& xy, void idwbuildmodifiedshepard(const ap::real_2d_array& xy,
int n, int n,
int nx, int nx,
int d, int d,
int nq, int nq,
int nw, int nw,
idwinterpolant& z); idwinterpolant& z);
/************************************************************************* /*************************************************************************
IDW interpolant using modified Shepard method for non-uniform datasets.
This type of model uses constant nodal functions and interpolates using
all nodes which are closer than user-specified radius R. It may be used
when points distribution is non-uniform at the small scale, but it is at
the distances as large as R.
INPUT PARAMETERS:
XY - X and Y values, array[0..N-1,0..NX].
First NX columns contain X-values, last column contain
Y-values.
N - number of nodes, N>0.
NX - space dimension, NX>=1.
R - radius, R>0
OUTPUT PARAMETERS:
Z - IDW interpolant.
NOTES:
* if there is less than IDWKMin points within R-ball, algorithm selects
IDWKMin closest ones, so that continuity properties of interpolant are
preserved even far from points.
-- ALGLIB PROJECT --
Copyright 11.04.2010 by Bochkanov Sergey
*************************************************************************/
void idwbuildmodifiedshepardr(const ap::real_2d_array& xy,
int n,
int nx,
double r,
idwinterpolant& z);
/*************************************************************************
IDW model for noisy data. IDW model for noisy data.
This subroutine may be used to handle noisy data, i.e. data with noise in This subroutine may be used to handle noisy data, i.e. data with noise in
OUTPUT values. It differs from IDWBuildModifiedShepard() in the following OUTPUT values. It differs from IDWBuildModifiedShepard() in the following
aspects: aspects:
* nodal functions are not constrained to pass through nodes: Qi(xi)<>yi, * nodal functions are not constrained to pass through nodes: Qi(xi)<>yi,
i.e. we have fitting instead of interpolation. i.e. we have fitting instead of interpolation.
* weights which are used during least squares fitting stage are all equal * weights which are used during least squares fitting stage are all equal
to 1.0 (independently of distance) to 1.0 (independently of distance)
* "fast"-linear or constant nodal functions are not supported (either not * "fast"-linear or constant nodal functions are not supported (either not
 End of changes. 6 change blocks. 
5 lines changed or deleted 44 lines changed or added


 lsfit.h   lsfit.h 
skipping to change at line 47 skipping to change at line 47
#include "rcond.h" #include "rcond.h"
#include "matinv.h" #include "matinv.h"
#include "hblas.h" #include "hblas.h"
#include "sblas.h" #include "sblas.h"
#include "ortfac.h" #include "ortfac.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "xblas.h" #include "xblas.h"
#include "densesolver.h" #include "densesolver.h"
#include "lbfgs.h" #include "linmin.h"
#include "minlbfgs.h"
#include "minlm.h" #include "minlm.h"
/************************************************************************* /*************************************************************************
Least squares fitting report: Least squares fitting report:
TaskRCond reciprocal of task's condition number TaskRCond reciprocal of task's condition number
RMSError RMS error RMSError RMS error
AvgError average error AvgError average error
AvgRelError average relative error (for non-zero Y[I]) AvgRelError average relative error (for non-zero Y[I])
MaxError maximum error MaxError maximum error
*************************************************************************/ *************************************************************************/
skipping to change at line 75 skipping to change at line 76
}; };
struct lsfitstate struct lsfitstate
{ {
int n; int n;
int m; int m;
int k; int k;
double epsf; double epsf;
double epsx; double epsx;
int maxits; int maxits;
double stpmax;
ap::real_2d_array taskx; ap::real_2d_array taskx;
ap::real_1d_array tasky; ap::real_1d_array tasky;
ap::real_1d_array w; ap::real_1d_array w;
bool cheapfg; bool cheapfg;
bool havehess; bool havehess;
bool needf; bool needf;
bool needfg; bool needfg;
bool needfgh; bool needfgh;
int pointindex; int pointindex;
ap::real_1d_array x; ap::real_1d_array x;
ap::real_1d_array c; ap::real_1d_array c;
double f; double f;
ap::real_1d_array g; ap::real_1d_array g;
ap::real_2d_array h; ap::real_2d_array h;
int repterminationtype; int repterminationtype;
double reprmserror; double reprmserror;
double repavgerror; double repavgerror;
double repavgrelerror; double repavgrelerror;
double repmaxerror; double repmaxerror;
lmstate optstate; minlmstate optstate;
lmreport optrep; minlmreport optrep;
ap::rcommstate rstate; ap::rcommstate rstate;
}; };
/************************************************************************* /*************************************************************************
Weighted linear least squares fitting. Weighted linear least squares fitting.
QR decomposition is used to reduce task to MxM, then triangular solver or QR decomposition is used to reduce task to MxM, then triangular solver or
SVD-based solver is used depending on condition number of the system. It SVD-based solver is used depending on condition number of the system. It
allows to maximize speed and retain decent accuracy. allows to maximize speed and retain decent accuracy.
skipping to change at line 271 skipping to change at line 273
This subroutine uses only f(x[i],c) and its gradient. This subroutine uses only f(x[i],c) and its gradient.
INPUT PARAMETERS: INPUT PARAMETERS:
X - array[0..N-1,0..M-1], points (one row = one point) X - array[0..N-1,0..M-1], points (one row = one point)
Y - array[0..N-1], function values. Y - array[0..N-1], function values.
W - weights, array[0..N-1] W - weights, array[0..N-1]
C - array[0..K-1], initial approximation to the solution, C - array[0..K-1], initial approximation to the solution,
N - number of points, N>1 N - number of points, N>1
M - dimension of space M - dimension of space
K - number of parameters being fitted K - number of parameters being fitted
EpsF - stopping criterion. Algorithm stops if
|F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}
EpsX - stopping criterion. Algorithm stops if
|X(k+1)-X(k)| <= EpsX*(1+|X(k)|)
MaxIts - stopping criterion. Algorithm stops after MaxIts iterations
.
MaxIts=0 means no stopping criterion.
CheapFG - boolean flag, which is: CheapFG - boolean flag, which is:
* True if both function and gradient calculation complexit y * True if both function and gradient calculation complexit y
are less than O(M^2). An improved algorithm can are less than O(M^2). An improved algorithm can
be used which allows to save O(N*M^2) operations be used which corresponds to FGJ scheme from
per iteration with additional cost of N function/ MINLM unit.
/gradient calculations.
* False otherwise. * False otherwise.
Standard Jacibian-bases Levenberg-Marquardt algo Standard Jacibian-bases Levenberg-Marquardt algo
will be used. will be used (FJ scheme).
OUTPUT PARAMETERS: OUTPUT PARAMETERS:
State - structure which stores algorithm state between subsequent State - structure which stores algorithm state between subsequent
calls of LSFitNonlinearIteration. Used for reverse calls of LSFitNonlinearIteration. Used for reverse
communication. This structure should be passed to communication. This structure should be passed to
LSFitNonlinearIteration subroutine. LSFitNonlinearIteration subroutine.
See also: See also:
LSFitNonlinearIteration LSFitNonlinearIteration
LSFitNonlinearResults LSFitNonlinearResults
LSFitNonlinearFG (fitting without weights) LSFitNonlinearFG (fitting without weights)
LSFitNonlinearWFGH (fitting using Hessian) LSFitNonlinearWFGH (fitting using Hessian)
LSFitNonlinearFGH (fitting using Hessian, without weights) LSFitNonlinearFGH (fitting using Hessian, without weights)
NOTE
Passing EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic
stopping criterion selection (small EpsX).
-- ALGLIB -- -- ALGLIB --
Copyright 17.08.2009 by Bochkanov Sergey Copyright 17.08.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void lsfitnonlinearwfg(const ap::real_2d_array& x, void lsfitnonlinearwfg(const ap::real_2d_array& x,
const ap::real_1d_array& y, const ap::real_1d_array& y,
const ap::real_1d_array& w, const ap::real_1d_array& w,
const ap::real_1d_array& c, const ap::real_1d_array& c,
int n, int n,
int m, int m,
int k, int k,
const double& epsf,
const double& epsx,
const int& maxits,
bool cheapfg, bool cheapfg,
lsfitstate& state); lsfitstate& state);
/************************************************************************* /*************************************************************************
Nonlinear least squares fitting, no individual weights. Nonlinear least squares fitting, no individual weights.
See LSFitNonlinearWFG for more information. See LSFitNonlinearWFG for more information.
-- ALGLIB -- -- ALGLIB --
Copyright 17.08.2009 by Bochkanov Sergey Copyright 17.08.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void lsfitnonlinearfg(const ap::real_2d_array& x, void lsfitnonlinearfg(const ap::real_2d_array& x,
const ap::real_1d_array& y, const ap::real_1d_array& y,
const ap::real_1d_array& c, const ap::real_1d_array& c,
int n, int n,
int m, int m,
int k, int k,
const double& epsf,
const double& epsx,
const int& maxits,
bool cheapfg, bool cheapfg,
lsfitstate& state); lsfitstate& state);
/************************************************************************* /*************************************************************************
Weighted nonlinear least squares fitting using gradient/Hessian. Weighted nonlinear least squares fitting using gradient/Hessian.
Nonlinear task min(F(c)) is solved, where Nonlinear task min(F(c)) is solved, where
F(c) = (w[0]*(f(x[0],c)-y[0]))^2 + ... + (w[n-1]*(f(x[n-1],c)-y[n-1]))^ 2, F(c) = (w[0]*(f(x[0],c)-y[0]))^2 + ... + (w[n-1]*(f(x[n-1],c)-y[n-1]))^ 2,
skipping to change at line 369 skipping to change at line 353
-- ALGLIB -- -- ALGLIB --
Copyright 17.08.2009 by Bochkanov Sergey Copyright 17.08.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void lsfitnonlinearwfgh(const ap::real_2d_array& x, void lsfitnonlinearwfgh(const ap::real_2d_array& x,
const ap::real_1d_array& y, const ap::real_1d_array& y,
const ap::real_1d_array& w, const ap::real_1d_array& w,
const ap::real_1d_array& c, const ap::real_1d_array& c,
int n, int n,
int m, int m,
int k, int k,
const double& epsf,
const double& epsx,
const int& maxits,
lsfitstate& state); lsfitstate& state);
/************************************************************************* /*************************************************************************
Nonlinear least squares fitting using gradient/Hessian without individual Nonlinear least squares fitting using gradient/Hessian without individual
weights. See LSFitNonlinearWFGH() for more information. weights. See LSFitNonlinearWFGH() for more information.
-- ALGLIB -- -- ALGLIB --
Copyright 17.08.2009 by Bochkanov Sergey Copyright 17.08.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void lsfitnonlinearfgh(const ap::real_2d_array& x, void lsfitnonlinearfgh(const ap::real_2d_array& x,
const ap::real_1d_array& y, const ap::real_1d_array& y,
const ap::real_1d_array& c, const ap::real_1d_array& c,
int n, int n,
int m, int m,
int k, int k,
const double& epsf,
const double& epsx,
const int& maxits,
lsfitstate& state); lsfitstate& state);
/************************************************************************* /*************************************************************************
Stopping conditions for nonlinear least squares fitting.
INPUT PARAMETERS:
State - structure which stores algorithm state between calls and
which is used for reverse communication. Must be initialize
d
with LSFitNonLinearCreate???()
EpsF - stopping criterion. Algorithm stops if
|F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}
EpsX - stopping criterion. Algorithm stops if
|X(k+1)-X(k)| <= EpsX*(1+|X(k)|)
MaxIts - stopping criterion. Algorithm stops after MaxIts iterations
.
MaxIts=0 means no stopping criterion.
NOTE
Passing EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic
stopping criterion selection (according to the scheme used by MINLM unit).
-- ALGLIB --
Copyright 17.08.2009 by Bochkanov Sergey
*************************************************************************/
void lsfitnonlinearsetcond(lsfitstate& state,
double epsf,
double epsx,
int maxits);
/*************************************************************************
This function sets maximum step length
INPUT PARAMETERS:
State - structure which stores algorithm state between calls and
which is used for reverse communication. Must be
initialized with LSFitNonLinearCreate???()
StpMax - maximum step length, >=0. Set StpMax to 0.0, if you don't
want to limit step length.
Use this subroutine when you optimize target function which contains exp()
or other fast growing functions, and optimization algorithm makes too
large steps which leads to overflow. This function allows us to reject
steps that are too large (and therefore expose us to the possible
overflow) without actually calculating function value at the x+stp*d.
NOTE: non-zero StpMax leads to moderate performance degradation because
intermediate step of preconditioned L-BFGS optimization is incompatible
with limits on step size.
-- ALGLIB --
Copyright 02.04.2010 by Bochkanov Sergey
*************************************************************************/
void lsfitnonlinearsetstpmax(lsfitstate& state, double stpmax);
/*************************************************************************
Nonlinear least squares fitting. Algorithm iteration. Nonlinear least squares fitting. Algorithm iteration.
Called after inialization of the State structure with LSFitNonlinearXXX() Called after inialization of the State structure with LSFitNonlinearXXX()
subroutine. See HTML docs for examples. subroutine. See HTML docs for examples.
INPUT PARAMETERS: INPUT PARAMETERS:
State - structure which stores algorithm state between subsequent State - structure which stores algorithm state between subsequent
calls and which is used for reverse communication. Must be calls and which is used for reverse communication. Must be
initialized with LSFitNonlinearXXX() call first. initialized with LSFitNonlinearXXX() call first.
 End of changes. 12 change blocks. 
31 lines changed or deleted 62 lines changed or added


 minlm.h   minlm.h 
skipping to change at line 47 skipping to change at line 47
#include "rcond.h" #include "rcond.h"
#include "matinv.h" #include "matinv.h"
#include "hblas.h" #include "hblas.h"
#include "sblas.h" #include "sblas.h"
#include "ortfac.h" #include "ortfac.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "xblas.h" #include "xblas.h"
#include "densesolver.h" #include "densesolver.h"
#include "lbfgs.h" #include "linmin.h"
#include "minlbfgs.h"
struct lmstate struct minlmstate
{ {
bool wrongparams; bool wrongparams;
int n; int n;
int m; int m;
double epsg;
double epsf; double epsf;
double epsx; double epsx;
int maxits; int maxits;
bool xrep;
double stpmax;
int flags; int flags;
int usermode; int usermode;
ap::real_1d_array x; ap::real_1d_array x;
double f; double f;
ap::real_1d_array fi; ap::real_1d_array fi;
ap::real_2d_array j; ap::real_2d_array j;
ap::real_2d_array h; ap::real_2d_array h;
ap::real_1d_array g; ap::real_1d_array g;
bool needf; bool needf;
bool needfg; bool needfg;
bool needfgh; bool needfgh;
bool needfij; bool needfij;
bool xupdated; bool xupdated;
lbfgsstate internalstate; minlbfgsstate internalstate;
lbfgsreport internalrep; minlbfgsreport internalrep;
ap::real_1d_array xprec; ap::real_1d_array xprec;
ap::real_1d_array xbase; ap::real_1d_array xbase;
ap::real_1d_array xdir; ap::real_1d_array xdir;
ap::real_1d_array gbase; ap::real_1d_array gbase;
ap::real_1d_array xprev;
double fprev;
ap::real_2d_array rawmodel; ap::real_2d_array rawmodel;
ap::real_2d_array model; ap::real_2d_array model;
ap::real_1d_array work; ap::real_1d_array work;
ap::rcommstate rstate; ap::rcommstate rstate;
int repiterationscount; int repiterationscount;
int repterminationtype; int repterminationtype;
int repnfunc; int repnfunc;
int repnjac; int repnjac;
int repngrad; int repngrad;
int repnhess; int repnhess;
int repncholesky; int repncholesky;
int solverinfo; int solverinfo;
densesolverreport solverrep; densesolverreport solverrep;
int invinfo; int invinfo;
matinvreport invrep; matinvreport invrep;
}; };
struct lmreport struct minlmreport
{ {
int iterationscount; int iterationscount;
int terminationtype; int terminationtype;
int nfunc; int nfunc;
int njac; int njac;
int ngrad; int ngrad;
int nhess; int nhess;
int ncholesky; int ncholesky;
}; };
skipping to change at line 122 skipping to change at line 128
F = F(x[0], ..., x[n-1]) F = F(x[0], ..., x[n-1])
EXAMPLE EXAMPLE
See HTML-documentation. See HTML-documentation.
INPUT PARAMETERS: INPUT PARAMETERS:
N - dimension, N>1 N - dimension, N>1
X - initial solution, array[0..N-1] X - initial solution, array[0..N-1]
EpsF - stopping criterion. Algorithm stops if
|F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}
EpsX - stopping criterion. Algorithm stops if
|X(k+1)-X(k)| <= EpsX*(1+|X(k)|)
MaxIts - stopping criterion. Algorithm stops after MaxIts iterations
.
MaxIts=0 means no stopping criterion.
OUTPUT PARAMETERS: OUTPUT PARAMETERS:
State - structure which stores algorithm state between subsequent State - structure which stores algorithm state between subsequent
calls of MinLMIteration. Used for reverse communication. calls of MinLMIteration. Used for reverse communication.
This structure should be passed to MinLMIteration subroutin e. This structure should be passed to MinLMIteration subroutin e.
See also MinLMIteration, MinLMResults. See also MinLMIteration, MinLMResults.
NOTE NOTES:
Passing EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic 1. you may tune stopping conditions with MinLMSetCond() function
stopping criterion selection (small EpsX). 2. if target function contains exp() or other fast growing functions, and
optimization algorithm makes too large steps which leads to overflow,
use MinLMSetStpMax() function to bound algorithm's steps.
-- ALGLIB -- -- ALGLIB --
Copyright 30.03.2009 by Bochkanov Sergey Copyright 30.03.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void minlmfgh(const int& n, void minlmcreatefgh(const int& n,
const ap::real_1d_array& x, const ap::real_1d_array& x,
const double& epsf, minlmstate& state);
const double& epsx,
const int& maxits,
lmstate& state);
/************************************************************************* /*************************************************************************
LEVENBERG-MARQUARDT-LIKE METHOD FOR NON-LINEAR OPTIMIZATION LEVENBERG-MARQUARDT-LIKE METHOD FOR NON-LINEAR OPTIMIZATION
Optimization using function gradient and Jacobian. Algorithm - Levenberg- Optimization using function gradient and Jacobian. Algorithm - Levenberg-
Marquardt modification with L-BFGS pre-optimization and internal Marquardt modification with L-BFGS pre-optimization and internal
pre-conditioned L-BFGS optimization after each Levenberg-Marquardt step. pre-conditioned L-BFGS optimization after each Levenberg-Marquardt step.
Function F is represented as sum of squares: Function F is represented as sum of squares:
F = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1]) F = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1])
EXAMPLE EXAMPLE
See HTML-documentation. See HTML-documentation.
INPUT PARAMETERS: INPUT PARAMETERS:
N - dimension, N>1 N - dimension, N>1
M - number of functions f[i] M - number of functions f[i]
X - initial solution, array[0..N-1] X - initial solution, array[0..N-1]
EpsF - stopping criterion. Algorithm stops if
|F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}
EpsX - stopping criterion. Algorithm stops if
|X(k+1)-X(k)| <= EpsX*(1+|X(k)|)
MaxIts - stopping criterion. Algorithm stops after MaxIts iterations
.
MaxIts=0 means no stopping criterion.
OUTPUT PARAMETERS: OUTPUT PARAMETERS:
State - structure which stores algorithm state between subsequent State - structure which stores algorithm state between subsequent
calls of MinLMIteration. Used for reverse communication. calls of MinLMIteration. Used for reverse communication.
This structure should be passed to MinLMIteration subroutin e. This structure should be passed to MinLMIteration subroutin e.
See also MinLMIteration, MinLMResults. See also MinLMIteration, MinLMResults.
NOTE NOTES:
Passing EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic 1. you may tune stopping conditions with MinLMSetCond() function
stopping criterion selection (small EpsX). 2. if target function contains exp() or other fast growing functions, and
optimization algorithm makes too large steps which leads to overflow,
use MinLMSetStpMax() function to bound algorithm's steps.
-- ALGLIB -- -- ALGLIB --
Copyright 30.03.2009 by Bochkanov Sergey Copyright 30.03.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void minlmfgj(const int& n, void minlmcreatefgj(const int& n,
const int& m, const int& m,
const ap::real_1d_array& x, const ap::real_1d_array& x,
const double& epsf, minlmstate& state);
const double& epsx,
const int& maxits,
lmstate& state);
/************************************************************************* /*************************************************************************
CLASSIC LEVENBERG-MARQUARDT METHOD FOR NON-LINEAR OPTIMIZATION CLASSIC LEVENBERG-MARQUARDT METHOD FOR NON-LINEAR OPTIMIZATION
Optimization using Jacobi matrix. Algorithm - classic Levenberg-Marquardt Optimization using Jacobi matrix. Algorithm - classic Levenberg-Marquardt
method. method.
Function F is represented as sum of squares: Function F is represented as sum of squares:
F = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1]) F = f[0]^2(x[0],...,x[n-1]) + ... + f[m-1]^2(x[0],...,x[n-1])
EXAMPLE EXAMPLE
See HTML-documentation. See HTML-documentation.
INPUT PARAMETERS: INPUT PARAMETERS:
N - dimension, N>1 N - dimension, N>1
M - number of functions f[i] M - number of functions f[i]
X - initial solution, array[0..N-1] X - initial solution, array[0..N-1]
EpsF - stopping criterion. Algorithm stops if
|F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}
EpsX - stopping criterion. Algorithm stops if
|X(k+1)-X(k)| <= EpsX*(1+|X(k)|)
MaxIts - stopping criterion. Algorithm stops after MaxIts iterations
.
MaxIts=0 means no stopping criterion.
OUTPUT PARAMETERS: OUTPUT PARAMETERS:
State - structure which stores algorithm state between subsequent State - structure which stores algorithm state between subsequent
calls of MinLMIteration. Used for reverse communication. calls of MinLMIteration. Used for reverse communication.
This structure should be passed to MinLMIteration subroutin e. This structure should be passed to MinLMIteration subroutin e.
See also MinLMIteration, MinLMResults. See also MinLMIteration, MinLMResults.
NOTE NOTES:
Passing EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to automatic 1. you may tune stopping conditions with MinLMSetCond() function
stopping criterion selection (small EpsX). 2. if target function contains exp() or other fast growing functions, and
optimization algorithm makes too large steps which leads to overflow,
use MinLMSetStpMax() function to bound algorithm's steps.
-- ALGLIB -- -- ALGLIB --
Copyright 30.03.2009 by Bochkanov Sergey Copyright 30.03.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void minlmfj(const int& n, void minlmcreatefj(const int& n,
const int& m, const int& m,
const ap::real_1d_array& x, const ap::real_1d_array& x,
const double& epsf, minlmstate& state);
const double& epsx,
const int& maxits, /*************************************************************************
lmstate& state); This function sets stopping conditions for Levenberg-Marquardt optimization
algorithm.
INPUT PARAMETERS:
State - structure which stores algorithm state between calls and
which is used for reverse communication. Must be initialize
d
with MinLMCreate???()
EpsG - >=0
The subroutine finishes its work if the condition
||G||<EpsG is satisfied, where ||.|| means Euclidian norm,
G - gradient.
EpsF - >=0
The subroutine finishes its work if on k+1-th iteration
the condition |F(k+1)-F(k)|<=EpsF*max{|F(k)|,|F(k+1)|,1}
is satisfied.
EpsX - >=0
The subroutine finishes its work if on k+1-th iteration
the condition |X(k+1)-X(k)| <= EpsX is fulfilled.
MaxIts - maximum number of iterations. If MaxIts=0, the number of
iterations is unlimited. Only Levenberg-Marquardt
iterations are counted (L-BFGS/CG iterations are NOT
counted because their cost is very low copared to that of
LM).
Passing EpsG=0, EpsF=0, EpsX=0 and MaxIts=0 (simultaneously) will lead to
automatic stopping criterion selection (small EpsX).
-- ALGLIB --
Copyright 02.04.2010 by Bochkanov Sergey
*************************************************************************/
void minlmsetcond(minlmstate& state,
double epsg,
double epsf,
double epsx,
int maxits);
/*************************************************************************
This function turns on/off reporting.
INPUT PARAMETERS:
State - structure which stores algorithm state between calls and
which is used for reverse communication. Must be
initialized with MinLMCreate???()
NeedXRep- whether iteration reports are needed or not
Usually algorithm returns from MinLMIteration() only when it needs
function/gradient/Hessian. However, with this function we can let it stop
after each iteration (one iteration may include more than one function
evaluation), which is indicated by XUpdated field.
Both Levenberg-Marquardt and L-BFGS iterations are reported.
-- ALGLIB --
Copyright 02.04.2010 by Bochkanov Sergey
*************************************************************************/
void minlmsetxrep(minlmstate& state, bool needxrep);
/*************************************************************************
This function sets maximum step length
INPUT PARAMETERS:
State - structure which stores algorithm state between calls and
which is used for reverse communication. Must be
initialized with MinCGCreate???()
StpMax - maximum step length, >=0. Set StpMax to 0.0, if you don't
want to limit step length.
Use this subroutine when you optimize target function which contains exp()
or other fast growing functions, and optimization algorithm makes too
large steps which leads to overflow. This function allows us to reject
steps that are too large (and therefore expose us to the possible
overflow) without actually calculating function value at the x+stp*d.
NOTE: non-zero StpMax leads to moderate performance degradation because
intermediate step of preconditioned L-BFGS optimization is incompatible
with limits on step size.
-- ALGLIB --
Copyright 02.04.2010 by Bochkanov Sergey
*************************************************************************/
void minlmsetstpmax(minlmstate& state, double stpmax);
/************************************************************************* /*************************************************************************
One Levenberg-Marquardt iteration. One Levenberg-Marquardt iteration.
Called after inialization of State structure with MinLMXXX subroutine. Called after inialization of State structure with MinLMXXX subroutine.
See HTML docs for examples. See HTML docs for examples.
Input parameters: Input parameters:
State - structure which stores algorithm state between subsequent State - structure which stores algorithm state between subsequent
calls and which is used for reverse communication. Must be calls and which is used for reverse communication. Must be
skipping to change at line 270 skipping to change at line 338
If subroutine returned True, then: If subroutine returned True, then:
* if State.NeedF=True, - function value F at State.X[0..N-1] * if State.NeedF=True, - function value F at State.X[0..N-1]
is required is required
* if State.NeedFG=True - function value F and gradient G * if State.NeedFG=True - function value F and gradient G
are required are required
* if State.NeedFiJ=True - function vector f[i] and Jacobi matrix J * if State.NeedFiJ=True - function vector f[i] and Jacobi matrix J
are required are required
* if State.NeedFGH=True - function value F, gradient G and Hesian H * if State.NeedFGH=True - function value F, gradient G and Hesian H
are required are required
* if State.XUpdated=True - algorithm reports about new iteration,
State.X contains current point,
State.F contains function value.
One and only one of this fields can be set at time. One and only one of this fields can be set at time.
Results are stored: Results are stored:
* function value - in LMState.F * function value - in MinLMState.F
* gradient - in LMState.G[0..N-1] * gradient - in MinLMState.G[0..N-1]
* Jacobi matrix - in LMState.J[0..M-1,0..N-1] * Jacobi matrix - in MinLMState.J[0..M-1,0..N-1]
* Hessian - in LMState.H[0..N-1,0..N-1] * Hessian - in MinLMState.H[0..N-1,0..N-1]
-- ALGLIB -- -- ALGLIB --
Copyright 10.03.2009 by Bochkanov Sergey Copyright 10.03.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
bool minlmiteration(lmstate& state); bool minlmiteration(minlmstate& state);
/************************************************************************* /*************************************************************************
Levenberg-Marquardt algorithm results Levenberg-Marquardt algorithm results
Called after MinLMIteration returned False. Called after MinLMIteration returned False.
Input parameters: Input parameters:
State - algorithm state (used by MinLMIteration). State - algorithm state (used by MinLMIteration).
Output parameters: Output parameters:
X - array[0..N-1], solution X - array[0..N-1], solution
Rep - optimization report: Rep - optimization report:
* Rep.TerminationType completetion code: * Rep.TerminationType completetion code:
* -1 incorrect parameters were specified * -1 incorrect parameters were specified
* 1 relative function improvement is no more than * 1 relative function improvement is no more than
EpsF. EpsF.
* 2 relative step is no more than EpsX. * 2 relative step is no more than EpsX.
* 4 gradient is no more than EpsG.
* 5 MaxIts steps was taken * 5 MaxIts steps was taken
* 7 stopping conditions are too stringent,
further improvement is impossible
* Rep.IterationsCount contains iterations count * Rep.IterationsCount contains iterations count
* Rep.NFunc - number of function calculations * Rep.NFunc - number of function calculations
* Rep.NJac - number of Jacobi matrix calculations * Rep.NJac - number of Jacobi matrix calculations
* Rep.NGrad - number of gradient calculations * Rep.NGrad - number of gradient calculations
* Rep.NHess - number of Hessian calculations * Rep.NHess - number of Hessian calculations
* Rep.NCholesky - number of Cholesky decomposition calculat ions * Rep.NCholesky - number of Cholesky decomposition calculat ions
-- ALGLIB -- -- ALGLIB --
Copyright 10.03.2009 by Bochkanov Sergey Copyright 10.03.2009 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void minlmresults(const lmstate& state, ap::real_1d_array& x, lmreport& rep void minlmresults(const minlmstate& state,
); ap::real_1d_array& x,
minlmreport& rep);
#endif #endif
 End of changes. 28 change blocks. 
57 lines changed or deleted 130 lines changed or added


 mlpe.h   mlpe.h 
skipping to change at line 39 skipping to change at line 39
#include "creflections.h" #include "creflections.h"
#include "hqrnd.h" #include "hqrnd.h"
#include "matgen.h" #include "matgen.h"
#include "ablasf.h" #include "ablasf.h"
#include "ablas.h" #include "ablas.h"
#include "trfac.h" #include "trfac.h"
#include "trlinsolve.h" #include "trlinsolve.h"
#include "safesolve.h" #include "safesolve.h"
#include "rcond.h" #include "rcond.h"
#include "matinv.h" #include "matinv.h"
#include "lbfgs.h" #include "linmin.h"
#include "minlbfgs.h"
#include "hblas.h" #include "hblas.h"
#include "sblas.h" #include "sblas.h"
#include "ortfac.h" #include "ortfac.h"
#include "blas.h" #include "blas.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "xblas.h" #include "xblas.h"
#include "densesolver.h" #include "densesolver.h"
#include "mlptrain.h" #include "mlptrain.h"
 End of changes. 1 change blocks. 
1 lines changed or deleted 2 lines changed or added


 mlptrain.h   mlptrain.h 
skipping to change at line 39 skipping to change at line 39
#include "creflections.h" #include "creflections.h"
#include "hqrnd.h" #include "hqrnd.h"
#include "matgen.h" #include "matgen.h"
#include "ablasf.h" #include "ablasf.h"
#include "ablas.h" #include "ablas.h"
#include "trfac.h" #include "trfac.h"
#include "trlinsolve.h" #include "trlinsolve.h"
#include "safesolve.h" #include "safesolve.h"
#include "rcond.h" #include "rcond.h"
#include "matinv.h" #include "matinv.h"
#include "lbfgs.h" #include "linmin.h"
#include "minlbfgs.h"
#include "hblas.h" #include "hblas.h"
#include "sblas.h" #include "sblas.h"
#include "ortfac.h" #include "ortfac.h"
#include "blas.h" #include "blas.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "xblas.h" #include "xblas.h"
#include "densesolver.h" #include "densesolver.h"
 End of changes. 1 change blocks. 
1 lines changed or deleted 2 lines changed or added


 nearestneighbor.h   nearestneighbor.h 
skipping to change at line 261 skipping to change at line 261
KDT - KD-tree KDT - KD-tree
X - pre-allocated array, at least K rows, at least NX columns X - pre-allocated array, at least K rows, at least NX columns
OUTPUT PARAMETERS OUTPUT PARAMETERS
X - K rows are filled with X-values X - K rows are filled with X-values
K - number of points K - number of points
NOTE NOTE
points are ordered by distance from the query point (first = closest) points are ordered by distance from the query point (first = closest)
SEE ALSO
* KDTreeQueryResultsXY() X- and Y-values
* KDTreeQueryResultsTags() tag values
* KDTreeQueryResultsDistances() distances
-- ALGLIB -- -- ALGLIB --
Copyright 28.02.2010 by Bochkanov Sergey Copyright 28.02.2010 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void kdtreequeryresultsx(const kdtree& kdt, ap::real_2d_array& x, int& k); void kdtreequeryresultsx(const kdtree& kdt, ap::real_2d_array& x, int& k);
/************************************************************************* /*************************************************************************
X- and Y-values from last query X- and Y-values from last query
INPUT PARAMETERS INPUT PARAMETERS
KDT - KD-tree KDT - KD-tree
XY - pre-allocated array, at least K rows, at least NX+NY column s XY - pre-allocated array, at least K rows, at least NX+NY column s
OUTPUT PARAMETERS OUTPUT PARAMETERS
X - K rows are filled with points: first NX columns with X - K rows are filled with points: first NX columns with
X-values, next NY columns - with Y-values. X-values, next NY columns - with Y-values.
K - number of points K - number of points
NOTE NOTE
points are ordered by distance from the query point (first = closest) points are ordered by distance from the query point (first = closest)
SEE ALSO
* KDTreeQueryResultsX() X-values
* KDTreeQueryResultsTags() tag values
* KDTreeQueryResultsDistances() distances
-- ALGLIB -- -- ALGLIB --
Copyright 28.02.2010 by Bochkanov Sergey Copyright 28.02.2010 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void kdtreequeryresultsxy(const kdtree& kdt, ap::real_2d_array& xy, int& k) ; void kdtreequeryresultsxy(const kdtree& kdt, ap::real_2d_array& xy, int& k) ;
/************************************************************************* /*************************************************************************
point tags from last query point tags from last query
INPUT PARAMETERS INPUT PARAMETERS
KDT - KD-tree KDT - KD-tree
Tags - pre-allocated array, at least K elements Tags - pre-allocated array, at least K elements
OUTPUT PARAMETERS OUTPUT PARAMETERS
Tags - first K elements are filled with tags associated with point s, Tags - first K elements are filled with tags associated with point s,
or, when no tags were supplied, with zeros or, when no tags were supplied, with zeros
K - number of points K - number of points
NOTE NOTE
points are ordered by distance from the query point (first = closest) points are ordered by distance from the query point (first = closest)
SEE ALSO
* KDTreeQueryResultsX() X-values
* KDTreeQueryResultsXY() X- and Y-values
* KDTreeQueryResultsDistances() distances
-- ALGLIB -- -- ALGLIB --
Copyright 28.02.2010 by Bochkanov Sergey Copyright 28.02.2010 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void kdtreequeryresultstags(const kdtree& kdt, void kdtreequeryresultstags(const kdtree& kdt,
ap::integer_1d_array& tags, ap::integer_1d_array& tags,
int& k); int& k);
/************************************************************************* /*************************************************************************
Distances from last query Distances from last query
skipping to change at line 323 skipping to change at line 338
R - pre-allocated array, at least K elements R - pre-allocated array, at least K elements
OUTPUT PARAMETERS OUTPUT PARAMETERS
R - first K elements are filled with distances R - first K elements are filled with distances
(in corresponding norm) (in corresponding norm)
K - number of points K - number of points
NOTE NOTE
points are ordered by distance from the query point (first = closest) points are ordered by distance from the query point (first = closest)
SEE ALSO
* KDTreeQueryResultsX() X-values
* KDTreeQueryResultsXY() X- and Y-values
* KDTreeQueryResultsTags() tag values
-- ALGLIB -- -- ALGLIB --
Copyright 28.02.2010 by Bochkanov Sergey Copyright 28.02.2010 by Bochkanov Sergey
*************************************************************************/ *************************************************************************/
void kdtreequeryresultsdistances(const kdtree& kdt, void kdtreequeryresultsdistances(const kdtree& kdt,
ap::real_1d_array& r, ap::real_1d_array& r,
int& k); int& k);
#endif #endif
 End of changes. 4 change blocks. 
0 lines changed or deleted 20 lines changed or added


 polint.h   polint.h 
skipping to change at line 49 skipping to change at line 49
#include "rcond.h" #include "rcond.h"
#include "matinv.h" #include "matinv.h"
#include "hblas.h" #include "hblas.h"
#include "sblas.h" #include "sblas.h"
#include "ortfac.h" #include "ortfac.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "xblas.h" #include "xblas.h"
#include "densesolver.h" #include "densesolver.h"
#include "lbfgs.h" #include "linmin.h"
#include "minlbfgs.h"
#include "minlm.h" #include "minlm.h"
#include "lsfit.h" #include "lsfit.h"
#include "ratint.h" #include "ratint.h"
#include "taskgen.h" #include "taskgen.h"
/************************************************************************* /*************************************************************************
Polynomial fitting report: Polynomial fitting report:
TaskRCond reciprocal of task's condition number TaskRCond reciprocal of task's condition number
RMSError RMS error RMSError RMS error
AvgError average error AvgError average error
 End of changes. 1 change blocks. 
1 lines changed or deleted 2 lines changed or added


 ratint.h   ratint.h 
skipping to change at line 49 skipping to change at line 49
#include "rcond.h" #include "rcond.h"
#include "matinv.h" #include "matinv.h"
#include "hblas.h" #include "hblas.h"
#include "sblas.h" #include "sblas.h"
#include "ortfac.h" #include "ortfac.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "xblas.h" #include "xblas.h"
#include "densesolver.h" #include "densesolver.h"
#include "lbfgs.h" #include "linmin.h"
#include "minlbfgs.h"
#include "minlm.h" #include "minlm.h"
#include "lsfit.h" #include "lsfit.h"
/************************************************************************* /*************************************************************************
Barycentric interpolant. Barycentric interpolant.
*************************************************************************/ *************************************************************************/
struct barycentricinterpolant struct barycentricinterpolant
{ {
int n; int n;
double sy; double sy;
 End of changes. 1 change blocks. 
1 lines changed or deleted 2 lines changed or added


 spline1d.h   spline1d.h 
skipping to change at line 48 skipping to change at line 48
#include "rcond.h" #include "rcond.h"
#include "matinv.h" #include "matinv.h"
#include "hblas.h" #include "hblas.h"
#include "sblas.h" #include "sblas.h"
#include "ortfac.h" #include "ortfac.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "xblas.h" #include "xblas.h"
#include "densesolver.h" #include "densesolver.h"
#include "lbfgs.h" #include "linmin.h"
#include "minlbfgs.h"
#include "minlm.h" #include "minlm.h"
#include "lsfit.h" #include "lsfit.h"
/************************************************************************* /*************************************************************************
1-dimensional spline inteprolant 1-dimensional spline inteprolant
*************************************************************************/ *************************************************************************/
struct spline1dinterpolant struct spline1dinterpolant
{ {
int n; int n;
int k; int k;
 End of changes. 1 change blocks. 
1 lines changed or deleted 2 lines changed or added


 spline2d.h   spline2d.h 
skipping to change at line 48 skipping to change at line 48
#include "rcond.h" #include "rcond.h"
#include "matinv.h" #include "matinv.h"
#include "hblas.h" #include "hblas.h"
#include "sblas.h" #include "sblas.h"
#include "ortfac.h" #include "ortfac.h"
#include "rotations.h" #include "rotations.h"
#include "bdsvd.h" #include "bdsvd.h"
#include "svd.h" #include "svd.h"
#include "xblas.h" #include "xblas.h"
#include "densesolver.h" #include "densesolver.h"
#include "lbfgs.h" #include "linmin.h"
#include "minlbfgs.h"
#include "minlm.h" #include "minlm.h"
#include "lsfit.h" #include "lsfit.h"
#include "spline1d.h" #include "spline1d.h"
/************************************************************************* /*************************************************************************
2-dimensional spline inteprolant 2-dimensional spline inteprolant
*************************************************************************/ *************************************************************************/
struct spline2dinterpolant struct spline2dinterpolant
{ {
int k; int k;
 End of changes. 1 change blocks. 
1 lines changed or deleted 2 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/