# HG changeset patch
# User jwe
# Date 1114094811 0
# Node ID 41273fff034ded172de6f338ecce565cd265244a
# Parent  b98bf1d70a0ae821f87fd52cfc11b5638468b751
[project @ 2005-04-21 14:46:50 by jwe]

diff -r b98bf1d70a0a -r 41273fff034d src/ChangeLog
--- a/src/ChangeLog	Thu Apr 21 14:43:20 2005 +0000
+++ b/src/ChangeLog	Thu Apr 21 14:46:51 2005 +0000
@@ -1,3 +1,13 @@
+2005-04-21  John W. Eaton  <jwe@octave.org>
+
+	* DLD-FUNCTIONS/__qp__.cc: New file.
+	* Makefile.in (DLD_XSRC): Add it to the list.
+
+2005-04-20  John W. Eaton  <jwe@octave.org>
+
+	* lex.l (IDENT): Allow $ in identifiers.
+	* utils.cc (valid_identifier): Likewise.
+
 2005-04-19  John W. Eaton  <jwe@octave.org>
 
 	* toplev.cc (Fsystem): Move enum exec_type declaration to file
diff -r b98bf1d70a0a -r 41273fff034d src/DLD-FUNCTIONS/__qp__.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/__qp__.cc	Thu Apr 21 14:46:51 2005 +0000
@@ -0,0 +1,531 @@
+/*
+
+Copyright (C) 2000, 2001, 2004, 2005 Gabriele Pannocchia
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 2, or (at your option) any
+later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+#include <cfloat>
+
+#include <iostream>
+
+#include <octave/oct.h>
+#include <octave/parse.h>
+#include <octave/EIG.h>
+#include <octave/dbleCHOL.h>
+#include <octave/dbleSVD.h>
+
+static inline double
+ABS (double x)
+{
+  return x < 0 ? -x : x;
+}
+
+static Matrix
+null (const Matrix& A, int& rank)
+{
+  Matrix retval;
+
+  rank = 0;
+
+  if (! A.is_empty ())
+    {
+      SVD A_svd (A);
+
+      DiagMatrix S = A_svd.singular_values ();
+
+      ColumnVector s = S.diag ();
+
+      Matrix V = A_svd.right_singular_matrix ();
+
+      int A_nr = A.rows ();
+      int A_nc = A.cols ();
+
+      int tmp = A_nr > A_nc ? A_nr : A_nc;
+
+      double tol = tmp * s(0) * DBL_EPSILON;
+
+      int n = s.length ();
+
+      for (int i = 0; i < n; i++)
+	{
+	  if (s(i) > tol)
+	    rank++;
+	}
+
+      if (rank < A_nc)
+	retval = V.extract (0, rank, A_nc-1, A_nc-1);
+      else
+	retval.resize (A_nc, 0);
+    }
+
+  return retval;
+}
+
+static Matrix
+cholinv (const Matrix& R)
+{
+  // R should be the result of calling chol on a symmetric positive
+  // definite matrix.
+  int n = R.rows ();
+  Matrix L = R.transpose ();
+  ColumnVector d = L.diag ();
+  ColumnVector tmp (n);
+  for (int k = 0; k < n; k++)
+    {
+      for (int j = 0; j < n; j++)
+	L(j,k) = L(j,k) / d(k);
+      tmp(k) = 1.0/(d(k)*d(k));
+    }
+  DiagMatrix Dinv (tmp);
+  Matrix invL = L.inverse ();
+  return invL.transpose () * Dinv * invL;
+}
+
+static int
+qp (const Matrix& H, const ColumnVector& q,
+    const Matrix& Aeq, const ColumnVector& beq,
+    const Matrix& Ain, const ColumnVector& bin,
+    int maxit,
+    ColumnVector& x, ColumnVector& lambda, int& iter)
+{
+  int info = 0;
+
+  iter = 0;
+
+  double rtol = sqrt (DBL_EPSILON);
+
+  // Problem dimension.
+  int n = x.length ();
+
+  // Dimension of constraints.
+  int n_eq = beq.length ();
+  int n_in = bin.length ();
+
+  // Filling the current active set.
+
+  int n_act = n_eq;
+
+  int n_tot = n_eq + n_in;
+
+  // Equality constraints come first.  We won't check the sign of the
+  // Lagrange multiplier for those.
+
+  Matrix Aact = Aeq;
+  ColumnVector bact = beq;
+  ColumnVector Wact;
+
+  if (n_in > 0)
+    {
+      ColumnVector res = Ain*x - bin;
+
+      for (int i = 0; i < n_in; i++)
+	{
+	  res(i) /= (1.0 + ABS (bin(i)));
+
+	  if (res(i) < rtol)
+	    {
+	      n_act++;
+	      Aact = Aact.stack (Ain.row (i));
+	      bact.resize (n_act, bin(i));
+	      Wact.resize (n_act, i);
+	    }
+	}
+    }
+
+  // Computing the ???
+
+  EIG eigH (H);
+  ColumnVector eigenvalH = real (eigH.eigenvalues ());
+  Matrix eigenvecH = real (eigH.eigenvectors ());
+  double minReal = eigenvalH.min ();
+  int indminR = 0;
+  for (int i = 0; i < n; i++)
+    {
+      if (minReal == eigenvalH(i))
+	{
+	  indminR = i;
+	  break;
+	}
+    }
+
+  bool done = false;
+
+  double alpha = 0.0;
+
+  Matrix R;
+  Matrix Y (n, 0, 0.0);
+
+  ColumnVector g (n, 0.0);
+  ColumnVector p (n, 0.0);
+
+  ColumnVector lambda_tmp (n_in, 0.0);
+
+  while (! done)
+    {
+      iter++;
+
+      // Current Gradient
+      // g = q + H * x;
+
+      g = q + H * x;
+
+      if (n_act == 0)
+	{
+	  // There are no active constraints.
+
+	  if (minReal > 0.0)
+	    {
+	      // Inverting the Hessian.  Using the Cholesky
+	      // factorization since the Hessian is positive
+	      // definite.
+
+	      CHOL cholH (H);
+
+	      R = cholH.chol_matrix ();
+
+	      Matrix Hinv = cholinv (R);
+
+	      // Computing the unconstrained step.
+	      // p = -Hinv * g;
+
+	      p = -Hinv * g;
+
+	      info = 0;
+	    }
+	  else
+	    {
+	      // Finding the negative curvature of H.
+
+	      p = eigenvecH.column (indminR);
+
+	      // Following the negative curvature of H.
+
+	      if (p.transpose () * g > DBL_EPSILON)
+	        p = -p;
+
+	      info = 1;
+	    }
+
+	  // Multipliers are zero.
+          lambda_tmp.fill (0.0);
+	}
+      else
+        {
+	  // There are active constraints.
+
+	  // Computing the null space.
+
+	  int rank;
+
+	  Matrix Z = null (Aact, rank);
+
+	  int dimZ = n - rank;
+
+	  // XXX FIXME XXX -- still remain to handle the case of
+	  // non-full rank active set matrix.
+
+	  // Computing the Y matrix (orthogonal to Z)
+	  Y = Aact.pseudo_inverse ();
+
+	  // Reduced Hessian
+	  Matrix Zt = Z.transpose ();
+	  Matrix rH = Zt * H * Z;
+
+	  int pR = 0;
+
+	  Matrix tR;
+
+	  if (dimZ > 0)
+	    {
+	      // Computing the Cholesky factorization (pR = 0 means
+	      // that the reduced Hessian was positive definite).
+
+	      CHOL cholrH (rH, pR);
+	      Matrix tR = cholrH.chol_matrix ();
+	      if (pR == 0)
+		R = tR;
+	    }
+
+	  if (pR == 0)
+	    {
+	      info = 0;
+
+	      // Computing the step pz. 
+	      if (dimZ > 0)
+		{
+		  // Using the Cholesky factorization to invert rH
+
+		  Matrix rHinv = cholinv (R);
+
+		  ColumnVector pz = -rHinv * Zt * g;
+
+		  // Global step.
+		  p = Z * pz;
+		}
+	      else
+		{
+		  // Global step.
+		  p.fill (0.0);
+		}
+	    }
+	  else
+	    {
+	      info = 1;
+
+	      // Searching for the most negative curvature.
+
+	      EIG eigrH (rH);
+	      ColumnVector eigenvalrH = real (eigrH.eigenvalues ());
+	      Matrix eigenvecrH = real (eigrH.eigenvectors ());
+	      double mRrH = eigenvalrH.min ();
+	      indminR = 0;
+	      for (int i = 0; i < n; i++)
+		{
+		  if (mRrH == eigenvalH(i))
+		    {
+		      indminR = i;
+		      break;
+		    }
+		}
+
+	      ColumnVector eVrH = eigenvecrH.column (indminR);
+
+	      // Computing the step pz.
+	      p = Z * eVrH;
+
+	      if (p.transpose () * g > DBL_EPSILON)
+		p = -p;
+	    }
+	}
+
+      // Checking the step-size.
+      ColumnVector abs_p (n);
+      for (int i = 0; i < n; i++)
+	abs_p(i) = ABS (p(i));
+      double max_p = abs_p.max ();
+
+      if (max_p < rtol)
+	{
+	  // The step is null.  Checking constraints.
+	  if (n_act - n_eq == 0)
+	    // Solution is found because no inequality
+	    // constraints are active.
+	    done = true;
+	  else
+	    {
+	      // Computing the multipliers only for the inequality
+	      // constraints that are active.  We do NOT compute
+	      // multipliers for the equality constraints.
+	      Matrix Yt = Y.transpose ();
+	      Yt = Yt.extract_n (n_eq, 0, n_act-n_eq, n);
+	      lambda_tmp = Yt * (g + H * p);
+	      if (n_act - n_eq < n_in)
+		{
+		  lambda_tmp.resize (n_in, 0.0);
+
+		  for (int i = n_act-n_eq; i < n_in; i++)
+		    lambda_tmp(i) = 0;
+		}
+
+	      // Checking the multipliers.  We remove the most
+	      // negative from the set (if any).
+	      double min_lambda = lambda_tmp.min ();
+	      if (min_lambda >= 0)
+		{
+		  // Solution is found.
+		  done = true;
+		}
+	      else
+		{
+		  int which_eig = 0;
+		  for (int i = 0; i < n_act; i++)
+		    {
+		      if (lambda_tmp(i) == min_lambda)
+			{
+			  which_eig = i;
+			  break;
+			}
+		    }
+
+		  // At least one multiplier is negative, we
+		  // remove it from the set.
+
+		  for (int i = which_eig; i < n_act - n_eq; i++)
+		    {
+		      Wact(i) = Wact(i+1);
+		      for (int j = 0; j < n; j++)
+			Aact(n_eq+i,j) = Aact(n_eq+i+1,j);
+		      bact(n_eq+i) = bact(n_eq+i+1);
+		    }
+		  n_act--;
+
+		  // Resizing the active set.
+		  Wact.resize (n_act-n_eq);
+		  bact.resize (n_act);
+		  Aact.resize (n_act, n);
+		}
+	    }
+	}
+      else
+	{
+	  // The step is not null.
+	  if (n_act - n_eq == n_in)
+	    {
+	      // All inequality constraints were active.  We can
+	      // add the whole step.
+	      x += p;
+	    }
+	  else
+	    {
+	      // Some constraints were not active.  Checking if
+	      // there is a blocking constraint.
+	      alpha = 1.0;
+	      int is_block = -1;
+
+	      for (int i = 0; i < n_in; i++)
+		{
+		  bool found = false;
+
+		  for (int j = 0; j < n_act-n_eq; j++)
+		    {
+		      if (Wact(j) == i)
+			{
+			  found = true;
+			  break;
+			}
+		    }
+
+		  if (! found)
+		    {
+		      // The i-th constraint was not in the set.  Is it a
+		      // blocking constraint?
+
+		      RowVector tmp_row = Ain.row (i);
+		      double tmp = tmp_row * p;
+		      double res = tmp_row * x;
+
+		      if (tmp < 0.0)
+		        {
+			  double alpha_tmp = (bin(i) - res) / tmp;
+
+			  if (alpha_tmp < alpha)
+			    {
+			      alpha = alpha_tmp;
+			      is_block = i;
+			    }
+			}
+		    }
+		}
+
+	      // In is_block there is the index of the blocking
+	      // constraint (if any).
+	      if (is_block >= 0)
+		{
+		  // There is a blocking constraint (index in
+		  // is_block) which is added to the active set.
+		  n_act++;
+		  Aact = Aact.stack (Ain.row (is_block));
+		  bact.resize (n_act, bin(is_block));
+		  Wact.resize (n_act, is_block);
+
+		  // Adding the reduced step
+		  x += alpha * p;
+		}
+	      else
+		{
+		  // There are no blocking constraints.  Adding the
+		  // whole step.
+		  x += alpha * p;
+		}
+	    }
+	}
+
+      if (iter == maxit)
+	{
+	  done = true;
+	  // warning ("qp_main: maximum number of iteration reached");
+	  info = 3;
+	}
+    }
+
+  lambda_tmp = Y.transpose () * (g + H * p);
+
+  // Reordering the Lagrange multipliers.
+
+  lambda.resize (n_tot);
+  lambda.fill (0.0);
+  for (int i = 0; i < n_eq; i++)
+    lambda(i) = lambda_tmp(i);
+
+  for (int i = n_eq; i < n_tot; i++)
+    {
+      for (int j = 0; j < n_act-n_eq; j++)
+	{
+	  if (Wact(j) == i)
+	    {
+	      lambda(i) = lambda_tmp(n_eq+j);
+	      break;
+	    }
+	}
+    }
+
+  return info;
+}
+
+DEFUN_DLD (__qp__, args, ,
+  "[x, lambda, info, iter] = __qp__ (x0, H, q, Aeq, beq, Ain, bin, maxit)")
+{
+  octave_value_list retval;
+
+  if (args.length () == 8)
+    {
+      const ColumnVector x0  (args(0) . vector_value ());
+      const Matrix H         (args(1) . matrix_value ());
+      const ColumnVector q   (args(2) . vector_value ());
+      const Matrix Aeq       (args(3) . matrix_value ());
+      const ColumnVector beq (args(4) . vector_value ());
+      const Matrix Ain       (args(5) . matrix_value ());
+      const ColumnVector bin (args(6) . vector_value ());
+      const int maxit        (args(7) . int_value ());
+
+      if (! error_state)
+	{
+	  int iter = 0;
+
+	  // Copying the initial guess in the working variable
+	  ColumnVector x = x0;
+
+	  // Reordering the Lagrange multipliers
+	  ColumnVector lambda;
+
+	  int info = qp (H, q, Aeq, beq, Ain, bin, maxit, x, lambda, iter);
+
+	  retval(3) = iter;
+	  retval(2) = info;
+	  retval(1) = lambda;
+	  retval(0) = x;
+	}
+      else
+	error ("__qp__: invalid arguments");
+    }
+  else
+    print_usage ("__qp__");
+
+  return retval;
+}
diff -r b98bf1d70a0a -r 41273fff034d src/Makefile.in
--- a/src/Makefile.in	Thu Apr 21 14:43:20 2005 +0000
+++ b/src/Makefile.in	Thu Apr 21 14:46:51 2005 +0000
@@ -47,10 +47,10 @@
 	eig.cc expm.cc fft.cc fft2.cc fftn.cc fftw_wisdom.cc \
 	filter.cc find.cc fsolve.cc gammainc.cc gcd.cc getgrent.cc \
 	getpwent.cc getrusage.cc givens.cc hess.cc inv.cc kron.cc \
-	lpsolve.cc lsode.cc lu.cc minmax.cc pinv.cc qr.cc \
+	lpsolve.cc lsode.cc lu.cc luinc.cc minmax.cc pinv.cc qr.cc \
 	quad.cc qz.cc rand.cc schur.cc sort.cc sparse.cc spdet.cc \
 	splu.cc spparms.cc sqrtm.cc svd.cc syl.cc time.cc gplot.l \
-	__glpk__.cc luinc.cc
+	__glpk__.cc __qp__.cc
 
 DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))
 
diff -r b98bf1d70a0a -r 41273fff034d src/lex.l
--- a/src/lex.l	Thu Apr 21 14:43:20 2005 +0000
+++ b/src/lex.l	Thu Apr 21 14:46:51 2005 +0000
@@ -286,7 +286,7 @@
 NOT	((\~)|(\!))
 POW     ((\*\*)|(\^))
 EPOW    (\.{POW})
-IDENT	([_a-zA-Z][_a-zA-Z0-9]*)
+IDENT	([_$a-zA-Z][_$a-zA-Z0-9]*)
 EXPON	([DdEe][+-]?{D}+)
 NUMBER	(({D}+\.?{D}*{EXPON}?)|(\.{D}+{EXPON}?)|(0[xX][0-9a-fA-F]+))
 %%
diff -r b98bf1d70a0a -r 41273fff034d src/utils.cc
--- a/src/utils.cc	Thu Apr 21 14:43:20 2005 +0000
+++ b/src/utils.cc	Thu Apr 21 14:46:51 2005 +0000
@@ -83,11 +83,11 @@
 bool
 valid_identifier (const char *s)
 {
-  if (! s || ! (isalpha (*s) || *s == '_'))
+  if (! s || ! (isalpha (*s) || *s == '_' || *s == '$'))
      return false;
 
   while (*++s != '\0')
-    if (! (isalnum (*s) || *s == '_'))
+    if (! (isalnum (*s) || *s == '_' || *s == '$'))
       return false;
 
   return true;