[Scipy-svn] r3185 - in trunk/Lib/optimize: . tnc
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Jul 24 05:30:21 EDT 2007
Author: dmitrey.kroshko
Date: 2007-07-24 04:29:41 -0500 (Tue, 24 Jul 2007)
New Revision: 3185
Modified:
trunk/Lib/optimize/tnc.py
trunk/Lib/optimize/tnc/HISTORY
trunk/Lib/optimize/tnc/LICENSE
trunk/Lib/optimize/tnc/README
trunk/Lib/optimize/tnc/example.c
trunk/Lib/optimize/tnc/moduleTNC.c
trunk/Lib/optimize/tnc/tnc.c
trunk/Lib/optimize/tnc/tnc.h
Log:
tnc 1.3 connected
Modified: trunk/Lib/optimize/tnc/HISTORY
===================================================================
--- trunk/Lib/optimize/tnc/HISTORY 2007-07-23 13:01:27 UTC (rev 3184)
+++ trunk/Lib/optimize/tnc/HISTORY 2007-07-24 09:29:41 UTC (rev 3185)
@@ -1,13 +1,16 @@
# TNC Release History
-# $Jeannot: HISTORY,v 1.12 2004/04/14 21:04:42 js Exp $
+# $Jeannot: HISTORY,v 1.15 2005/01/28 18:27:31 js Exp $
+01/28/2005, V1.3 : Fix a bug in the anti-zigzaging mechanism (many thanks
+ to S.G. NASH).
+ Warning: API changes: refined stopping criterions (xtol,
+ pgtol). ftol is no more relative to the value of f.
+ new parameter offset to translate the coordinates
+04/18/2004, V1.2.5 : n==0 is now valid
04/14/2004, V1.2.4 : Fix a potential bug in the Python interface (infinity==0)
04/14/2004, V1.2.3 : Fix a bug in the Python interface (reference counting)
04/13/2004, V1.2.2 : Fix a bug in the Python interface (memory allocation)
04/13/2004, V1.2.1 : Fix a bug in the Python interface (scaling ignored)
-
-04/08/2004 Included into SciPy
-
04/03/2004, V1.2 : Memory allocation checks
04/02/2004, V1.1 : Setup script for the python module
Ability to abort the minimization at any time
Modified: trunk/Lib/optimize/tnc/LICENSE
===================================================================
--- trunk/Lib/optimize/tnc/LICENSE 2007-07-23 13:01:27 UTC (rev 3184)
+++ trunk/Lib/optimize/tnc/LICENSE 2007-07-24 09:29:41 UTC (rev 3185)
@@ -1,4 +1,4 @@
-Copyright (c) 2002-2004, Jean-Sebastien Roy (js at jeannot.org)
+Copyright (c) 2002-2005, Jean-Sebastien Roy (js at jeannot.org)
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the
Modified: trunk/Lib/optimize/tnc/README
===================================================================
--- trunk/Lib/optimize/tnc/README 2007-07-23 13:01:27 UTC (rev 3184)
+++ trunk/Lib/optimize/tnc/README 2007-07-24 09:29:41 UTC (rev 3185)
@@ -1,8 +1,8 @@
# TNC : truncated newton bound contrained minimization in C
-# Version 1.2
-# Copyright J.S. Roy (js at jeannot.org), 2002-2004
+# Version 1.3
+# Copyright J.S. Roy (js at jeannot.org), 2002-2005
# See the LICENSE file for copyright information.
-# $Jeannot: README,v 1.26 2004/04/03 16:01:48 js Exp $
+# $Jeannot: README,v 1.32 2005/01/28 15:12:09 js Exp $
This software is a C implementation of TNBC, a truncated newton minimization
package originally developed by Stephen G. Nash in Fortran.
Modified: trunk/Lib/optimize/tnc/example.c
===================================================================
--- trunk/Lib/optimize/tnc/example.c 2007-07-23 13:01:27 UTC (rev 3184)
+++ trunk/Lib/optimize/tnc/example.c 2007-07-24 09:29:41 UTC (rev 3185)
@@ -1,5 +1,5 @@
/* TNC : Minimization example */
-/* $Jeannot: example.c,v 1.17 2004/04/02 18:51:04 js Exp $ */
+/* $Jeannot: example.c,v 1.19 2005/01/28 18:27:31 js Exp $ */
#include <stdio.h>
#include <stdlib.h>
@@ -26,14 +26,16 @@
xopt[2] = {0.0, 1.0},
low[2], up[2],
eta = -1.0, stepmx = -1.0,
- accuracy = -1.0, fmin = 0.0, ftol = -1.0, rescale = -1.0;
-
+ accuracy = -1.0, fmin = 0.0, ftol = -1.0, xtol = -1.0, pgtol = -1.0,
+ rescale = -1.0;
+
low[0] = - HUGE_VAL; low[1] = 1.0;
up[0] = HUGE_VAL; up[1] = HUGE_VAL;
-
- rc = tnc(2, x, &f, g, function, NULL, low, up, NULL, TNC_MSG_ALL,
- maxCGit, maxnfeval, eta, stepmx, accuracy, fmin, ftol, rescale, &nfeval);
+ rc = tnc(2, x, &f, g, function, NULL, low, up, NULL, NULL, TNC_MSG_ALL,
+ maxCGit, maxnfeval, eta, stepmx, accuracy, fmin, ftol, xtol, pgtol,
+ rescale, &nfeval);
+
printf("After %d function evaluations, TNC returned:\n%s\n", nfeval,
tnc_rc_string[rc - TNC_MINRC]);
Modified: trunk/Lib/optimize/tnc/moduleTNC.c
===================================================================
--- trunk/Lib/optimize/tnc/moduleTNC.c 2007-07-23 13:01:27 UTC (rev 3184)
+++ trunk/Lib/optimize/tnc/moduleTNC.c 2007-07-24 09:29:41 UTC (rev 3185)
@@ -1,8 +1,8 @@
/* Python TNC module */
/*
- * Copyright (c) 2004, Jean-Sebastien Roy (js at jeannot.org)
- *
+ * Copyright (c) 2004-2005, Jean-Sebastien Roy (js at jeannot.org)
+ *
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
@@ -10,10 +10,10 @@
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
- *
+ *
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
- *
+ *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
@@ -24,7 +24,7 @@
*/
static char const rcsid[] =
- "@(#) $Jeannot: moduleTNC.c,v 1.8 2004/04/14 18:16:03 js Exp $";
+ "@(#) $Jeannot: moduleTNC.c,v 1.12 2005/01/28 18:27:31 js Exp $";
#include "Python.h"
#include <stdlib.h>
@@ -52,10 +52,9 @@
PyObject *py_float;
py_float = PyNumber_Float(py_obj);
-
- if (py_float == NULL)
- return -1;
-
+
+ if (py_float == NULL) return -1;
+
*x = PyFloat_AsDouble(py_float);
Py_DECREF(py_float);
@@ -66,15 +65,18 @@
{
int i;
double *x;
-
+
if (!PyList_Check(py_list))
+ {
+ *size = -1;
return NULL;
-
+ }
+
*size = PyList_Size(py_list);
+ if (*size <= 0) return NULL;
x = malloc((*size)*sizeof(*x));
- if (x == NULL)
- return NULL;
-
+ if (x == NULL) return NULL;
+
for (i=0; i<(*size); i++)
{
PyObject *py_float = PyList_GetItem(py_list, i);
@@ -84,30 +86,27 @@
return NULL;
}
}
-
+
return x;
}
int PyList_IntoDoubleArray(PyObject *py_list, double *x, int size)
{
int i;
-
- if (py_list == NULL)
- return 1;
-
- if (!PyList_Check(py_list))
- return 1;
-
- if (size != PyList_Size(py_list))
- return 1;
-
+
+ if (py_list == NULL) return 1;
+
+ if (!PyList_Check(py_list)) return 1;
+
+ if (size != PyList_Size(py_list)) return 1;
+
for (i=0; i<size; i++)
{
PyObject *py_float = PyList_GetItem(py_list, i);
if (py_float == NULL || PyObject_AsDouble(py_float, &(x[i])))
return 1;
}
-
+
return 0;
}
@@ -115,11 +114,10 @@
{
int i;
PyObject *py_list;
-
+
py_list = PyList_New(size);
- if (py_list == NULL)
- return NULL;
-
+ if (py_list == NULL) return NULL;
+
for (i=0; i<size; i++)
{
PyObject *py_float;
@@ -130,7 +128,7 @@
return NULL;
}
}
-
+
return py_list;
}
@@ -142,10 +140,10 @@
py_list = PyDoubleArray_AsList(py_state->n, x);
if (py_list == NULL)
{
- PyErr_SetString(PyExc_MemoryError, "Not enough memory for a list.");
+ PyErr_SetString(PyExc_MemoryError, "tnc: memory allocation failed.");
goto failure;
}
-
+
arglist = Py_BuildValue("(N)", py_list);
result = PyEval_CallObject(py_state->py_function, arglist);
Py_DECREF(arglist);
@@ -162,7 +160,7 @@
if (!PyArg_ParseTuple(result, "dO!", f, &PyList_Type, &py_grad))
{
PyErr_SetString(PyExc_ValueError,
- "Bad return value from minimized function.");
+ "tnc: invalid return value from minimized function.");
goto failure;
}
@@ -172,7 +170,7 @@
Py_DECREF(result);
return 0;
-
+
failure:
py_state->failed = 1;
Py_XDECREF(result);
@@ -181,114 +179,123 @@
PyObject *moduleTNC_minimize(PyObject *self, PyObject *args)
{
- PyObject *py_x0, *py_low, *py_up, *py_list, *py_scale;
+ PyObject *py_x0, *py_low, *py_up, *py_list, *py_scale, *py_offset;
PyObject *py_function = NULL;
pytnc_state py_state;
- int n, n1, n2, n3;
+ int n, n1, n2, n3, n4;
int rc, msg, maxCGit, maxnfeval, nfeval = 0;
- double *x, *low, *up, *scale = NULL;
- double f, eta, stepmx, accuracy, fmin, ftol, rescale;
-
- if (!PyArg_ParseTuple(args, "OO!O!O!O!iiidddddd",
+ double *x, *low, *up, *scale = NULL, *offset = NULL;
+ double f, eta, stepmx, accuracy, fmin, ftol, xtol, pgtol, rescale;
+
+ if (!PyArg_ParseTuple(args, "OO!O!O!O!O!iiidddddddd",
&py_function,
&PyList_Type, &py_x0,
&PyList_Type, &py_low,
&PyList_Type, &py_up,
&PyList_Type, &py_scale,
+ &PyList_Type, &py_offset,
&msg, &maxCGit, &maxnfeval, &eta, &stepmx, &accuracy, &fmin, &ftol,
+ &xtol, &pgtol,
&rescale
))
return NULL;
if (!PyCallable_Check(py_function))
{
- PyErr_SetString(PyExc_TypeError, "function must be callable");
+ PyErr_SetString(PyExc_TypeError, "tnc: function must be callable");
return NULL;
}
-
- if (PyList_Size(py_scale) != 0)
- {
- scale = PyList_AsDoubleArray(py_scale, &n3);
- if (scale == NULL)
- {
- PyErr_SetString(PyExc_ValueError, "Invalid parameters.");
- return NULL;
- }
+
+ scale = PyList_AsDoubleArray(py_scale, &n3);
+ if (n3 != 0 && scale == NULL)
+ {
+ PyErr_SetString(PyExc_ValueError, "tnc: invalid scaling parameters.");
+ return NULL;
}
+ offset = PyList_AsDoubleArray(py_offset, &n4);
+ if (n4 != 0 && offset == NULL)
+ {
+ PyErr_SetString(PyExc_ValueError, "tnc: invalid offset parameters.");
+ return NULL;
+ }
+
x = PyList_AsDoubleArray(py_x0, &n);
- if (x != NULL && n == 0)
+ if (n != 0 && x == NULL)
{
- free(x);
-
- PyErr_SetString(PyExc_ValueError,
- "Vector size must be greater than 0.");
+ if (scale) free(scale);
+
+ PyErr_SetString(PyExc_ValueError, "tnc: invalid initial vector.");
return NULL;
}
low = PyList_AsDoubleArray(py_low, &n1);
up = PyList_AsDoubleArray(py_up, &n2);
- if (x == NULL || low == NULL || up == NULL)
+ if ((n1 != 0 && low == NULL) || (n2 != 0 && up == NULL))
{
- if (x != NULL) free(x);
- if (low != NULL) free(low);
- if (up != NULL) free(up);
-
- PyErr_SetString(PyExc_ValueError, "Invalid parameters.");
+ if (scale) free(scale);
+ if (x) free(x);
+ if (low) free(low);
+ if (up) free(up);
+
+ PyErr_SetString(PyExc_ValueError, "tnc: invalid bounds.");
return NULL;
}
-
- if (n1 != n2 || n != n1 || (scale != NULL && n != n3))
+
+ if (n1 != n2 || n != n1 || (scale != NULL && n != n3)
+ || (offset != NULL && n != n4))
{
- free(x);
- free(low);
- free(up);
if (scale) free(scale);
-
- PyErr_SetString(PyExc_ValueError, "Vector sizes must be equal.");
+ if (offset) free(offset);
+ if (x) free(x);
+ if (low) free(low);
+ if (up) free(up);
+
+ PyErr_SetString(PyExc_ValueError, "tnc: vector sizes must be equal.");
return NULL;
}
-
+
py_state.py_function = py_function;
py_state.n = n;
py_state.failed = 0;
Py_INCREF(py_function);
-
- rc = tnc(n, x, &f, NULL, function, &py_state, low, up, scale, msg,
- maxCGit, maxnfeval, eta, stepmx, accuracy, fmin, ftol, rescale,
+
+ rc = tnc(n, x, &f, NULL, function, &py_state, low, up, scale, offset, msg,
+ maxCGit, maxnfeval, eta, stepmx, accuracy, fmin, ftol, xtol, pgtol, rescale,
&nfeval);
Py_DECREF(py_function);
-
- free(low);
- free(up);
+
+ if (low) free(low);
+ if (up) free(up);
if (scale) free(scale);
+ if (offset) free(offset);
if (py_state.failed)
{
- free(x);
+ if (x) free(x);
return NULL;
}
-
+
if (rc == TNC_ENOMEM)
{
- PyErr_SetString(PyExc_MemoryError, "Not enough memory for TNC.");
- free(x);
+ PyErr_SetString(PyExc_MemoryError, "tnc: memory allocation failed.");
+ if (x) free(x);
return NULL;
}
py_list = PyDoubleArray_AsList(n, x);
- free(x);
+ if (x) free(x);
if (py_list == NULL)
{
- PyErr_SetString(PyExc_MemoryError, "Not enough memory for a list.");
+ PyErr_SetString(PyExc_MemoryError, "tnc: memory allocation failed.");
return NULL;
}
- return Py_BuildValue("(Nii)", py_list, nfeval, rc);
+ return Py_BuildValue("(iiN)", rc, nfeval, py_list);;
}
static PyMethodDef moduleTNC_methods[] =
Modified: trunk/Lib/optimize/tnc/tnc.c
===================================================================
--- trunk/Lib/optimize/tnc/tnc.c 2007-07-23 13:01:27 UTC (rev 3184)
+++ trunk/Lib/optimize/tnc/tnc.c 2007-07-24 09:29:41 UTC (rev 3185)
@@ -1,9 +1,9 @@
-/* tnc : truncated newton bound contrained minimization
+/* tnc : truncated newton bound constrained minimization
using gradient information, in C */
/*
- * Copyright (c) 2002-2004, Jean-Sebastien Roy (js at jeannot.org)
- *
+ * Copyright (c) 2002-2005, Jean-Sebastien Roy (js at jeannot.org)
+ *
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
@@ -11,10 +11,10 @@
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
- *
+ *
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
- *
+ *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
@@ -27,12 +27,12 @@
/*
* This software is a C implementation of TNBC, a truncated newton minimization
* package originally developed by Stephen G. Nash in Fortran.
- *
+ *
* The original source code can be found at :
* http://iris.gmu.edu/~snash/nash/software/software.html
- *
+ *
* Copyright for the original TNBC fortran routines:
- *
+ *
* TRUNCATED-NEWTON METHOD: SUBROUTINES
* WRITTEN BY: STEPHEN G. NASH
* SCHOOL OF INFORMATION TECHNOLOGY & ENGINEERING
@@ -46,7 +46,7 @@
*/
static char const rcsid[] =
- "@(#) $Jeannot: tnc.c,v 1.201 2004/04/02 22:36:25 js Exp $";
+ "@(#) $Jeannot: tnc.c,v 1.205 2005/01/28 18:27:31 js Exp $";
static char const copyright[] =
"(c) 2002-2003, Jean-Sebastien Roy (js at jeannot.org)";
@@ -67,13 +67,14 @@
* Return code strings
*/
-char *tnc_rc_string[10] =
+char *tnc_rc_string[11] =
{
"Memory allocation failed",
- "Invalid parameters (less than 1 dimension)",
+ "Invalid parameters (n<0)",
"Infeasible (low bound > up bound)",
"Local minima reach (|pg| ~= 0)",
"Converged (|f_n-f_(n-1)| ~= 0)",
+ "Converged (|x_n-x_(n-1)| ~= 0)",
"Maximum number of function evaluations reached",
"Linear search failed",
"All lower bounds are equal to the upper bounds",
@@ -108,30 +109,30 @@
* Prototypes
*/
static tnc_rc tnc_minimize(int n, double x[], double *f, double g[],
- tnc_function *function, void *state,
- double xscale[], double *fscale,
+ tnc_function *function, void *state,
+ double xscale[], double xoffset[], double *fscale,
double low[], double up[], tnc_message messages,
int maxCGit, int maxnfeval, int *nfeval,
double eta, double stepmx, double accuracy,
- double fmin, double ftol, double rescale);
+ double fmin, double ftol, double xtol, double pgtol, double rescale);
-static getptc_rc getptcInit(double *reltol, double *abstol, double tnytol,
- double eta, double rmu, double xbnd,
+static getptc_rc getptcInit(double *reltol, double *abstol, double tnytol,
+ double eta, double rmu, double xbnd,
double *u, double *fu, double *gu, double *xmin,
- double *fmin, double *gmin, double *xw, double *fw,
- double *gw, double *a, double *b, double *oldf,
- double *b1, double *scxbnd, double *e, double *step,
- double *factor, logical *braktd, double *gtest1,
+ double *fmin, double *gmin, double *xw, double *fw,
+ double *gw, double *a, double *b, double *oldf,
+ double *b1, double *scxbnd, double *e, double *step,
+ double *factor, logical *braktd, double *gtest1,
double *gtest2, double *tol);
-static getptc_rc getptcIter(double big, double
- rtsmll, double *reltol, double *abstol, double tnytol,
- double fpresn, double xbnd,
+static getptc_rc getptcIter(double big, double
+ rtsmll, double *reltol, double *abstol, double tnytol,
+ double fpresn, double xbnd,
double *u, double *fu, double *gu, double *xmin,
- double *fmin, double *gmin, double *xw, double *fw,
- double *gw, double *a, double *b, double *oldf,
- double *b1, double *scxbnd, double *e, double *step,
- double *factor, logical *braktd, double *gtest1,
+ double *fmin, double *gmin, double *xw, double *fw,
+ double *gw, double *a, double *b, double *oldf,
+ double *b1, double *scxbnd, double *e, double *step,
+ double *factor, logical *braktd, double *gtest1,
double *gtest2, double *tol);
static void printCurrentIteration(int n, double f, double g[], int niter,
@@ -141,63 +142,63 @@
static ls_rc linearSearch(int n, tnc_function *function, void *state,
double low[], double up[],
- double xscale[], double fscale, int pivot[],
- double eta, double ftol, double xbnd,
- double p[], double x[], double *f,
+ double xscale[], double xoffset[], double fscale, int pivot[],
+ double eta, double ftol, double xbnd,
+ double p[], double x[], double *f,
double *alpha, double gfull[], int maxnfeval, int *nfeval);
static int tnc_direction(double *zsol, double *diagb,
double *x, double *g, int n,
- int maxCGit, int maxnfeval, int *nfeval,
+ int maxCGit, int maxnfeval, int *nfeval,
logical upd1, double yksk, double yrsr,
double *sk, double *yk, double *sr, double *yr,
- logical lreset, tnc_function *function, void *state,
- double xscale[], double fscale,
+ logical lreset, tnc_function *function, void *state,
+ double xscale[], double xoffset[], double fscale,
int *pivot, double accuracy,
double gnorm, double xnorm, double *low, double *up);
-static double stepMax(double step, int n, double x[], double p[], int pivot[],
- double low[], double up[], double xscale[]);
+static double stepMax(double step, int n, double x[], double p[], int pivot[],
+ double low[], double up[], double xscale[], double xoffset[]);
/* Active set of constraints */
-static void setContraints(int n, double x[], int pivot[], double xscale[],
- double low[], double up[]);
+static void setConstraints(int n, double x[], int pivot[], double xscale[],
+ double xoffset[], double low[], double up[]);
static logical addConstraint(int n, double x[], double p[], int pivot[],
- double low[], double up[], double xscale[]);
+ double low[], double up[], double xscale[], double xoffset[]);
-static logical removeConstraint(double gtpnew, double f,
- double *fLastConstraint, double g[], int pivot[], int n);
+static logical removeConstraint(double gtpnew, double gnorm, double pgtolfs,
+ double f, double fLastConstraint, double g[], int pivot[], int n);
static void project(int n, double x[], int pivot[]);
-static int hessianTimesVector(double v[], double gv[], int n,
- double x[], double g[], tnc_function *function, void *state,
- double xscale[], double fscale,
+static int hessianTimesVector(double v[], double gv[], int n,
+ double x[], double g[], tnc_function *function, void *state,
+ double xscale[], double xoffset[], double fscale,
double accuracy, double xnorm, double low[], double up[]);
-static int msolve(double g[], double *y, int n,
- double sk[], double yk[], double diagb[], double sr[],
- double yr[], logical upd1, double yksk, double yrsr,
+static int msolve(double g[], double *y, int n,
+ double sk[], double yk[], double diagb[], double sr[],
+ double yr[], logical upd1, double yksk, double yrsr,
logical lreset);
static void diagonalScaling(int n, double e[], double v[], double gv[],
double r[]);
static void ssbfgs(int n, double gamma, double sj[], double *hjv,
- double hjyj[], double yjsj,
+ double hjyj[], double yjsj,
double yjhyj, double vsj, double vhyj, double hjp1v[]);
-static int initPreconditioner(double diagb[], double emat[], int n,
+static int initPreconditioner(double diagb[], double emat[], int n,
logical lreset, double yksk, double yrsr,
- double sk[], double yk[], double sr[], double yr[],
+ double sk[], double yk[], double sr[], double yr[],
logical upd1);
/* Scaling */
static void coercex(int n, double x[], double low[], double up[]);
-static void unscalex(int n, double x[], double xscale[]);
+static void unscalex(int n, double x[], double xscale[], double xoffset[]);
static void scaleg(int n, double g[], double xscale[], double fscale);
-static void scalex(int n, double x[], double xscale[]);
+static void scalex(int n, double x[], double xscale[], double xoffset[]);
static void projectConstants(int n, double x[], double xscale[]);
/* Machine precision */
@@ -212,15 +213,14 @@
/* additionnal blas-like functions */
static void dneg1(int n, double v[]);
-static double dnrmi1(int n, double v[]);
/*
* This routine solves the optimization problem
- *
+ *
* minimize f(x)
* x
* subject to low <= x <= up
- *
+ *
* where x is a vector of n real variables. The method used is
* a truncated-newton algorithm (see "newton-type minimization via
* the lanczos algorithm" by s.g. nash (technical report 378, math.
@@ -230,17 +230,18 @@
* global solution), but does assume that the function is bounded below.
* it can solve problems having any number of variables, but it is
* especially useful when the number of variables (n) is large.
- *
+ *
*/
extern int tnc(int n, double x[], double *f, double g[], tnc_function *function,
- void *state, double low[], double up[], double scale[], int messages,
- int maxCGit, int maxnfeval, double eta, double stepmx,
- double accuracy, double fmin, double ftol, double rescale, int *nfeval)
+ void *state, double low[], double up[], double scale[], double offset[],
+ int messages, int maxCGit, int maxnfeval, double eta, double stepmx,
+ double accuracy, double fmin, double ftol, double xtol, double pgtol,
+ double rescale, int *nfeval)
{
int rc, frc, i, nc, nfeval_local,
free_low = TNC_FALSE, free_up = TNC_FALSE,
free_g = TNC_FALSE;
- double *xscale = NULL, fscale, epsmch, rteps;
+ double *xscale = NULL, fscale, epsmch, rteps, *xoffset = NULL;
if(nfeval==NULL)
{
@@ -248,7 +249,7 @@
nfeval = &nfeval_local;
}
*nfeval = 0;
-
+
/* Version info */
if (messages & TNC_MSG_VERS)
{
@@ -257,12 +258,18 @@
}
/* Check for errors in the input parameters */
- if (n < 1)
+ if (n == 0)
{
+ rc = TNC_CONSTANT;
+ goto cleanup;
+ }
+
+ if (n < 0)
+ {
rc = TNC_EINVAL;
goto cleanup;
}
-
+
/* Check bounds arrays */
if (low == NULL)
{
@@ -305,7 +312,7 @@
rc = TNC_MAXFUN;
goto cleanup;
}
-
+
/* Allocate g if necessary */
if(g == NULL)
{
@@ -318,7 +325,7 @@
free_g = TNC_TRUE;
}
- /* Initial function evaluation */
+ /* Initial function evaluation */
frc = function(x, f, g, state);
(*nfeval) ++;
if (frc)
@@ -326,7 +333,7 @@
rc = TNC_USERABORT;
goto cleanup;
}
-
+
/* Constant problem ? */
for (nc = 0, i = 0 ; i < n ; i++)
if ((low[i] == up[i]) || (scale != NULL && scale[i] == 0.0))
@@ -338,13 +345,19 @@
goto cleanup;
}
- /* Scaling parameters */
+ /* Scaling parameters */
xscale = malloc(sizeof(*xscale)*n);
if (xscale == NULL)
{
rc = TNC_ENOMEM;
goto cleanup;
}
+ xoffset = malloc(sizeof(*xoffset)*n);
+ if (xoffset == NULL)
+ {
+ rc = TNC_ENOMEM;
+ goto cleanup;
+ }
fscale = 1.0;
for (i = 0 ; i < n ; i++)
@@ -353,12 +366,20 @@
{
xscale[i] = fabs(scale[i]);
if (xscale[i] == 0.0)
- low[i] = up[i] = x[i];
+ xoffset[i] = low[i] = up[i] = x[i];
}
else if (low[i] != -HUGE_VAL && up[i] != HUGE_VAL)
+ {
xscale[i] = up[i] - low[i];
+ xoffset[i] = (up[i]+low[i])*0.5;
+ }
else
+ {
xscale[i] = 1.0+fabs(x[i]);
+ xoffset[i] = x[i];
+ }
+ if (offset != NULL)
+ xoffset[i] = offset[i];
}
/* Default values for parameters */
@@ -375,13 +396,16 @@
else if (maxCGit > 50) maxCGit = 50;
}
if (maxCGit > n) maxCGit = n;
- if (ftol < 0.0) ftol = 0.0;
if (accuracy <= epsmch) accuracy = rteps;
+ if (ftol < 0.0) ftol = accuracy;
+ if (pgtol < 0.0) pgtol = 1e-2 * sqrt(accuracy);
+ if (xtol < 0.0) xtol = rteps;
/* Optimisation */
rc = tnc_minimize(n, x, f, g, function, state,
- xscale, &fscale, low, up, messages,
- maxCGit, maxnfeval, nfeval, eta, stepmx, accuracy, fmin, ftol, rescale);
+ xscale, xoffset, &fscale, low, up, messages,
+ maxCGit, maxnfeval, nfeval, eta, stepmx, accuracy, fmin, ftol, xtol, pgtol,
+ rescale);
cleanup:
if (messages & TNC_MSG_EXIT)
@@ -391,7 +415,8 @@
if (free_low) free(low);
if (free_up) free(up);
if (free_g) free(g);
-
+ if (xoffset) free(xoffset);
+
return rc;
}
@@ -399,7 +424,7 @@
static void coercex(int n, double x[], double low[], double up[])
{
int i;
-
+
for (i = 0 ; i < n ; i++)
{
if (x[i]<low[i]) x[i] = low[i];
@@ -408,20 +433,20 @@
}
/* Unscale x */
-static void unscalex(int n, double x[], double xscale[])
+static void unscalex(int n, double x[], double xscale[], double xoffset[])
{
int i;
for (i = 0 ; i < n ; i++)
- x[i] *= xscale[i];
+ x[i] = x[i]*xscale[i]+xoffset[i];
}
/* Scale x */
-static void scalex(int n, double x[], double xscale[])
+static void scalex(int n, double x[], double xscale[], double xoffset[])
{
int i;
for (i = 0 ; i < n ; i++)
if (xscale[i]>0.0)
- x[i] /= xscale[i];
+ x[i] = (x[i]-xoffset[i])/xscale[i];
}
/* Scale g */
@@ -433,18 +458,16 @@
}
/* Caculate the pivot vector */
-static void setContraints(int n, double x[], int pivot[], double xscale[],
- double low[], double up[])
+static void setConstraints(int n, double x[], int pivot[], double xscale[],
+ double xoffset[], double low[], double up[])
{
int i;
double epsmch;
-
+
epsmch = mchpr1();
for (i = 0; i < n; i++)
{
- double tol;
-
/* tolerances should be better ajusted */
if (xscale[i] == 0.0)
{
@@ -452,13 +475,13 @@
}
else
{
- tol = epsmch * 10.0 * (fabs(low[i]) + 1.0);
- if ((x[i]*xscale[i] - low[i] <= tol) && low[i] != - HUGE_VAL)
+ if (low[i] != - HUGE_VAL &&
+ (x[i]*xscale[i]+xoffset[i] - low[i] <= epsmch * 10.0 * (fabs(low[i]) + 1.0)))
pivot[i] = -1;
else
{
- tol = epsmch * 10.0 * (fabs(up[i]) + 1.0);
- if ((x[i]*xscale[i] - up[i] >= tol) && up[i] != HUGE_VAL)
+ if (up[i] != HUGE_VAL &&
+ (x[i]*xscale[i]+xoffset[i] - up[i] >= epsmch * 10.0 * (fabs(up[i]) + 1.0)))
pivot[i] = 1;
else
pivot[i] = 0;
@@ -474,15 +497,16 @@
* in this routine) with a further diagonal scaling
* (see routine diagonalscaling).
*/
-static tnc_rc tnc_minimize(int n, double x[],
- double *f, double gfull[], tnc_function *function, void *state,
- double xscale[], double *fscale,
- double low[], double up[], tnc_message messages,
- int maxCGit, int maxnfeval, int *nfeval, double eta, double stepmx,
- double accuracy, double fmin, double ftol, double rescale)
+static tnc_rc tnc_minimize(int n, double x[],
+ double *f, double gfull[], tnc_function *function, void *state,
+ double xscale[], double xoffset[], double *fscale,
+ double low[], double up[], tnc_message messages,
+ int maxCGit, int maxnfeval, int *nfeval, double eta, double stepmx,
+ double accuracy, double fmin, double ftol, double xtol, double pgtol,
+ double rescale)
{
double fLastReset, difnew, epsmch, epsred, oldgtp,
- difold, oldf, rteps, xnorm, newscale,
+ difold, oldf, xnorm, newscale,
gnorm, ustpmax, fLastConstraint, spe, yrsr, yksk,
*temp = NULL, *sk = NULL, *yk = NULL, *diagb = NULL, *sr = NULL,
*yr = NULL, *oldg = NULL, *pk = NULL, *g = NULL;
@@ -491,7 +515,7 @@
logical lreset, newcon, upd1, remcon;
tnc_rc rc = TNC_ENOMEM; /* Default error */
- /* Allocate temporary vectors */
+ /* Allocate temporary vectors */
oldg = malloc(sizeof(*oldg)*n);
if (oldg == NULL) goto cleanup;
g = malloc(sizeof(*g)*n);
@@ -517,25 +541,24 @@
/* Initialize variables */
epsmch = mchpr1();
- rteps = sqrt(epsmch);
difnew = 0.0;
epsred = 0.05;
upd1 = TNC_TRUE;
icycle = n - 1;
newcon = TNC_TRUE;
-
+
/* Uneeded initialisations */
lreset = TNC_FALSE;
yrsr = 0.0;
yksk = 0.0;
-
+
/* Initial scaling */
- scalex(n, x, xscale);
+ scalex(n, x, xscale, xoffset);
(*f) *= *fscale;
/* initial pivot calculation */
- setContraints(n, x, pivot, xscale, low, up);
+ setConstraints(n, x, pivot, xscale, xoffset, low, up);
dcopy1(n, gfull, g);
scaleg(n, g, xscale, *fscale);
@@ -564,14 +587,14 @@
/* Start of main iterative loop */
while(TNC_TRUE)
{
- /* Tolerance should be user modifiable */
- if (dnrmi1(n, g) <= 1.0e-2*rteps*fabs(*f))
+ /* Local minimum test */
+ if (dnrm21(n, g) <= pgtol * (*fscale))
{
/* |PG| == 0.0 => local minimum */
dcopy1(n, gfull, g);
project(n, g, pivot);
if (messages & TNC_MSG_INFO) fprintf(stderr,
- "tnc: |pg| = %g -> local minimum\n",dnrmi1(n, g));
+ "tnc: |pg| = %g -> local minimum\n", dnrm21(n, g) / (*fscale));
rc = TNC_LOCALMINIMUM;
break;
}
@@ -584,11 +607,11 @@
}
/* Rescale function if necessary */
- newscale = dnrmi1(n, g);
+ newscale = dnrm21(n, g);
if ((newscale > epsmch) && (fabs(log10(newscale)) > rescale))
{
newscale = 1.0/newscale;
-
+
*f *= newscale;
*fscale *= newscale;
gnorm *= newscale;
@@ -603,7 +626,7 @@
icycle = n - 1;
newcon = TNC_TRUE;
- if (messages & TNC_MSG_INFO) fprintf(stderr,
+ if (messages & TNC_MSG_INFO) fprintf(stderr,
"tnc: fscale = %g\n", *fscale);
}
@@ -615,7 +638,7 @@
/* Compute the new search direction */
frc = tnc_direction(pk, diagb, x, g, n, maxCGit, maxnfeval, nfeval,
upd1, yksk, yrsr, sk, yk, sr, yr,
- lreset, function, state, xscale, *fscale,
+ lreset, function, state, xscale, xoffset, *fscale,
pivot, accuracy, gnorm, xnorm, low, up);
if (frc == -1)
@@ -658,8 +681,8 @@
ustpmax = stepmx / (dnrm21(n, pk) + epsmch);
/* Maximum constrained step length */
- spe = stepMax(ustpmax, n, x, pk, pivot, low, up, xscale);
-
+ spe = stepMax(ustpmax, n, x, pk, pivot, low, up, xscale, xoffset);
+
if (spe > 0.0)
{
ls_rc lsrc;
@@ -668,7 +691,7 @@
/* Perform the linear search */
lsrc = linearSearch(n, function, state, low, up,
- xscale, *fscale, pivot,
+ xscale, xoffset, *fscale, pivot,
eta, ftol, spe, pk, x, f, &alpha, gfull, maxnfeval, nfeval);
if (lsrc == LS_ENOMEM)
@@ -683,6 +706,12 @@
break;
}
+ if (lsrc == LS_FAIL)
+ {
+ rc = TNC_LSFAIL;
+ break;
+ }
+
/* If we went up to the maximum unconstrained step, increase it */
if (alpha >= 0.9 * ustpmax)
{
@@ -717,7 +746,7 @@
if (newcon)
{
- if(!addConstraint(n, x, pk, pivot, low, up, xscale))
+ if(!addConstraint(n, x, pk, pivot, low, up, xscale, xoffset))
{
if(*nfeval == oldnfeval)
{
@@ -725,6 +754,7 @@
break;
}
}
+
fLastConstraint = *f;
}
@@ -750,18 +780,36 @@
gnorm = dnrm21(n, temp);
/* Reset pivot */
- remcon = removeConstraint(oldgtp, *f, &fLastConstraint, g, pivot, n);
+ remcon = removeConstraint(oldgtp, gnorm, pgtol * (*fscale), *f,
+ fLastConstraint, g, pivot, n);
+ /* If a constraint is removed */
+ if (remcon)
+ {
+ /* Recalculate gnorm and reset fLastConstraint */
+ dcopy1(n, g, temp);
+ project(n, temp, pivot);
+ gnorm = dnrm21(n, temp);
+ fLastConstraint = *f;
+ }
+
if (!remcon && !newcon)
{
- /* No constraint removed & no new constraint : test for convergence */
- if (fabs(difnew) <= ftol*epsmch*0.5*(fabs(oldf)+fabs(*f)))
+ /* No constraint removed & no new constraint : tests for convergence */
+ if (fabs(difnew) <= ftol * (*fscale))
{
- if (messages & TNC_MSG_INFO) fprintf(stderr,
- "tnc: |fn-fn-1] = %g -> convergence\n",fabs(difnew));
- rc = TNC_CONVERGED;
+ if (messages & TNC_MSG_INFO) fprintf(stderr,
+ "tnc: |fn-fn-1] = %g -> convergence\n", fabs(difnew) / (*fscale));
+ rc = TNC_FCONVERGED;
break;
}
+ if (alpha * dnrm21(n, pk) <= xtol)
+ {
+ if (messages & TNC_MSG_INFO) fprintf(stderr,
+ "tnc: |xn-xn-1] = %g -> convergence\n", alpha * dnrm21(n, pk));
+ rc = TNC_XCONVERGED;
+ break;
+ }
}
project(n, g, pivot);
@@ -781,7 +829,7 @@
/* Set up parameters used in updating the preconditioning strategy */
yksk = ddot1(n, yk, sk);
-
+
if (icycle == (n - 1) || difnew < epsred * (fLastReset - *f))
lreset = TNC_TRUE;
else
@@ -798,11 +846,11 @@
niter, *nfeval, pivot);
/* Unscaling */
- unscalex(n, x, xscale);
+ unscalex(n, x, xscale, xoffset);
coercex(n, x, low, up);
(*f) /= *fscale;
-cleanup:
+cleanup:
if (oldg) free(oldg);
if (g) free(g);
if (temp) free(temp);
@@ -813,7 +861,7 @@
if (yk) free(yk);
if (sr) free(sr);
if (yr) free(yr);
-
+
if (pivot) free(pivot);
return rc;
@@ -835,7 +883,7 @@
}
/*
- * Set x[i] = 0.0 if direction i is currently constrained
+ * Set x[i] = 0.0 if direction i is currently constrained
*/
static void project(int n, double x[], int pivot[])
{
@@ -860,7 +908,7 @@
* Compute the maximum allowable step length
*/
static double stepMax(double step, int n, double x[], double dir[],
- int pivot[], double low[], double up[], double xscale[])
+ int pivot[], double low[], double up[], double xscale[], double xoffset[])
{
int i;
double t;
@@ -872,17 +920,17 @@
{
if (dir[i] < 0.0)
{
- t = low[i]/xscale[i] - x[i];
+ t = (low[i]-xoffset[i])/xscale[i] - x[i];
if (t > step * dir[i]) step = t / dir[i];
}
else
{
- t = up[i]/xscale[i] - x[i];
+ t = (up[i]-xoffset[i])/xscale[i] - x[i];
if (t < step * dir[i]) step = t / dir[i];
}
}
}
-
+
return step;
}
@@ -890,7 +938,7 @@
* Update the constraint vector pivot if a new constraint is encountered
*/
static logical addConstraint(int n, double x[], double p[], int pivot[],
- double low[], double up[], double xscale[])
+ double low[], double up[], double xscale[], double xoffset[])
{
int i, newcon = TNC_FALSE;
double tol, epsmch;
@@ -904,20 +952,20 @@
if (p[i] < 0.0 && low[i] != - HUGE_VAL)
{
tol = epsmch * 10.0 * (fabs(low[i]) + 1.0);
- if (x[i]*xscale[i] - low[i] <= tol)
+ if (x[i]*xscale[i]+xoffset[i] - low[i] <= tol)
{
pivot[i] = -1;
- x[i] = low[i]/xscale[i];
+ x[i] = (low[i]-xoffset[i])/xscale[i];
newcon = TNC_TRUE;
}
}
else if (up[i] != HUGE_VAL)
{
tol = epsmch * 10.0 * (fabs(up[i]) + 1.0);
- if (up[i] - x[i]*xscale[i] <= tol)
+ if (up[i] - (x[i]*xscale[i]+xoffset[i]) <= tol)
{
pivot[i] = 1;
- x[i] = up[i]/xscale[i];
+ x[i] = (up[i]-xoffset[i])/xscale[i];
newcon = TNC_TRUE;
}
}
@@ -929,36 +977,33 @@
/*
* Check if a constraint is no more active
*/
-static logical removeConstraint(double gtpnew, double f,
- double *fLastConstraint, double g[], int pivot[], int n)
+static logical removeConstraint(double gtpnew, double gnorm, double pgtolfs,
+ double f, double fLastConstraint, double g[], int pivot[], int n)
{
double cmax, t;
int imax, i;
- logical ltest;
+ if (((fLastConstraint - f) <= (gtpnew * -0.5)) && (gnorm > pgtolfs))
+ return TNC_FALSE;
+
imax = -1;
cmax = 0.0;
- ltest = (*fLastConstraint - f) <= (gtpnew * -0.5);
+
for (i = 0; i < n; i++)
{
- if (pivot[i] != 2)
+ if (pivot[i] == 2)
+ continue;
+ t = -pivot[i] * g[i];
+ if (t < cmax)
{
- t = -pivot[i] * g[i];
- if (t < 0.0)
- {
- if ((!ltest) && (cmax > t))
- {
- cmax = t;
- imax = i;
- }
- }
+ cmax = t;
+ imax = i;
}
}
if (imax != -1)
{
pivot[imax] = 0;
- *fLastConstraint = f;
return TNC_TRUE;
}
else
@@ -981,11 +1026,11 @@
*/
static int tnc_direction(double *zsol, double *diagb,
double *x, double g[], int n,
- int maxCGit, int maxnfeval, int *nfeval,
+ int maxCGit, int maxnfeval, int *nfeval,
logical upd1, double yksk, double yrsr,
double *sk, double *yk, double *sr, double *yr,
- logical lreset, tnc_function *function, void *state,
- double xscale[], double fscale,
+ logical lreset, tnc_function *function, void *state,
+ double xscale[], double xoffset[], double fscale,
int *pivot, double accuracy,
double gnorm, double xnorm, double low[], double up[])
{
@@ -1063,7 +1108,7 @@
project(n, v, pivot);
frc = hessianTimesVector(v, gv, n, x, g, function, state,
- xscale, fscale, accuracy, xnorm, low, up);
+ xscale, xoffset, fscale, accuracy, xnorm, low, up);
++(*nfeval);
if (frc) goto cleanup;
project(n, gv, pivot);
@@ -1123,7 +1168,7 @@
return frc;
}
-/*
+/*
* Update the preconditioning matrix based on a diagonal version
* of the bfgs quasi-newton update.
*/
@@ -1161,14 +1206,14 @@
/*
* Hessian vector product through finite differences
*/
-static int hessianTimesVector(double v[], double gv[], int n,
- double x[], double g[], tnc_function *function, void *state,
- double xscale[], double fscale,
+static int hessianTimesVector(double v[], double gv[], int n,
+ double x[], double g[], tnc_function *function, void *state,
+ double xscale[], double xoffset[], double fscale,
double accuracy, double xnorm, double low[], double up[])
{
double dinv, f, delta, *xv;
int i, frc;
-
+
xv = malloc(sizeof(*xv)*n);
if (xv == NULL) return -1;
@@ -1176,7 +1221,7 @@
for (i = 0; i < n; i++)
xv[i] = x[i] + delta * v[i];
- unscalex(n, xv, xscale);
+ unscalex(n, xv, xscale, xoffset);
coercex(n, xv, low, up);
frc = function(xv, &f, gv, state);
free(xv);
@@ -1186,22 +1231,22 @@
dinv = 1.0 / delta;
for (i = 0; i < n; i++)
gv[i] = (gv[i] - g[i]) * dinv;
-
+
projectConstants(n, gv, xscale);
return 0;
}
/*
- * This routine acts as a preconditioning step for the
- * linear conjugate-gradient routine. It is also the
- * method of computing the search direction from the
- * gradient for the non-linear conjugate-gradient code.
- * It represents a two-step self-scaled bfgs formula.
+ * This routine acts as a preconditioning step for the
+ * linear conjugate-gradient routine. It is also the
+ * method of computing the search direction from the
+ * gradient for the non-linear conjugate-gradient code.
+ * It represents a two-step self-scaled bfgs formula.
*/
-static int msolve(double g[], double y[], int n,
- double sk[], double yk[], double diagb[], double sr[],
- double yr[], logical upd1, double yksk, double yrsr,
+static int msolve(double g[], double y[], int n,
+ double sk[], double yk[], double diagb[], double sr[],
+ double yr[], logical upd1, double yksk, double yrsr,
logical lreset)
{
double ghyk, ghyr, yksr, ykhyk, ykhyr, yrhyr, rdiagb, gsr, gsk;
@@ -1257,12 +1302,12 @@
ghyk = ddot1(n, hyk, g);
ssbfgs(n, 1.0, sk, hg, hyk, yksk, ykhyk, gsk, ghyk, y);
}
-
+
cleanup:
if (hg) free(hg);
if (hyk) free(hyk);
if (hyr) free(hyr);
-
+
return frc;
}
@@ -1270,7 +1315,7 @@
* Self-scaled BFGS
*/
static void ssbfgs(int n, double gamma, double sj[], double hjv[],
- double hjyj[], double yjsj,
+ double hjyj[], double yjsj,
double yjhyj, double vsj, double vhyj, double hjp1v[])
{
double beta, delta;
@@ -1294,9 +1339,9 @@
/*
* Initialize the preconditioner
*/
-static int initPreconditioner(double diagb[], double emat[], int n,
+static int initPreconditioner(double diagb[], double emat[], int n,
logical lreset, double yksk, double yrsr,
- double sk[], double yk[], double sr[], double yr[],
+ double sk[], double yk[], double sr[], double yr[],
logical upd1)
{
double srds, yrsk, td, sds;
@@ -1308,11 +1353,11 @@
dcopy1(n, diagb, emat);
return 0;
}
-
+
bsk = malloc(sizeof(*bsk)*n);
if (bsk == NULL) return -1;
-
- if (lreset)
+
+ if (lreset)
{
for (i = 0; i < n; i++) bsk[i] = diagb[i] * sk[i];
sds = ddot1(n, sk, bsk);
@@ -1344,7 +1389,7 @@
for (i = 0; i < n; i++)
emat[i] = emat[i] - bsk[i] * bsk[i] / sds + yk[i] * yk[i] / yksk;
}
-
+
free(bsk);
return 0;
}
@@ -1355,7 +1400,7 @@
*/
static ls_rc linearSearch(int n, tnc_function *function, void *state,
double low[], double up[],
- double xscale[], double fscale, int pivot[],
+ double xscale[], double xoffset[], double fscale, int pivot[],
double eta, double ftol, double xbnd,
double p[], double x[], double *f,
double *alpha, double gfull[], int maxnfeval, int *nfeval)
@@ -1368,7 +1413,7 @@
ls_rc rc;
getptc_rc itest;
logical braktd;
-
+
rc = LS_ENOMEM;
temp = malloc(sizeof(*temp)*n);
if (temp == NULL) goto cleanup;
@@ -1401,7 +1446,7 @@
itcnt = 0;
/* Set the estimated relative precision in f(x). */
- fpresn = epsmch * ftol;
+ fpresn = ftol;
u = *alpha;
fu = *f;
@@ -1409,11 +1454,11 @@
rmu = 1e-4;
/* Setup */
- itest = getptcInit(&reltol, &abstol, tnytol, eta, rmu,
+ itest = getptcInit(&reltol, &abstol, tnytol, eta, rmu,
xbnd, &u, &fu, &gu, alpha, &fmin, &gmin, &xw, &fw, &gw, &a, &b,
&oldf, &b1, &scxbnd, &e, &step, &factor, &braktd, >est1, >est2, &tol);
- /* If itest == GETPTC_EVAL, the algorithm requires the function value to be
+ /* If itest == GETPTC_EVAL, the algorithm requires the function value to be
calculated */
while(itest == GETPTC_EVAL)
{
@@ -1425,7 +1470,7 @@
temp[i] = x[i] + ualpha * p[i];
/* Function evaluation */
- unscalex(n, temp, xscale);
+ unscalex(n, temp, xscale, xoffset);
coercex(n, temp, low, up);
frc = function(temp, &fu, tempgfull, state);
@@ -1445,7 +1490,7 @@
itest = getptcIter(big, rtsmll, &reltol, &abstol, tnytol, fpresn,
xbnd, &u, &fu, &gu, alpha, &fmin, &gmin, &xw, &fw, &gw, &a, &b,
&oldf, &b1, &scxbnd, &e, &step, &factor, &braktd, >est1, >est2, &tol);
-
+
/* New best point ? */
if (*alpha == ualpha)
dcopy1(n, tempgfull, newgfull);
@@ -1481,17 +1526,18 @@
* in which a lower point is to be found and from this getptc computes a
* point at which the function can be evaluated by the calling program.
*/
-static getptc_rc getptcInit(double *reltol, double *abstol, double tnytol,
- double eta, double rmu, double xbnd,
+static getptc_rc getptcInit(double *reltol, double *abstol, double tnytol,
+ double eta, double rmu, double xbnd,
double *u, double *fu, double *gu, double *xmin,
- double *fmin, double *gmin, double *xw, double *fw,
- double *gw, double *a, double *b, double *oldf,
- double *b1, double *scxbnd, double *e, double *step,
- double *factor, logical *braktd, double *gtest1,
+ double *fmin, double *gmin, double *xw, double *fw,
+ double *gw, double *a, double *b, double *oldf,
+ double *b1, double *scxbnd, double *e, double *step,
+ double *factor, logical *braktd, double *gtest1,
double *gtest2, double *tol)
{
/* Check input parameters */
- if (*u <= 0.0 || xbnd <= tnytol || *gu > 0.0) return GETPTC_EINVAL;
+ if (*u <= 0.0 || xbnd <= tnytol || *gu > 0.0)
+ return GETPTC_EINVAL;
if (xbnd < *abstol) *abstol = xbnd;
*tol = *abstol;
@@ -1530,7 +1576,7 @@
/* If the step is too large, replace by the scaled bound (so as to */
/* compute the new point on the boundary). */
if (*step >= *scxbnd)
- {
+ {
*step = *scxbnd;
/* Move sxbd to the left so that sbnd + tol(xbnd) = xbnd. */
*scxbnd -= (*reltol * fabs(xbnd) + *abstol) / (1.0 + *reltol);
@@ -1541,17 +1587,17 @@
return GETPTC_EVAL;
}
-static getptc_rc getptcIter(double big, double
- rtsmll, double *reltol, double *abstol, double tnytol,
+static getptc_rc getptcIter(double big, double
+ rtsmll, double *reltol, double *abstol, double tnytol,
double fpresn, double xbnd,
double *u, double *fu, double *gu, double *xmin,
- double *fmin, double *gmin, double *xw, double *fw,
- double *gw, double *a, double *b, double *oldf,
- double *b1, double *scxbnd, double *e, double *step,
- double *factor, logical *braktd, double *gtest1,
+ double *fmin, double *gmin, double *xw, double *fw,
+ double *gw, double *a, double *b, double *oldf,
+ double *b1, double *scxbnd, double *e, double *step,
+ double *factor, logical *braktd, double *gtest1,
double *gtest2, double *tol)
{
- double abgw, absr, p, q, r, s, scale, denom,
+ double abgw, absr, p, q, r, s, scale, denom,
a1, d1, d2, sumsq, abgmin, chordm, chordu,
xmidpt, twotol;
logical convrg;
@@ -1625,7 +1671,7 @@
xmidpt = 0.5 * (*a + *b);
/* Check termination criteria */
- convrg = (fabs(xmidpt) <= twotol - 0.5 * (*b - *a)) ||
+ convrg = (fabs(xmidpt) <= twotol - 0.5 * (*b - *a)) ||
(fabs(*gmin) <= *gtest2 && *fmin < *oldf && ((fabs(*xmin - xbnd) > *tol) ||
(! (*braktd))));
if (convrg)
@@ -1638,7 +1684,7 @@
* unimodality constant, tol. If the change in f(x) is larger than
* expected, reduce the value of tol.
*/
- if (fabs(*oldf - *fw) <= fpresn * 0.5 * (fabs(*fw) + fabs(*oldf)))
+ if (fabs(*oldf - *fw) <= fpresn)
return GETPTC_FAIL;
*tol = 0.1 * *tol;
if (*tol < tnytol) return GETPTC_FAIL;
@@ -1664,7 +1710,7 @@
abgw = fabs(*gw);
abgmin = fabs(*gmin);
s = sqrt(abgmin) * sqrt(abgw);
- if (*gw / abgw * *gmin > 0.0)
+ if (*gw / abgw * *gmin > 0.0)
{
if (r >= s || r <= -s)
{
@@ -1777,7 +1823,7 @@
/* If the step is too large, replace by the scaled bound (so as to */
/* compute the new point on the boundary). */
if (*step >= *scxbnd)
- {
+ {
*step = *scxbnd;
/* Move sxbd to the left so that sbnd + tol(xbnd) = xbnd. */
*scxbnd -= (*reltol * fabs(xbnd) + *abstol) / (1.0 + *reltol);
@@ -1795,7 +1841,7 @@
static double mchpr1(void)
{
static double epsmch = 0.0;
-
+
if (epsmch == 0.0)
{
double eps = 1.0;
@@ -1852,16 +1898,6 @@
return dtemp;
}
-/* Infinity norm */
-static double dnrmi1(int n, double v[])
-{
- int i;
- double dtemp, dmax;
- for (dmax = fabs(v[0]), i = 1; i < n; i++)
- if ((dtemp = fabs(v[i])) > dmax) dmax = dtemp;
- return dmax;
-}
-
/* Euclidian norm */
static double dnrm21(int n, double dx[])
{
Modified: trunk/Lib/optimize/tnc/tnc.h
===================================================================
--- trunk/Lib/optimize/tnc/tnc.h 2007-07-23 13:01:27 UTC (rev 3184)
+++ trunk/Lib/optimize/tnc/tnc.h 2007-07-24 09:29:41 UTC (rev 3185)
@@ -1,9 +1,9 @@
-/* tnc : truncated newton bound contrained minimization
+/* tnc : truncated newton bound constrained minimization
using gradient information, in C */
/*
- * Copyright (c) 2002-2004, Jean-Sebastien Roy (js at jeannot.org)
- *
+ * Copyright (c) 2002-2005, Jean-Sebastien Roy (js at jeannot.org)
+ *
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
@@ -11,10 +11,10 @@
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
- *
+ *
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
- *
+ *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
@@ -27,12 +27,12 @@
/*
* This software is a C implementation of TNBC, a truncated newton minimization
* package originally developed by Stephen G. Nash in Fortran.
- *
+ *
* The original source code can be found at :
* http://iris.gmu.edu/~snash/nash/software/software.html
- *
+ *
* Copyright for the original TNBC fortran routines:
- *
+ *
* TRUNCATED-NEWTON METHOD: SUBROUTINES
* WRITTEN BY: STEPHEN G. NASH
* SCHOOL OF INFORMATION TECHNOLOGY & ENGINEERING
@@ -40,12 +40,12 @@
* FAIRFAX, VA 22030
*/
-/* $Jeannot: tnc.h,v 1.48 2004/04/03 16:01:49 js Exp $ */
+/* $Jeannot: tnc.h,v 1.55 2005/01/28 18:27:31 js Exp $ */
#ifndef _TNC_
#define _TNC_
-#define TNC_VERSION "1.2"
+#define TNC_VERSION "1.3"
#ifdef __cplusplus
extern "C" {
@@ -71,16 +71,17 @@
typedef enum
{
TNC_MINRC = -3, /* Constant to add to get the rc_string */
- TNC_ENOMEM = -3, /* Memory allocation failed */
- TNC_EINVAL = -2, /* Invalid parameters (les than 1 dimension) */
+ TNC_ENOMEM = -3, /* Memory allocation failed */
+ TNC_EINVAL = -2, /* Invalid parameters (n<0) */
TNC_INFEASIBLE = -1, /* Infeasible (low bound > up bound) */
TNC_LOCALMINIMUM = 0, /* Local minima reach (|pg| ~= 0) */
- TNC_CONVERGED = 1, /* Converged (|f_n-f_(n-1)| ~= 0) */
- TNC_MAXFUN = 2, /* Max. number of function evaluations reach */
- TNC_LSFAIL = 3, /* Linear search failed */
- TNC_CONSTANT = 4, /* All lower bounds are equal to the upper bounds */
- TNC_NOPROGRESS = 5, /* Unable to progress */
- TNC_USERABORT = 6 /* User requested end of minization */
+ TNC_FCONVERGED = 1, /* Converged (|f_n-f_(n-1)| ~= 0) */
+ TNC_XCONVERGED = 2, /* Converged (|x_n-x_(n-1)| ~= 0) */
+ TNC_MAXFUN = 3, /* Max. number of function evaluations reach */
+ TNC_LSFAIL = 4, /* Linear search failed */
+ TNC_CONSTANT = 5, /* All lower bounds are equal to the upper bounds */
+ TNC_NOPROGRESS = 6, /* Unable to progress */
+ TNC_USERABORT = 7 /* User requested end of minization */
} tnc_rc;
/*
@@ -89,7 +90,7 @@
* return code rc.
*/
-extern char *tnc_rc_string[10];
+extern char *tnc_rc_string[11];
/*
* A function as required by tnc
@@ -109,7 +110,7 @@
* tnc : minimize a function with variables subject to bounds, using
* gradient information.
*
- * n : number of variables (must be > 0)
+ * n : number of variables (must be >= 0)
* x : on input, initial estimate ; on output, the solution
* f : on output, the function value at the solution
* g : on output, the gradient value at the solution
@@ -124,7 +125,10 @@
* if up == NULL, the upper bounds are removed.
* scale : scaling factors to apply to each variable
* if NULL, the factors are up-low for interval bounded variables
- * and 1+|x] fo the others.
+ * and 1+|x] for the others.
+ * offset : constant to substract to each variable
+ * if NULL, the constant are (up+low)/2 for interval bounded
+ * variables and x for the others.
* messages : see the tnc_message enum
* maxCGit : max. number of hessian*vector evaluation per main iteration
* if maxCGit == 0, the direction chosen is -gradient
@@ -137,9 +141,15 @@
* if <= machine_precision, set to sqrt(machine_precision)
* fmin : minimum function value estimate
* ftol : precision goal for the value of f in the stoping criterion
- * relative to the machine precision and the value of f.
- * if ftol < 0.0, ftol is set to 0.0
- * rescale : Scaling factor (in log10) used to trigger rescaling
+ * if ftol < 0.0, ftol is set to accuracy
+ * xtol : precision goal for the value of x in the stopping criterion
+ * (after applying x scaling factors)
+ * if xtol < 0.0, xtol is set to sqrt(machine_precision)
+ * pgtol : precision goal for the value of the projected gradient in the
+ * stopping criterion (after applying x scaling factors)
+ * if pgtol < 0.0, pgtol is set to 1e-2 * sqrt(accuracy)
+ * setting it to 0.0 is not recommended
+ * rescale : f scaling factor (in log10) used to trigger f value rescaling
* if 0, rescale at each iteration
* if a big value, never rescale
* if < 0, rescale is set to 1.3
@@ -152,9 +162,10 @@
*/
extern int tnc(int n, double x[], double *f, double g[],
tnc_function *function, void *state,
- double low[], double up[], double scale[],
+ double low[], double up[], double scale[], double offset[],
int messages, int maxCGit, int maxnfeval, double eta, double stepmx,
- double accuracy, double fmin, double ftol, double rescale, int *nfeval);
+ double accuracy, double fmin, double ftol, double xtol, double pgtol,
+ double rescale, int *nfeval);
#ifdef __cplusplus
}
Modified: trunk/Lib/optimize/tnc.py
===================================================================
--- trunk/Lib/optimize/tnc.py 2007-07-23 13:01:27 UTC (rev 3184)
+++ trunk/Lib/optimize/tnc.py 2007-07-24 09:29:41 UTC (rev 3185)
@@ -1,9 +1,7 @@
-## Automatically adapted for scipy Oct 07, 2005 by convertcode.py
-
# TNC Python interface
-# @(#) $Jeannot: tnc.py,v 1.7 2004/04/02 20:40:21 js Exp $
+# @(#) $Jeannot: tnc.py,v 1.11 2005/01/28 18:27:31 js Exp $
-# Copyright (c) 2004, Jean-Sebastien Roy (js at jeannot.org)
+# Copyright (c) 2004-2005, Jean-Sebastien Roy (js at jeannot.org)
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the
@@ -33,9 +31,8 @@
value of the function, and whose second argument is the gradient of the function
(as a list of values); or None, to abort the minimization.
"""
-
from scipy.optimize import moduleTNC
-from numpy import asarray, inf
+from numpy import asarray
MSG_NONE = 0 # No messages
MSG_ITER = 1 # One line per iteration
@@ -53,21 +50,24 @@
MSG_ALL : "All messages"
}
-EINVAL = -2 # Invalid parameters (n<1)
+HUGE_VAL=1e200*1e200 # No standard representation of Infinity in Python 2.3.3
+ # FIXME: can we use inf now that we have numpy and IEEE floats?
+
INFEASIBLE = -1 # Infeasible (low > up)
LOCALMINIMUM = 0 # Local minima reach (|pg| ~= 0)
-CONVERGED = 1 # Converged (|f_n-f_(n-1)| ~= 0)
-MAXFUN = 2 # Max. number of function evaluations reach
-LSFAIL = 3 # Linear search failed
-CONSTANT = 4 # All lower bounds are equal to the upper bounds
-NOPROGRESS = 5 # Unable to progress
-USERABORT = 6 # User requested end of minimization
+FCONVERGED = 1 # Converged (|f_n-f_(n-1)| ~= 0)
+XCONVERGED = 2 # Converged (|x_n-x_(n-1)| ~= 0)
+MAXFUN = 3 # Max. number of function evaluations reach
+LSFAIL = 4 # Linear search failed
+CONSTANT = 5 # All lower bounds are equal to the upper bounds
+NOPROGRESS = 6 # Unable to progress
+USERABORT = 7 # User requested end of minimization
RCSTRINGS = {
- EINVAL : "Invalid parameters (n<1)",
INFEASIBLE : "Infeasible (low > up)",
LOCALMINIMUM : "Local minima reach (|pg| ~= 0)",
- CONVERGED : "Converged (|f_n-f_(n-1)| ~= 0)",
+ FCONVERGED : "Converged (|f_n-f_(n-1)| ~= 0)",
+ XCONVERGED : "Converged (|x_n-x_(n-1)| ~= 0)",
MAXFUN : "Max. number of function evaluations reach",
LSFAIL : "Linear search failed",
CONSTANT : "All lower bounds are equal to the upper bounds",
@@ -75,106 +75,103 @@
USERABORT : "User requested end of minimization"
}
-
# Changes to interface made by Travis Oliphant, Apr. 2004 for inclusion in
# SciPy
import optimize
approx_fprime = optimize.approx_fprime
-
def fmin_tnc(func, x0, fprime=None, args=(), approx_grad=0, bounds=None, epsilon=1e-8,
- scale=None, messages=MSG_ALL, maxCGit=-1, maxfun=None, eta=-1,
- stepmx=0, accuracy=0, fmin=0, ftol=0, rescale=-1):
+ scale=None, offset=None, messages=MSG_ALL, maxCGit=-1, maxfun=None, eta=-1,
+ stepmx=0, accuracy=0, fmin=0, ftol=-1, xtol=-1, pgtol=-1, rescale=-1):
"""Minimize a function with variables subject to bounds, using gradient
information.
- :Parameters:
+ returns (rc, nfeval, x).
- func : callable func(x, *args)
- Function to minimize.
- x0 : float
- Initial guess to minimum.
- fprime : callable fprime(x, *args)
- Gradient of func. If None, then func must return the
- function value and the gradient, e.g.
- f,g = func(x,*args).
- args : tuple
- Arguments to pass to function.
- approx_grad : bool
- If true, approximate the gradient numerically.
- bounds : list
- (min, max) pairs for each element in x, defining the
- bounds on that parameter. Use None for one of min or max
- when there is no bound in that direction
- scale : list of floats
- Scaling factors to apply to each variable. If None, the
- factors are up-low for interval bounded variables and
- 1+|x] fo the others. Defaults to None.
- messages :
- Bit mask used to select messages display during
- minimization values defined in the optimize.tnc.MSGS dict.
- defaults to optimize.tnc.MGS_ALL.
- maxCGit : int
- Maximum number of hessian*vector evaluation per main
- iteration. If maxCGit == 0, the direction chosen is
- -gradient. If maxCGit < 0, maxCGit is set to
- max(1,min(50,n/2)). Defaults to -1.
- maxfun : int
- Maximum number of function evaluation. If None, maxfun
- is set to max(1000, 100*len(x0)). Defaults to None.
- eta : float
- Severity of the line search. if < 0 or > 1, set to 0.25.
- Defaults to -1.
- stepmx : float
- Maximum step for the line search. May be increased during
- call. If too small, will be set to 10.0. Defaults to 0.
- accuracy : float
- Relative precision for finite difference calculations. If
- <= machine_precision, set to sqrt(machine_precision).
- Defaults to 0.
- fmin : float
- Minimum function value estimate. Defaults to 0.
- ftol : float
- Precision goal for the value of f in the stoping criterion
- relative to the machine precision and the value of f. If
- ftol < 0.0, ftol is set to 0.0. Defaults to 0.
- rescale : float
- Scaling factor (in log10) used to trigger rescaling. If
- 0, rescale at each iteration. If a large value, never
- rescale. If < 0, rescale is set to 1.3.
+ Inputs:
+ func : the function to minimize. Must take one argument, x and return
+ f and g, where f is the value of the function and g its
+ gradient (a list of floats).
+ if the function returns None, the minimization is aborted.
+ x0 : initial estimate (a list of floats)
+ fprime : gradient of func. If None, then func returns the function
+ value and the gradient ( f, g = func(x, *args) ).
+ Called as fprime(x, *args)
+ args : arguments to pass to function
+ approx_grad : if true, approximate the gradient numerically
+ bounds : a list of (min, max) pairs for each element in x, defining
+ the bounds on that parameter. Use None for one of min or max
+ when there is no bound in that direction
+ scale : scaling factors to apply to each variable (a list of floats)
+ if None, the factors are up-low for interval bounded variables
+ and 1+|x] fo the others.
+ defaults to None
+ offset : constant to substract to each variable
+ if None, the constant are (up+low)/2 for interval bounded
+ variables and x for the others.
+ messages : bit mask used to select messages display during minimization
+ values defined in the MSGS dict.
+ defaults to MGS_ALL
+ maxCGit : max. number of hessian*vector evaluation per main iteration
+ if maxCGit == 0, the direction chosen is -gradient
+ if maxCGit < 0, maxCGit is set to max(1,min(50,n/2))
+ defaults to -1
+ maxfun : max. number of function evaluation
+ if None, maxfun is set to max(100, 10*len(x0))
+ defaults to None
+ eta : severity of the line search. if < 0 or > 1, set to 0.25
+ defaults to -1
+ stepmx : maximum step for the line search. may be increased during call
+ if too small, will be set to 10.0
+ defaults to 0
+ accuracy : relative precision for finite difference calculations
+ if <= machine_precision, set to sqrt(machine_precision)
+ defaults to 0
+ fmin : minimum function value estimate
+ defaults to 0
+ ftol : precision goal for the value of f in the stoping criterion
+ if ftol < 0.0, ftol is set to 0.0
+ defaults to -1
+ xtol : precision goal for the value of x in the stopping criterion
+ (after applying x scaling factors)
+ if xtol < 0.0, xtol is set to sqrt(machine_precision)
+ defaults to -1
+ pgtol : precision goal for the value of the projected gradient in the
+ stopping criterion (after applying x scaling factors)
+ if pgtol < 0.0, pgtol is set to 1e-2 * sqrt(accuracy)
+ setting it to 0.0 is not recommended.
+ defaults to -1
+ rescale : f scaling factor (in log10) used to trigger f value rescaling
+ if 0, rescale at each iteration
+ if a large value, never rescale
+ if < 0, rescale is set to 1.3
- :Returns:
+ Outputs:
+ x : the solution (a list of floats)
+ nfeval : the number of function evaluations
+ rc : return code as defined in the RCSTRINGS dict
- x : list of floats
- The solution.
- nfeval : int
- The number of function evaluations.
- rc :
- Return code (corresponding message in optimize.tnc.RCSTRINGS).
+See also:
- :SeeAlso:
+ fmin, fmin_powell, fmin_cg,
+ fmin_bfgs, fmin_ncg -- multivariate local optimizers
+ leastsq -- nonlinear least squares minimizer
- - fmin, fmin_powell, fmin_cg, fmin_bfgs, fmin_ncg : multivariate
- local optimizers
+ fmin_l_bfgs_b, fmin_tnc,
+ fmin_cobyla -- constrained multivariate optimizers
- - leastsq : nonlinear least squares minimizer
+ anneal, brute -- global optimizers
- - fmin_l_bfgs_b, fmin_tnc, fmin_cobyla : constrained
- multivariate optimizers
+ fminbound, brent, golden, bracket -- local scalar minimizers
- - anneal, brute : global optimizers
+ fsolve -- n-dimenstional root-finding
- - fminbound, brent, golden, bracket : local scalar minimizers
+ brentq, brenth, ridder, bisect, newton -- one-dimensional root-finding
- - fsolve : n-dimenstional root-finding
+ fixed_point -- scalar fixed-point finder
+"""
- - brentq, brenth, ridder, bisect, newton : one-dimensional root-finding
-
- - fixed_point : scalar fixed-point finder
-
- """
-
n = len(x0)
if bounds is None:
@@ -200,31 +197,44 @@
g = fprime(x, *args)
return f, list(g)
+ """
+ low, up : the bounds (lists of floats)
+ set low[i] to -HUGE_VAL to remove the lower bound
+ set up[i] to HUGE_VAL to remove the upper bound
+ if low == None, the lower bounds are removed.
+ if up == None, the upper bounds are removed.
+ low and up defaults to None
+ """
low = [0]*n
up = [0]*n
for i in range(n):
- l,u = bounds[i]
- if l is None:
- low[i] = -inf
+ if bounds[i] is None: l, u = -HUGE_VAL, HUGE_VAL
else:
- low[i] = l
- if u is None:
- up[i] = inf
- else:
- up[i] = u
+ l,u = bounds[i]
+ if l is None:
+ low[i] = -HUGE_VAL
+ else:
+ low[i] = l
+ if u is None:
+ up[i] = HUGE_VAL
+ else:
+ up[i] = u
if scale == None:
scale = []
+ if offset == None:
+ offset = []
+
if maxfun == None:
- maxfun = max(1000, 100*len(x0))
+ maxfun = max(100, 10*len(x0))
- return moduleTNC.minimize(func_and_grad, x0, low, up, scale, messages,
- maxCGit, maxfun, eta, stepmx, accuracy,
- fmin, ftol, rescale)
+ return moduleTNC.minimize(func_and_grad, x0, low, up, scale, offset,
+ messages, maxCGit, maxfun, eta, stepmx, accuracy,
+ fmin, ftol, xtol, pgtol, rescale)
if __name__ == '__main__':
- # Examples for TNC
+ # Examples for TNC
def example():
print "Example"
@@ -239,7 +249,7 @@
return f, g
# Optimizer call
- x, nf, rc = fmin_tnc(function, [-7, 3], bounds=([-10, 10], [1, 10]))
+ rc, nf, x = fmin_tnc(function, [-7, 3], bounds=([-10, 1], [10, 10]))
print "After", nf, "function evaluations, TNC returned:", RCSTRINGS[rc]
print "x =", x
@@ -247,3 +257,100 @@
print
example()
+
+ # Tests
+ # These tests are taken from Prof. K. Schittkowski test examples for
+ # constrained nonlinear programming.
+ # http://www.uni-bayreuth.de/departments/math/~kschittkowski/home.htm
+ tests = []
+ def test1fg(x):
+ f = 100.0*pow((x[1]-pow(x[0],2)),2)+pow(1.0-x[0],2)
+ dif = [0,0]
+ dif[1] = 200.0*(x[1]-pow(x[0],2))
+ dif[0] = -2.0*(x[0]*(dif[1]-1.0)+1.0)
+ return f, dif
+ tests.append ((test1fg, [-2,1], ([-HUGE_VAL, None], [-1.5, None]), [1,1]))
+
+ def test2fg(x):
+ f = 100.0*pow((x[1]-pow(x[0],2)),2)+pow(1.0-x[0],2)
+ dif = [0,0]
+ dif[1] = 200.0*(x[1]-pow(x[0],2))
+ dif[0] = -2.0*(x[0]*(dif[1]-1.0)+1.0)
+ return f, dif
+ tests.append ((test2fg, [-2,1], ([-HUGE_VAL, None], [1.5,None]), [-1.2210262419616387,1.5]))
+
+ def test3fg(x):
+ f = x[1]+pow(x[1]-x[0],2)*1.0e-5
+ dif = [0,0]
+ dif[0] = -2.0*(x[1]-x[0])*1.0e-5
+ dif[1] = 1.0-dif[0]
+ return f, dif
+ tests.append ((test3fg, [10,1], ([-HUGE_VAL, 0.0], None), [0,0]))
+
+ def test4fg(x):
+ f = pow(x[0]+1.0,3)/3.0+x[1]
+ dif = [0,0]
+ dif[0] = pow(x[0]+1.0,2)
+ dif[1] = 1.0
+ return f, dif
+ tests.append ((test4fg, [1.125,0.125], ([1, None], [0,None]), [1,0]))
+
+ from math import *
+
+ def test5fg(x):
+ f = sin(x[0]+x[1])+pow(x[0]-x[1],2)-1.5*x[0]+2.5*x[1]+1.0
+ dif = [0,0]
+ v1 = cos(x[0]+x[1]);
+ v2 = 2.0*(x[0]-x[1]);
+
+ dif[0] = v1+v2-1.5;
+ dif[1] = v1-v2+2.5;
+ return f, dif
+ tests.append ((test5fg, [0,0], ([-1.5, 3], [-3,4]), [-0.54719755119659763, -1.5471975511965976]))
+
+ def test38fg(x):
+ f = (100.0*pow(x[1]-pow(x[0],2),2)+pow(1.0-x[0],2)+90.0*pow(x[3]-pow(x[2],2),2) \
+ +pow(1.0-x[2],2)+10.1*(pow(x[1]-1.0,2)+pow(x[3]-1.0,2)) \
+ +19.8*(x[1]-1.0)*(x[3]-1.0))*1.0e-5
+ dif = [0,0,0,0]
+ dif[0] = (-400.0*x[0]*(x[1]-pow(x[0],2))-2.0*(1.0-x[0]))*1.0e-5
+ dif[1] = (200.0*(x[1]-pow(x[0],2))+20.2*(x[1]-1.0)+19.8*(x[3]-1.0))*1.0e-5
+ dif[2] = (-360.0*x[2]*(x[3]-pow(x[2],2))-2.0*(1.0-x[2]))*1.0e-5
+ dif[3] = (180.0*(x[3]-pow(x[2],2))+20.2*(x[3]-1.0)+19.8*(x[1]-1.0))*1.0e-5
+ return f, dif
+ tests.append ((test38fg, [-3,-1,-3,-1], (([-10, 10],)*4), [1]*4))
+
+ def test45fg(x):
+ f = 2.0-x[0]*x[1]*x[2]*x[3]*x[4]/120.0
+ dif = [0]*5
+ dif[0] = -x[1]*x[2]*x[3]*x[4]/120.0
+ dif[1] = -x[0]*x[2]*x[3]*x[4]/120.0
+ dif[2] = -x[0]*x[1]*x[3]*x[4]/120.0
+ dif[3] = -x[0]*x[1]*x[2]*x[4]/120.0
+ dif[4] = -x[0]*x[1]*x[2]*x[3]/120.0
+ return f, dif
+ tests.append ((test45fg, [2]*5, ([0,1], [0,2], [0,3], [0,4], [0,5],), [1,2,3,4,5]))
+
+ def test(fg, x, bounds, xopt):
+ print "** Test", fg.__name__
+ rc, nf, x = fmin_tnc(fg, x, bounds=bounds, messages = MSG_NONE, maxfun = 200)
+ print "After", nf, "function evaluations, TNC returned:", RCSTRINGS[rc]
+ print "x =", x
+ print "exact value =", xopt
+ enorm = 0.0
+ norm = 1.0
+ for y,yo in zip(x, xopt):
+ enorm += (y-yo)*(y-yo)
+ norm += yo*yo
+ ex = pow(enorm/norm, 0.5)
+ print "X Error =", ex
+ ef = abs(fg(xopt)[0] - fg(x)[0])
+ print "F Error =", ef
+ if ef > 1e-8:
+ raise "Test "+fg.__name__+" failed"
+
+ for fg, x, bounds, xopt in tests:
+ test(fg, x, bounds, xopt)
+
+ print
+ print "** All TNC tests passed."
More information about the Scipy-svn
mailing list