Index: ThirdParty/mpqc_open/src/lib/math/optimize/qnewton.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/optimize/qnewton.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/optimize/qnewton.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,414 @@
+//
+// qnewton.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit 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 Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <math.h>
+#include <float.h>
+
+#include <util/state/stateio.h>
+#include <math/optimize/qnewton.h>
+#include <util/keyval/keyval.h>
+#include <util/misc/formio.h>
+
+using namespace std;
+using namespace sc;
+
+/////////////////////////////////////////////////////////////////////////
+// QNewtonOpt
+
+static ClassDesc QNewtonOpt_cd(
+  typeid(QNewtonOpt),"QNewtonOpt",2,"public Optimize",
+  0, create<QNewtonOpt>, create<QNewtonOpt>);
+
+QNewtonOpt::QNewtonOpt(const Ref<KeyVal>&keyval):
+  Optimize(keyval)
+{
+
+  if (function_.null()) {
+      ExEnv::err0() << "QNewtonOpt requires a function keyword" << endl;
+      abort();
+  }
+
+  init();
+
+  update_ << keyval->describedclassvalue("update");
+  if (update_.nonnull()) update_->set_inverse();
+  lineopt_ << keyval->describedclassvalue("lineopt");
+  accuracy_ = keyval->doublevalue("accuracy");
+  if (keyval->error() != KeyVal::OK) accuracy_ = 0.0001;
+  print_x_ = keyval->booleanvalue("print_x");
+  print_hessian_ = keyval->booleanvalue("print_hessian");
+  print_gradient_ = keyval->booleanvalue("print_gradient");
+  linear_ = keyval->booleanvalue("linear");
+  if (keyval->error() != KeyVal::OK) linear_ = 0;
+  restrict_ = keyval->booleanvalue("restrict");
+  if (keyval->error() != KeyVal::OK) restrict_ = 1;
+  dynamic_grad_acc_ = keyval->booleanvalue("dynamic_grad_acc");
+  if (keyval->error() != KeyVal::OK) dynamic_grad_acc_ = 1;
+  force_search_ = keyval->booleanvalue("force_search");
+  if (keyval->error() != KeyVal::OK) force_search_ = 0;
+  restart_ = keyval->booleanvalue("restart");
+  if (keyval->error() != KeyVal::OK) restart_ = 1;
+
+  RefSymmSCMatrix hessian(dimension(),matrixkit());
+  // get a guess hessian from the function
+  function()->guess_hessian(hessian);
+  
+  // see if any hessian matrix elements have been given in the input
+  if (keyval->exists("hessian")) {
+      int n = hessian.n();
+      for (int i=0; i<n; i++) {
+          if (keyval->exists("hessian",i)) {
+              for (int j=0; j<=i; j++) {
+                  double tmp = keyval->doublevalue("hessian",i,j);
+                  if (keyval->error() == KeyVal::OK) hessian(i,j) = tmp;
+                }
+            }
+        }
+    }
+  ihessian_ = function()->inverse_hessian(hessian);
+}
+
+QNewtonOpt::QNewtonOpt(StateIn&s):
+  SavableState(s),
+  Optimize(s)
+{
+  ihessian_ = matrixkit()->symmmatrix(dimension());
+  ihessian_.restore(s);
+  update_ << SavableState::restore_state(s);
+  s.get(accuracy_);
+  s.get(take_newton_step_);
+  s.get(maxabs_gradient);
+  if (s.version(::class_desc<QNewtonOpt>()) > 1) {
+    s.get(print_hessian_);
+    s.get(print_x_);
+    s.get(print_gradient_);
+  }
+  else {
+    print_hessian_ = 0;
+    print_x_ = 0;
+    print_gradient_ = 0;
+  }
+  lineopt_ << SavableState::restore_state(s);
+}
+
+QNewtonOpt::~QNewtonOpt()
+{
+}
+
+void
+QNewtonOpt::save_data_state(StateOut&s)
+{
+  Optimize::save_data_state(s);
+  ihessian_.save(s);
+  SavableState::save_state(update_.pointer(),s);
+  s.put(accuracy_);
+  s.put(take_newton_step_);
+  s.put(maxabs_gradient);
+  s.put(print_hessian_);
+  s.put(print_x_);
+  s.put(print_gradient_);
+  SavableState::save_state(lineopt_.pointer(),s);
+}
+
+void
+QNewtonOpt::init()
+{
+  Optimize::init();
+  take_newton_step_ = 1;
+  maxabs_gradient = -1.0;
+}
+
+int
+QNewtonOpt::update()
+{
+  // these are good candidates to be input options
+  const double maxabs_gradient_to_desired_accuracy = 0.05;
+  const double maxabs_gradient_to_next_desired_accuracy = 0.005;
+  const double roundoff_error_factor = 1.1;
+  
+  // the gradient convergence criterion.
+  double old_maxabs_gradient = maxabs_gradient;
+  RefSCVector xcurrent;
+  RefSCVector gcurrent;
+
+  if( dynamic_grad_acc_ ) {
+    // get the next gradient at the required level of accuracy.
+    // usually only one pass is needed, unless we happen to find
+    // that the accuracy was set too low.
+    int accurate_enough;
+    do {
+        // compute the current point
+        function()->set_desired_gradient_accuracy(accuracy_);
+        xcurrent = function()->get_x();
+        gcurrent = function()->gradient().copy();
+      
+        // compute the gradient convergence criterion now so i can see if
+        // the accuracy needs to be tighter
+        maxabs_gradient = gcurrent.maxabs();
+        // compute the required accuracy
+        accuracy_ = maxabs_gradient * maxabs_gradient_to_desired_accuracy;
+      
+        if (accuracy_ < DBL_EPSILON) accuracy_ = DBL_EPSILON;
+
+        // The roundoff_error_factor is thrown in to allow for round off making
+        // the current gcurrent.maxabs() a bit smaller than the previous,
+        // which would make the current required accuracy less than the
+        // gradient's actual accuracy and cause everything to be recomputed.
+        accurate_enough = (
+            function()->actual_gradient_accuracy()
+            <= accuracy_*roundoff_error_factor);
+
+        if (!accurate_enough) {
+          ExEnv::out0().unsetf(ios::fixed);
+          ExEnv::out0() << indent << 
+	       "NOTICE: function()->actual_gradient_accuracy() > accuracy_:\n"
+               << indent
+               << scprintf(
+                 "        function()->actual_gradient_accuracy() = %15.8e",
+                 function()->actual_gradient_accuracy()) << endl << indent
+               << scprintf(
+                 "                                     accuracy_ = %15.8e",
+                 accuracy_) << endl;
+        }
+     } while(!accurate_enough);
+    // increase accuracy, since the next gradient will be smaller
+    accuracy_ = maxabs_gradient * maxabs_gradient_to_next_desired_accuracy;
+  }
+  else {
+    xcurrent = function()->get_x();
+    gcurrent = function()->gradient().copy();
+  }
+
+  if (old_maxabs_gradient >= 0.0 && old_maxabs_gradient < maxabs_gradient) {
+    ExEnv::out0() << indent
+	 << scprintf("NOTICE: maxabs_gradient increased from %8.4e to %8.4e",
+                     old_maxabs_gradient, maxabs_gradient) << endl;
+  }
+
+  // update the hessian
+  if(update_.nonnull())
+    update_->update(ihessian_,function(),xcurrent,gcurrent);
+
+  conv_->reset();
+  conv_->get_grad(function());
+  conv_->get_x(function());
+
+  // compute the quadratic step
+  RefSCVector xdisp = -1.0*(ihessian_ * gcurrent);
+  RefSCVector xnext = xcurrent + xdisp;
+
+  // either do a lineopt or check stepsize
+  double tot;
+  if(lineopt_.nonnull()) {
+    if (dynamic_cast<Backtrack*>(lineopt_.pointer()) != 0) {
+      // The Backtrack line search is a special case.
+
+      // perform a search
+      double factor;
+      if( n_iterations_ == 0 && force_search_ ) 
+        factor = lineopt_->set_decrease_factor(1.0);
+      lineopt_->init(xdisp,function());
+      // reset value acc here so line search "precomputes" are 
+      // accurate enough for subsequent gradient evals
+      function()->set_desired_value_accuracy(accuracy_/100);
+      int acceptable = lineopt_->update();
+      if( n_iterations_ == 0 && force_search_ )
+        lineopt_->set_decrease_factor( factor );
+
+      if( !acceptable ) {
+        if( force_search_ ) factor = lineopt_->set_decrease_factor(1.0);
+
+        // try a new guess hessian
+        if( restart_ ) {
+          ExEnv::out0() << endl << indent << 
+            "Restarting Hessian approximation" << endl;
+          RefSymmSCMatrix hessian(dimension(),matrixkit());
+          function()->guess_hessian(hessian);
+          ihessian_ = function()->inverse_hessian(hessian);
+          xdisp = -1.0 * (ihessian_ * gcurrent);
+          lineopt_->init(xdisp,function());
+          acceptable = lineopt_->update();
+        }
+      
+        // try steepest descent direction
+        if( !acceptable ) {
+          ExEnv::out0() << endl << indent << 
+            "Trying steepest descent direction." << endl;
+          xdisp = -1.0 * gcurrent;
+          lineopt_->init(xdisp,function());
+          acceptable = lineopt_->update();
+        }
+
+        // give up and use steepest descent step
+        if( !acceptable ) {
+          ExEnv::out0() << endl << indent << 
+            "Resorting to unscaled steepest descent step." << endl;
+          function()->set_x(xcurrent + xdisp);
+          Ref<NonlinearTransform> t = function()->change_coordinates();
+          apply_transform(t);
+        }
+
+        if( force_search_ ) lineopt_->set_decrease_factor( factor );
+      }
+    }
+    else {
+      // All line searches other than Backtrack use this
+      ExEnv::out0() << indent
+                    << "......................................."
+                    << endl
+                    << indent
+                    << "Starting line optimization."
+                    << endl;
+      lineopt_->init(xdisp,function());
+      int nlineopt = 0;
+      int maxlineopt = 3;
+      for (int ilineopt=0; ilineopt<maxlineopt; ilineopt++) {
+        double maxabs_gradient = function()->gradient()->maxabs();
+
+        int converged = lineopt_->update();
+
+        ExEnv::out0() << indent
+                      << "Completed line optimization step " << ilineopt+1
+                      << (converged?" (converged)":" (not converged)")
+                      << endl
+                      << indent
+                      << "......................................."
+                      << endl;
+
+        if (converged) break;
+
+        // Improve accuracy, since we might be able to reuse the next
+        // gradient for the next quasi-Newton step.
+        if (dynamic_grad_acc_)  {
+          accuracy_ = maxabs_gradient*maxabs_gradient_to_next_desired_accuracy;
+          function()->set_desired_gradient_accuracy(accuracy_);
+        }
+      }
+    }
+    xnext = function()->get_x();
+    xdisp = xnext - xcurrent;
+    tot = sqrt(xdisp.scalar_product(xdisp));
+  }
+  else {
+
+    tot = sqrt(xdisp.scalar_product(xdisp));
+
+    if ( tot > max_stepsize_ ) {
+      if( restrict_ ) {
+	double scal = max_stepsize_/tot;
+	ExEnv::out0() << endl << indent <<
+	  scprintf("stepsize of %f is too big, scaling by %f",tot,scal)
+	  << endl;
+	xdisp.scale(scal);
+	tot *= scal;
+      }
+      else {
+	ExEnv::out0() << endl << indent <<
+	  scprintf("stepsize of %f is too big, but scaling is disabled",tot)
+	  << endl;
+      }
+    }
+    xnext = xcurrent + xdisp;
+  }
+
+  if (print_hessian_) {
+    RefSymmSCMatrix hessian = ihessian_.gi();
+    ExEnv::out0() << indent << "hessian = [" << endl;
+    ExEnv::out0() << incindent;
+    int n = hessian.n();
+    for (int i=0; i<n; i++) {
+      ExEnv::out0() << indent << "[";
+      for (int j=0; j<=i; j++) {
+        ExEnv::out0() << scprintf(" % 10.6f",double(hessian(i,j)));
+      }
+      ExEnv::out0() << " ]" << endl;
+    }
+    ExEnv::out0() << decindent;
+    ExEnv::out0() << indent << "]" << endl;
+  }
+  if (print_x_) {
+    int n = xcurrent.n();
+    ExEnv::out0() << indent << "x = [";
+    for (int i=0; i<n; i++) {
+      ExEnv::out0() << scprintf(" % 16.12f",double(xcurrent(i)));
+    }
+    ExEnv::out0() << " ]" << endl;
+  }
+  if (print_gradient_) {
+    int n = gcurrent.n();
+    ExEnv::out0() << indent << "gradient = [";
+    for (int i=0; i<n; i++) {
+      ExEnv::out0() << scprintf(" % 16.12f",double(gcurrent(i)));
+    }
+    ExEnv::out0() << " ]" << endl;
+  }
+
+  // check for convergence
+  conv_->set_nextx(xnext);
+  int converged = conv_->converged();
+  if (converged) return converged;
+
+  ExEnv::out0() << indent
+		<< scprintf("taking step of size %f", tot) << endl;
+  ExEnv::out0() << indent << "Optimization iteration " 
+		<< n_iterations_ + 1 << " complete" << endl;
+  ExEnv::out0() << indent
+    << "//////////////////////////////////////////////////////////////////////"
+    << endl;
+
+  if( lineopt_.null() ) {
+    function()->set_x(xnext);
+    Ref<NonlinearTransform> t = function()->change_coordinates();
+    apply_transform(t);
+  }
+
+  if( dynamic_grad_acc_ ) 
+    function()->set_desired_gradient_accuracy(accuracy_);
+
+  return converged;
+}
+
+void
+QNewtonOpt::apply_transform(const Ref<NonlinearTransform> &t)
+{
+  if (t.null()) return;
+  Optimize::apply_transform(t);
+  if (lineopt_.nonnull()) lineopt_->apply_transform(t);
+  if (ihessian_.nonnull()) t->transform_ihessian(ihessian_);
+  if (update_.nonnull()) update_->apply_transform(t);
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "ETS"
+// End:
