ExampleSOCP.cpp

/* This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/


/** @example ExampleSOCP.cpp
    @relates SOCP
    @ingroup g_examples
    @brief   Example for Second Order Cone Programming.
             Similar to Matlab example provided by M. S. Lobo et al.
             from http://www.stanford.edu/~boyd/old_software/socp/

             Minimize x + y subject to (x+1)^2 + (y-1)^2 <= 1.

    @author  esquivel 05/2012
*/

#include <bias_config.h>

#include <MathAlgo/SOCP.hh>
#include <Base/Math/Random.hh>

using namespace BIAS;
using namespace std;

int main(int args, char **arg)
{
  // Specify SOCP instance with single constraint
  const int m = 2; // number of variables
  const int n = 2; // dimension of constraint
  Vector<double> f(m);
  Matrix<double> A(n, m, MatrixZero);
  Vector<double> b(n), c(m);
  double d;
  f[0] = f[1] = 1;
  A[0][0] = A[1][1] = 1;
  b[0] = 1;
  b[1] = -1;
  c[0] = c[1] = 0;
  d = 1;

  // Create initial primal and dual solution
  Vector<double> x(m), z(n);
  double w;
  x[0] = -0.5;
  x[1] = 1;
  z[0] = 1;
  z[1] = 1;
  w = 4;

  // Solve SOCP instance
  SOCP socp;
  socp.SetNu(5.0);
  socp.SetMaxIterations(50);
  socp.SetRelativeTolerance(0.0);
  socp.SetAbsoluteTolerance(1e-4);
  socp.SetTargetValue(0.0);
  int res = socp.Compute(m, n, f, A, b, c, d, x, z, w);

  // Show results
  cout << "Return code of SOCP solver = " << res << endl;
  cout << "Primal problem : " << endl;
  cout << "  Final primal solution x = " << x << endl;
  cout << "  Target value of f'x = " << f.ScalarProduct(x) << endl;
  cout << "  Constraint evaluation for solution x : " << endl
       << "    |Ax + b| = " << BIAS::Vector<double>(A*x+b).NormL2()
       << " <= c'x + d = " << c.ScalarProduct(x) + d << endl;
  cout << "Dual problem : " << endl;
  cout << "  Final dual solution z = " << z << ", w = " << w << endl;
  cout << "  Target value of -(b'z + dw) = " << -(b.ScalarProduct(z) + d*w) << endl;
  cout << "  Constraint evaluation for solution z, w : " << endl
       << "    A'z + wc = " << BIAS::Vector<double>(A.Transpose()*z + w*c)
       << " = f = " << f << endl;
  cout << "    |z| = " << z.NormL2() << " <= w = " << w << endl;
  return 0;
}