source: trunk/Couenne/src/standardize/linStandardize.cpp @ 91

Last change on this file since 91 was 91, checked in by pbelotti, 11 years ago

fixed serious bug. Failed to find optimum due to duplicate variable...

File size: 3.4 KB
Line 
1/*
2 * Name:    linStandardize.cpp
3 * Author:  Pietro Belotti
4 * Purpose: standardize sum expressions (expr{Sum,Sub,Quad,Group,Opp})
5 *
6 * (C) Carnegie-Mellon University, 2007.
7 * This file is licensed under the Common Public License (CPL)
8 */
9
10#include "CouenneProblem.hpp"
11#include "exprSum.hpp"
12#include "exprSub.hpp"
13#include "exprOpp.hpp"
14#include "exprMul.hpp"
15#include "exprPow.hpp"
16#include "exprGroup.hpp"
17#include "exprQuad.hpp"
18#include "lqelems.hpp"
19
20//#define DEBUG
21
22/// standardization of linear exprOp's
23exprAux *CouenneProblem::linStandardize (bool addAux, 
24                                         CouNumber c0, 
25                                         LinMap  &lmap,
26                                         QuadMap &qmap) {
27
28  // check if quadratic forms are dense enough ///////////////////////////////////////////
29
30  analyzeSparsity (c0, lmap, qmap);
31
32  ////////////////////////////////////////////////////////////////////////////////////////
33
34  int  nq = qmap.Map().size (),     /// data for exprQuad
35      *qi = new int [nq+1], 
36      *qj = new int [nq+1];
37
38  CouNumber *qc = new CouNumber [nq];
39
40  int 
41     nl = lmap.Map().size(),      /// data for exprGroup
42    *li = new int [nl+1];
43
44  CouNumber *lc = new CouNumber [nl];
45
46  // terminate arrays with negative index
47  qi [nq] = li [nl] = -1; 
48
49  std::map <int, CouNumber>::iterator lit = lmap.Map().begin (); 
50
51  // fill in arrays for linear part
52  for (int i=0; i<nl; i++, lit++) {
53
54    li [i] = lit -> first;
55    lc [i] = lit -> second;
56  }
57
58  std::map <std::pair <int, int>, CouNumber>::iterator qit = qmap.Map().begin (); 
59
60  // fill in arrays for quadratic part
61  for (int i=0; i < nq; i++, qit++) {
62    qi [i] = qit -> first. first;
63    qj [i] = qit -> first. second;
64    qc [i] = qit -> second;
65  }
66
67  nl = lmap.Map().size ();
68  nq = qmap.Map().size ();
69
70  // particular cases ///////////////////////////////////////////////////////////
71
72  expression *ret;
73
74  // a constant
75  if ((nq==0) && (nl==0)) 
76
77    ret = new exprClone (addAuxiliary (new exprConst (c0))); // a constant auxiliary? FIX!
78
79  else if ((nq==0) && (fabs (c0) < COUENNE_EPS) && (nl==1)) { // a linear monomial, cx
80
81    if (fabs (*lc - 1) < COUENNE_EPS) 
82      ret    = new exprClone (Var (*li));
83    else ret = new exprMul (new exprConst (*lc), new exprClone (Var (*li)));
84
85  } else if ((nl==0) && (fabs (c0) < COUENNE_EPS) && (nq==1)) { 
86
87    // a bilinear/quadratic monomial, cx^2 or cxy
88
89    expression *quad;
90
91    if (*qi == *qj) quad = new exprPow (new exprClone (Var (*qi)), new exprConst (2.));
92    else            quad = new exprMul (new exprClone (Var (*qi)), 
93                                        new exprClone (Var (*qj)));
94
95    if (fabs (*qc - 1) < COUENNE_EPS) 
96      ret    = quad;
97    else {
98      quad = addAuxiliary (quad);
99      ret  = new exprMul (new exprConst (*qc), new exprClone (quad));
100    }
101
102  } else {
103
104    // general case ///////////////////////////////////////////////////////////////
105
106    std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
107    indcoe2vector (li, lc, lcoeff);
108
109    std::vector <quadElem> qcoeff;
110    indcoe2vector (qi, qj, qc, qcoeff);
111
112    if (!nq) ret = new exprGroup (c0, lcoeff); // create exprGroup
113    else     ret = new exprQuad  (c0, lcoeff, qcoeff); // create exprQuad
114  }
115
116  delete [] li;
117  delete [] lc;
118  delete [] qi;
119  delete [] qj;
120  delete [] qc;
121
122#ifdef DEBUG
123  printf ("\nlinstand ==> "); 
124  ret -> print (); printf ("\n"); 
125  //  ret -> Image () -> print (); printf ("\n");
126#endif
127
128  //if (ret -> Type () == AUX)
129  //return ret;
130
131  return (addAux ? addAuxiliary (ret) : new exprAux (ret, &domain_));
132}
Note: See TracBrowser for help on using the repository browser.