[Scipy-svn] r3185 - in trunk/Lib/optimize: . tnc

scipy-svn@scip... scipy-svn@scip...
Tue Jul 24 04:30:21 CDT 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@jeannot.org)
+Copyright (c) 2002-2005, Jean-Sebastien Roy (js@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@jeannot.org), 2002-2004
+# Version 1.3
+# Copyright J.S. Roy (js@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@jeannot.org)
- * 
+ * Copyright (c) 2004-2005, Jean-Sebastien Roy (js@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@jeannot.org)
- * 
+ * Copyright (c) 2002-2005, Jean-Sebastien Roy (js@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@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, &gtest1, &gtest2, &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, &gtest1, &gtest2, &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@jeannot.org)
- * 
+ * Copyright (c) 2002-2005, Jean-Sebastien Roy (js@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@jeannot.org)
+# Copyright (c) 2004-2005, Jean-Sebastien Roy (js@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