/home/coin/SVN-release/OS-2.4.0/Couenne/src/standardize/analyzeSparsity.cpp

Go to the documentation of this file.
00001 /* $Id: analyzeSparsity.cpp 732 2011-07-03 20:06:50Z pbelotti $
00002  *
00003  * Name:    analyzeSparsity.cpp
00004  * Author:  Pietro Belotti
00005  * Purpose: return one or more exprGroup/exprQuad based on sparsity of
00006  *          original one
00007  *
00008  * (C) Carnegie-Mellon University, 2007-10.
00009  * This file is licensed under the Eclipse Public License (EPL)
00010  */
00011 
00012 #include <stdio.h>
00013 #include <map>
00014 #include <set>
00015 
00016 #include "CouenneExprQuad.hpp"
00017 #include "CouenneTypes.hpp"
00018 #include "CouenneProblem.hpp"
00019 #include "CouenneExprAux.hpp"
00020 #include "CouenneExprMul.hpp"
00021 #include "CouenneExprPow.hpp"
00022 #include "CouenneLQelems.hpp"
00023 
00024 using namespace Couenne;
00025 
00026 #define MIN_DENSITY 0.5
00027 
00031 void CouenneProblem::analyzeSparsity (CouNumber c0, 
00032                                       LinMap &lmap,
00033                                       QuadMap &qmap) {
00034 
00035   if (qmap.Map().size () == 0) return;
00036 
00037   // simple technique: if number of elements in quadratic map is more
00038   // than a given fraction of n^2, then turn it into an exprQuad,
00039   // otherwise break it down.
00040 
00041   std::set <int> occur;
00042   unsigned int nsquares = 0;
00043 
00044   for (std::map <std::pair <int,int>, CouNumber>::iterator i = qmap.Map().begin ();
00045        i != qmap.Map().end (); ++i) {
00046 
00047     int
00048       first  = i -> first.first,
00049       second = i -> first.second;
00050 
00051     if (occur.find (first) == occur.end ()) 
00052       occur.insert (first);
00053 
00054     if (first != second) {
00055       if (occur.find (second) == occur.end ())
00056         occur.insert (second);
00057     } else nsquares++;
00058   }
00059 
00060   if (jnlst_ -> ProduceOutput (Ipopt::J_ALL, J_REFORMULATE)) {
00061     printf ("qmap has %d element, occur has %d, md*s*(s+1)/2 = %g\n", 
00062             qmap.Map().size (), 
00063             occur.size (),
00064             MIN_DENSITY * (double) (occur.size ()) * ((double) (occur.size ()) + 1.) / 2);
00065   }
00066 
00067   int nterms = occur.size ();
00068   
00069   if (useQuadratic_ &&
00070       (((qmap.Map().size () >= MIN_DENSITY * nterms * (nterms+1) / 2) && 
00071         (nterms >= 2))
00072       //|| (nsquares > nterms/2)
00073        || (nsquares >= occur.size ()))
00074       )
00075     return; // keep current exprQuad structure
00076 
00077   // flatten exprQuad to a sum of terms (disaggregate). This is while
00078   // we are testing exprQuad's
00079 
00080   for (std::map <std::pair <int,int>, CouNumber>::iterator i = qmap.Map().begin ();
00081        i != qmap.Map().end (); ++i) {
00082 
00083     int indI = i -> first.first,
00084         indJ = i -> first.second;
00085 
00086     exprAux *aux = (indI != indJ) ? 
00087       addAuxiliary 
00088       (new exprMul (new exprClone (Var (indI)),
00089                     new exprClone (Var (indJ)))) : 
00090       addAuxiliary 
00091       (new exprPow (new exprClone (Var (indI)),
00092                     new exprConst (2.)));
00093 
00094     //    aux -> print (); printf (" := "); aux -> Image () -> print (); printf ("\n");
00095 
00096     lmap.insert (aux -> Index (), i -> second);
00097   }
00098 
00099   if (qmap.Map().size () == 1) {
00100 
00101     // very simple case: we have a linear term plus a single bilinear
00102     // x*y (or square x^2) term. 
00103   }
00104 
00105   qmap.Map().erase (qmap.Map().begin (), qmap.Map().end ());
00106 
00107   // in general, decompose qmap+lmap into several (qmap+lmap)s so that
00108   // each corresponds to an exprQuad to be transformed into a single
00109   // auxiliary variable
00110 
00111   // build graph and look for components -- TODO
00112 }

Generated on Thu Sep 22 03:05:59 2011 by  doxygen 1.4.7