C++ LinearFunctionPtr::addTerm方法代码示例

VariablePtr Transformer::newVar_(VariablePtr iv, double d, ProblemPtr newp)
{
  if (fabs(d)<zTol_) {
    return iv;
  } else {
    LinearFunctionPtr lf = (LinearFunctionPtr) new LinearFunction();
    VariablePtr ov;
    FunctionPtr f;
    ConstraintPtr cnew;

    ov = yVars_->findY(iv, d);
    if (!ov) {
      ov = newp->newVariable(VarTran);
      yVars_->insert(ov, iv, d);
      lf->addTerm(iv, 1.0);
      lf->addTerm(ov, -1.0);
      f = (FunctionPtr) new Function(lf);
      cnew = newp->newConstraint(f, -d, -d);
#if SPEW
      logger_->msgStream(LogDebug) << me_ << "added new constraint "
                                   << std::endl;
      cnew->write(logger_->msgStream(LogDebug));
#endif 
      lHandler_->addConstraint(cnew);
    }
    return ov;
  }
}

// This chunk of code adds
// z_t = \sum_{k=1}^{2 |V_g|} \lambda_k^g, \Prod_{j \in J_t} \chi_j^{g,k}
//   \forall J_t \subseteq V_g
void
MultilinearTermsHandler::handleZDefConstraints_(RelaxationPtr relaxation, HandleCallingFunction wherefrom, ModVector &mods)
{
  for(ConstTermIterator it = termsR_.begin(); it != termsR_.end(); ++it) {
    ConstVariablePtr zt = it->first;
    SetOfVars const &jt = it->second;

    for (UInt gix = 0; gix < groups_.size(); ++gix) {
      SetOfVars &vg = groups_[gix];
       
      if (std::includes(vg.begin(), vg.end(), jt.begin(), jt.end())) { 
        // jt is a subset of vg, add constraint
        LinearFunctionPtr lf = (LinearFunctionPtr) new LinearFunction();
        lf->addTerm(zt, -1.0);
        int pix = 0;

        for (std::set<SetOfVars>::iterator it2 = points_[gix].begin(); 
             it2 != points_[gix].end(); ++it2) {
          double prodval = 1.0;
          VariablePtr lam = lambdavars_[gix][pix];

          for(SetOfVars::const_iterator jt_it = jt.begin(); jt_it != jt.end(); ++jt_it) {
            ConstVariablePtr xvar = *jt_it;
            double tmp = varIsAtLowerBoundAtPoint_(xvar, *it2) ? xvar->getLb() : xvar->getUb();
            prodval *= tmp;
          }

          lf->addTerm(lam, prodval);
          ++pix;
        }        
        FunctionPtr f = (FunctionPtr) new Function(lf);

        if (wherefrom == relaxInit_Call) {
          ConstraintPtr c = relaxation->newConstraint(f, 0.0, 0.0);
          zConMap_.insert(std::make_pair(IntVarPtrPair(gix, zt), c));
        }
        else { 
          IntVarPtrPairConstraintMap::iterator pos;
          pos = zConMap_.find(IntVarPtrPair(gix, zt));
          if (pos == zConMap_.end()) {
            assert(0);
          }
          ConstraintPtr c = pos->second;
          //XXX Here you should just check if the constraint really was going to
          //change and do nothing if it doesn't...  (will be faster).
          if (wherefrom == relaxNodeInc_Call) {
            relaxation->changeConstraint(c, lf, 0.0, 0.0);
          }
          LinConModPtr lcmod = (LinConModPtr) new LinConMod(c, lf, 0.0, 0.0); 
          mods.push_back(lcmod);
        }
      }
      
    }
  } 


}

void ParQGHandler::linearizeObj_()
{
  ObjectivePtr o = minlp_->getObjective();
  FunctionType fType = o->getFunctionType();
  if (!o) {
    assert(!"need objective in QG!");
  } else if (fType != Linear && fType != Constant) {
    oNl_ = true;
    FunctionPtr f;
    VariablePtr vPtr;
    ObjectiveType objType = o->getObjectiveType();
    LinearFunctionPtr lf = (LinearFunctionPtr) new LinearFunction();
    for (VariableConstIterator viter=rel_->varsBegin(); viter!=rel_->varsEnd();
         ++viter) {
      vPtr = *viter;
      if (vPtr->getName() == "eta") {
        assert(o->getObjectiveType()==Minimize);
        rel_->removeObjective();
        lf->addTerm(vPtr, 1.0);
        f = (FunctionPtr) new Function(lf);
        rel_->newObjective(f, 0.0, objType);
        objVar_ = vPtr;
        break;
      }
    }
  }
  return;
}

void Transformer::makeObjLin_()
{
  ObjectivePtr obj;
  assert(p_);
  assert(newp_);

  obj = p_->getObjective();
  if (!obj) {
    return;
  }
  
  if (obj->getFunctionType() != Linear && obj->getFunctionType() != Constant) {
    VariablePtr eta = newp_->newVariable(VarTran);
    FunctionPtr etaf;
    FunctionPtr f = obj->getFunction();
    LinearFunctionPtr lz = LinearFunctionPtr(new LinearFunction());
    LinearFunctionPtr lz2;
    ObjectiveType otype = obj->getObjectiveType();
    ConstraintPtr objcon;

    lz->addTerm(eta, -1.0);
    f->add(lz);
    
    if (otype == Minimize) {
      //XXX Do you want to keep track of the objective constraint?
      objcon = newp_->newConstraint(f, -INFINITY, 0.0);
      if (lHandler_) {
        lHandler_->addConstraint(objcon);
      }
    } else {
      objcon = newp_->newConstraint(f, 0.0, INFINITY);
      if (lHandler_) {
        lHandler_->addConstraint(objcon);
      }
    }

    lz2 = (LinearFunctionPtr) new LinearFunction();
    lz2->addTerm(eta, 1.0);
    etaf = (FunctionPtr) new Function(lz2);

    newp_->newObjective(etaf, obj->getConstant(), otype);
  }
}

void CxUnivarConstraintData::addLin(RelaxationPtr rel, ConstVariablePtr riv,
				    ConstVariablePtr rov, FunctionPtr fn,
                                    DoubleVector& tmpX, DoubleVector& grad,
                                    bool init, ModVector &mods)
{
 
  int error;
  ConstraintPtr cons; 
  double xlb = riv->getLb();
  double xub = riv->getUb();
  double fxlbval=0, fxubval=0, dfxlbval=0, dfxubval=0;
  double tmpxval, fxval, dfxval; 
  LinearFunctionPtr lf; 
  FunctionPtr f;

  // More sophisticated strategies hopefully could be obtained by simply
  // changing this array 
  int npts = 3;
  double xvals[] = {xlb, xub, (xub-xlb)/2.0};

#if defined(DEBUG_CXUNIVARHANDLER)
  std::cout << "Adding linearizations.  rix id: " << riv->getId() 
	    << " rix index: " << riv->getIndex() << " rov id: " << rov->getId() 
	    << " rov index: " << rov->getIndex()
	    << " xlb: " << xlb << " xub: " << xub << std::endl;
#endif
  
  for (int i = 0; i < npts; i++) {

    // Zero out tmpX and grad each time, or else bad things happen
    for (UInt j = 0; j < tmpX.size(); ++j) {
      tmpX[j] = 0.0;
      grad[j] = 0.0;
    }
    
    if (i == 2) {
      // Third linearization point taken to be where first two intersect:
      // x3 = (f'(xub)*xub - f'(xlb)*xlb + f(xlb) - f(xub))/(f'(xub) - f'(xlb))
      // Unless this would put it too close to one of the end points
      if (dfxubval - dfxlbval > 0.0001 || dfxubval - dfxlbval < -0.0001) {
        tmpxval = (dfxubval*xub - dfxlbval*xlb + fxlbval - fxubval)/
                  (dfxubval - dfxlbval);
        if (tmpxval < xlb + (xub-xlb)*0.05) {
          xvals[2] = xlb + (xub-xlb)*0.05;
        }
        else if (tmpxval > xub - (xub-xlb)*0.05) {
          xvals[2] = xub - (xub-xlb)*0.05;
        }
        else {
          xvals[2] = tmpxval;
        }
      }
    }
    tmpX[riv->getIndex()] = xvals[i];
    error = 0;
    fxval =  fn->eval(tmpX, &error);
    fn->evalGradient(&tmpX[0], &grad[0], &error);
#if defined(DEBUG_CXUNIVARHANDLER2)
    for (UInt j = 0; j < tmpX.size(); ++j) {
      std::cout << "x[" << j << "] = " << tmpX[j] << " dfdx[" << j << "] = "
                << grad[j] << std::endl;
    }
#endif
    dfxval = grad[riv->getIndex()];
    if (i == 0) {
       fxlbval = fxval;
       dfxlbval = dfxval; 
    }
    else if (i == 1) {
       fxubval = fxval;
       dfxubval = dfxval; 
    }
    // linearization:  rov >= f(xval) + f'(xval)(riv - xval) 
    //                 rov - f'(xval)*riv >= f(xval) - f'(xval)*xval
    lf = (LinearFunctionPtr) new LinearFunction();
    lf->addTerm(rov, 1.0);
    lf->addTerm(riv, -dfxval);
    if (init) {
        f = (FunctionPtr) new Function(lf);
    	cons = rel->newConstraint(f, fxval - dfxval*xvals[i], INFINITY);
	linCons_.push_back(cons);
    }
    else {
#if defined(DEBUG_CXUNIVARHANDLER)
       std::cout << "Will change 'linearization  ' constraint to have "
                 << "linear function: ";
       lf->write(std::cout);  
       std::cout << std::endl;
#endif

       rel->changeConstraint(linCons_[i], lf, fxval - dfxval*xvals[i], INFINITY); 
       LinConModPtr lcmod = (LinConModPtr) new LinConMod(linCons_[i], lf, 
                                                         fxval -
                                                         dfxval*xvals[i],
                                                         INFINITY); 
       mods.push_back(lcmod);
    }
  }
  tmpX[riv->getIndex()] = 0.0;
  grad[riv->getIndex()] = 0.0;
//.........这里部分代码省略.........

void
MultilinearFormulation::makeTermByTerm()
{

  // First we do some processing of the instance to determine how many 
  //  multilinear or quadratic constraints, and we store the indicies
  // in the instance for later    
  vector<int> lcid;
  vector<int> mlcid;
  for(ConstConstraintIterator it = originalInstance_->consBegin(); 
      it != originalInstance_->consEnd(); ++it) {
    FunctionType ft = (*it)->getFunctionType();
    if (ft == Multilinear) {
      mlcid.push_back((*it)->getId());
    }    
    else if (ft == Bilinear) {
      mlcid.push_back((*it)->getId());
    }
    else if (ft == Linear) {
      lcid.push_back((*it)->getId());
    }
  }

  // add x variables
  vector<double> lb;
  vector<double> ub;

  vector<VariablePtr> xvars;
  int nv = 0;
  for (ConstVariableIterator it = originalInstance_->varsBegin(); 
       it != originalInstance_->varsEnd(); ++it) {
    VariablePtr v = *it;
    VariablePtr vnew = VariablePtr(new Variable(nv, v->getLb(), v->getUb(), v->getType()));
    lb.push_back(v->getLb());
    ub.push_back(v->getUb());

    variableMapping_.insert(make_pair(vnew,v));    
    variables_.push_back(vnew);
    xvars.push_back(vnew);
    nv++;
  }

  // Add the linear constraints
  for(int i = 0; i < lcid.size(); i++) {
    const ConstraintPtr mlc = originalInstance_->getConstraint(lcid[i]);
    const LinearFunctionPtr olf = mlc->getLinearFunction();
    LinearFunctionPtr lf = LinearFunctionPtr(new LinearFunction());
    for(ConstVariableGroupIterator it = olf->varsBegin(); it != olf->varsEnd(); ++it) {
      lf->addTerm(xvars[it->first->getId()], it->second);
    }
    FunctionPtr f = (FunctionPtr) new Function(lf);
    ConstraintPtr c = (ConstraintPtr) new Constraint(f, mlc->getLb(), mlc->getUb());
    constraints_.push_back(c);    

#if defined(DEBUG_TERM_BY_TERM)
      c->display();
#endif

  }


  // The w variables
  vector<VariablePtr> wvars;
  // This holds a map between the 'w' variable added and indices of x vars in multilinear product
  map <VariablePtr, vector<int> > mlterms;

  // Go through multilinear rows.  Add constraints, and create maps
  for(int i = 0; i < mlcid.size(); i++) {

    const ConstraintPtr omlc = originalInstance_->getConstraint(mlcid[i]);
    const LinearFunctionPtr olf = omlc->getLinearFunction();
    const QuadraticFunctionPtr oqf = omlc->getQuadraticFunction();
    const NonlinearFunctionPtr onlf = omlc->getNonlinearFunction();
    //!!! Don't make this shared by boost, it will get confused in counting
    MultilinearFunction *omlf = dynamic_cast<MultilinearFunction *>(onlf.get());
    
    LinearFunctionPtr lf = LinearFunctionPtr(new LinearFunction());

    // Linear part of constraint remains the same
    for(ConstVariableGroupIterator it = olf->varsBegin(); it != olf->varsEnd(); ++it) {
      lf->addTerm(xvars[it->first->getId()], it->second);
    }

    // Quadratic part gets a new variable for every term
    for(ConstVariablePairGroupIterator it = oqf->begin(); it != oqf->end(); ++it) {

      vector<int> mlix;
      mlix.push_back(it->first.first->getId());
      mlix.push_back(it->first.second->getId());

      VariablePtr w = VariablePtr(new Variable(nv, -INFINITY,INFINITY, Continuous));
      nv++;
      variables_.push_back(w);
      wvars.push_back(w);
      //XXX Need to store term for evaluation
      mlterms.insert(make_pair(w,mlix));
      lf->addTerm(w,it->second);
    }

    // Multilinear part gets a new var for every term
//.........这里部分代码省略.........

void
MultilinearTermsHandler::relaxInitInc(RelaxationPtr relaxation, bool *)
{

  //  General notes...
  // Use Id to order variables
  // x[v->getIndex()] returns solution.
  
  /*
    0) Create "terms" map using the relaxation variables, not the original
       problem variables
    1) Create groups based on container of terms
    2) Add all lambda variables based on groups.
       (Be sure to keep a container of which constraints are associated with
       which groups)
   */

#if defined(DEBUG_MULTILINEARTERMS_HANDLER)
  cout << "In MultilinearTermHandler::relaxInitInc()" << endl << " Current Relaxation: " << endl;
  relaxation->write(cout);
  cout << "And current terms: " << endl;

  for(ConstTermIterator it = termsO_.begin(); it != termsO_.end(); ++it) {
    std::cout << "zvar: ";
    it->first->write(std::cout);
    std::cout << "Contains vars: " << std::endl;
    for(SetOfVars::const_iterator it2 = it->second.begin(); 
        it2 != it->second.end(); ++it2) {
      (*it2)->write(cout);
    }
  }
#endif  

  // Copy the original variable pointers into a data structure consisting of 
  //  relaxation variables   (It just makes likfe easier to keep things in 
  //  terms of relaxation variables).
  for (ConstTermIterator it = termsO_.begin(); it != termsO_.end(); ++it) {
    ConstVariablePtr ov = it->first;
    SetOfVars const &ovset = it->second;
    ConstVariablePtr rv = relaxation->getRelaxationVar(ov);
    SetOfVars rvset;
    for(SetOfVars::const_iterator it2 = ovset.begin(); it2 != ovset.end(); ++it2) {
      rvset.insert( relaxation->getRelaxationVar(*it2));
    }
    termsR_.insert(make_pair(rv, rvset));
  }
      
    
  // First we make the groups.

  makeGroups_();

#if defined(DEBUG_MULTILINEARTERMS_HANDLER)
  cout << "After making groups: " << endl;

  for (ConstGroupIterator it = groups_.begin(); it != groups_.end(); ++it) {  
    cout << "Group of: " << endl;
    for(set<ConstVariablePtr>::const_iterator it2 = it->begin(); it2 != it->end(); ++it2) {
      (*it2)->write(cout);
    }
  }

  // Just checking groups now.  Stop here.
  //exit(1);
#endif

  /* This chunk of code creates the (powerset) representation 
   * of the extreme points in the groups, adds the lambda variables,
   * and adds the convexity constraints
   */

  lambdavars_.resize(groups_.size());
  for(UInt gix = 0; gix < groups_.size(); ++gix) {    

    LinearFunctionPtr lf = (LinearFunctionPtr) new LinearFunction();

    std::set<SetOfVars> p = powerset_(groups_[gix]);
    int pix = 0;
    for (std::set<SetOfVars>::iterator it2 = p.begin(); it2 != p.end(); ++it2) {
      std::string name;
      std::stringstream name_stream;
      name_stream << "lam_" << gix << "_" << pix;
      name = name_stream.str();
      
      VariablePtr lam;
      lam = relaxation->newVariable(0.0, 1.0, Continuous, name);
      lf->addTerm(lam, 1.0);
#if defined(DEBUG_MULTILINEARTERMS_HANDLER2)
      std::cout << "Adding: " << name << std::endl;
#endif

      lambdavars_[gix].push_back(lam);
      pix++;
    }    
    points_.push_back(p);
    FunctionPtr f = (FunctionPtr) new Function(lf);
    relaxation->newConstraint(f, 1.0, 1.0);
  }  


//.........这里部分代码省略.........

ModificationPtr MultilinearTermsHandler::getBrMod(BrCandPtr cand, DoubleVector &xval, 
                                                  RelaxationPtr , BranchDirection dir)
{
  LinModsPtr linmods;

  //XXX Put (bool init) back in handle{x,z}def...

  BrVarCandPtr  vcand = boost::dynamic_pointer_cast <BrVarCand> (cand);
  VariablePtr v = vcand->getVar();
  
  double branching_value = xval[v->getIndex()];
  BoundType lu;
  VariableType vtype = v->getType();

  // Change bounds on the x var (called v here)
  if (dir == DownBranch) { 
    lu = Upper;    
    if (vtype != Continuous)  branching_value = floor(branching_value);
  }
  else {
    lu = Lower;
    if (vtype != Continuous)  branching_value = ceil(branching_value);
  }

  linmods = (LinModsPtr) new LinMods();

  VarBoundModPtr vmod = (VarBoundModPtr) new VarBoundMod(v, lu, branching_value);
  linmods->insert(vmod);
  
  
  // This chunk of code changes the
  // x_{V_g} = \sum_{k=1}^{2 |V_g|} \lambda_k^g \chi^{k,g} \forall g \in G


  for (UInt gix = 0; gix < groups_.size(); ++gix) {
    for(SetOfVars::const_iterator it = groups_[gix].begin(); it != groups_[gix].end(); ++it) {
      ConstVariablePtr xvar = *it;
      if (v != xvar) continue;

      LinearFunctionPtr lf = (LinearFunctionPtr) new LinearFunction();
      lf->addTerm(xvar, -1.0);

      UInt pix = 0;
      for (std::set<SetOfVars>::iterator it2 = points_[gix].begin(); it2 != points_[gix].end(); ++it2) {
        VariablePtr lam = lambdavars_[gix][pix];
        double val = -INFINITY;

        bool atLower = varIsAtLowerBoundAtPoint_(v, *it2);
        bool atUpper = !atLower;
        
        if (lu == Upper && atUpper) val = branching_value;
        else if (lu == Lower && atLower) val = branching_value;
        else val = (atLower ? v->getLb() : v->getUb());

        lf->addTerm(lam, val);
        ++pix;
      }       
      FunctionPtr f = (FunctionPtr) new Function(lf);

      IntVarPtrPairConstraintMap::iterator pos;
      pos = xConMap_.find(IntVarPtrPair(gix, xvar));
      if (pos == xConMap_.end()) {
        assert(0);
      }
      ConstraintPtr c = pos->second;

      LinConModPtr lcmod = (LinConModPtr) new LinConMod(c, lf, 0.0, 0.0);
#if defined(DEBUG_MULTILINEARTERMS_HANDLER)
      std::cout << "getBrMod().  Will change 'x =' constraint to have linear function ";
      lf->write(std::cout);        
      std::cout << std::endl;
#endif
      linmods->insert(lcmod);
    
    }
  }

  // This will change the z_t = sum \sum_{k=1}^{2|V_g} \lambda_k^g \chi^{k,g}.
  //  Probably not very efficient way to do this...
  for(ConstTermIterator it = termsR_.begin(); it != termsR_.end(); ++it) {
    SetOfVars const &jt = it->second;

    for (UInt gix = 0; gix < groups_.size(); ++gix) {
      SetOfVars &vg = groups_[gix];

      std::set<ConstVariablePtr>::iterator pos1;
      pos1 = jt.find(v);
      if (pos1 == jt.end()) continue; // J_t does not contain v, go to next group
      // J_t is not in V_g, go to next group
      if (! std::includes(vg.begin(), vg.end(), jt.begin(), jt.end())) continue;

      ConstVariablePtr zvar = it->first;
      LinearFunctionPtr lf = (LinearFunctionPtr) new LinearFunction();
      lf->addTerm(zvar, -1.0);

      // Get ConstraintToChange
      IntVarPtrPairConstraintMap::iterator pos2;
      pos2 = zConMap_.find(IntVarPtrPair(gix, zvar));
      if (pos2 == zConMap_.end()) {
        assert(0);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		  if(p_mult || pcf_constant || pvf_constant) {
		    p_rltf = p_mult + pcf_constant + pvf_constant;
		    q_rltf = q_lcf_lvf + qvf_constant + qcf_constant;
		    l_rltf = lvf_constant + lcf_constant;
		    
		    NonlinearFunctionPtr n_rltf = boost::dynamic_pointer_cast <NonlinearFunction> (p_rltf);
		    rltf = (FunctionPtr) new Function(l_rltf, q_rltf, n_rltf);
		  }
		  else {
		    q_rltf = q_lcf_lvf + qcf_constant + qvf_constant;
		    l_rltf = lcf_constant + lvf_constant;
		    rltf = (FunctionPtr) new Function(l_rltf, q_rltf);
		  }
		  
		  temp_v.push_back(rltf);
		  temp_v_lo.push_back(-1 * (*con_it)->getUb() * v_up[v_cntr]);
		  temp_v_up.push_back(INFINITY);
		}
	      }
	    }
	  }
	}
      }
    }

    // multiply variable bound constraints against regular constraints
   
    for(VariableConstIterator var_it = inst->varsBegin(); var_it != inst->varsEnd(); ++var_it) {
      v_cntr = -1;

      // If the variable has a lower bound
      if((*var_it)->getLb() > -1e-10) {
	LinearFunctionPtr lbf = LinearFunctionPtr(new LinearFunction());
	lbf->addTerm(*var_it, 1);

	// Iterate through the functions in v
	for(std::vector<FunctionPtr>::iterator vr_it = v.begin(); vr_it != v.end(); ++vr_it) {
	  v_cntr++;
	  
	  // If v is = then skip it for now 
	  if(v_up[v_cntr] - v_lo[v_cntr] > 1e-6) {
	    
	    FunctionPtr vf = *vr_it;

	    // Break down the v function
	    LinearFunctionPtr lvf = vf->getLinearFunction();
	    QuadraticFunctionPtr qvf = vf->getQuadraticFunction();
	    NonlinearFunctionPtr nvf = vf->getNonlinearFunction();
	    PolyFunPtr pvf;
	    if(nvf)
	      pvf = boost::dynamic_pointer_cast <PolynomialFunction> (nvf);
	    
	    // Multiply parts of the two functions against eachother
	    QuadraticFunctionPtr q_lbf_lvf;
	    PolyFunPtr p_lbf_qvf;
	    PolyFunPtr p_lbf_pvf;
	    
	    // Multiply the linear part of lbf to the linear part of vf
	    if(lbf && lvf)
	      q_lbf_lvf = lbf * lvf;
	    
	    // Multiply the linear part of lbf to the quadratic part of vf
	    if(lbf && qvf)
	      p_lbf_qvf = qvf * lbf;
	    
	    // Multiply the linear part of lbf to the polynomial part of vf

void CxUnivarConstraintData::addSecant(RelaxationPtr rel,
                                       ConstVariablePtr riv,
                                       ConstVariablePtr rov,
                                       FunctionPtr fn, DoubleVector& tmpX,
                                       bool init, ModVector &mods) 
{

  int error;
  double xlb, xub, fxlb, fxub, m, intercept;
  LinearFunctionPtr lf; 
  FunctionPtr f;

  // First add the secant inequalities based on variable bounds
  xlb = riv->getLb();
  xub = riv->getUb();
	
#if defined(DEBUG_CXUNIVARHANDLER)
  std::cout << "Adding secant on variable rix index: " << riv->getIndex() 
	    << " rov index: " << rov->getIndex()
	    << " xlb: " << xlb << " xub: " << xub << std::endl;
#endif 
  // no secant if unbounded either way
  if (xlb <= -0.9*INFINITY || xub >= 0.9*INFINITY) {
    std::cout << "Cannot add secant -- bound is infinite" << std::endl;
    return;
  }

  // TODO: Check the error value!
  tmpX[riv->getIndex()] = xlb;
  fxlb =  fn->eval(tmpX, &error);
  tmpX[riv->getIndex()] = xub;
  fxub =  fn->eval(tmpX, &error);
  tmpX[riv->getIndex()] = 0.0;

  // TODO: check/remedy numerical issues in this division
  if (xub - xlb > 10e-7) {
    m = (fxub - fxlb)/(xub - xlb);
  }
  else {
    m = 0.0;
  }

  intercept = fxlb - m*xlb;
  lf = (LinearFunctionPtr) new LinearFunction();
  lf->addTerm(rov, 1.0);
  lf->addTerm(riv, -m);

  // rovar <= m*rivar + intercept 
  if (init) {
     f = (FunctionPtr) new Function(lf);
     secCon_ = rel->newConstraint(f, -INFINITY, intercept);
  }
  else {
    rel->changeConstraint(secCon_, lf, -INFINITY, intercept);
    LinConModPtr lcmod = (LinConModPtr) new LinConMod(secCon_, lf, -INFINITY,
                                                      intercept);
    mods.push_back(lcmod);
  }

  
}

void ParQGHandler::oaCutToObj_(const double *nlpx, const double *lpx,
                            CutManager *, SeparationStatus *status)
{
  if (oNl_) {
    int error=0;
    FunctionPtr f;
    double c, vio, act;
    ConstraintPtr newcon;
    std::stringstream sstm;
    ObjectivePtr o = minlp_->getObjective();

    act = o->eval(lpx, &error);
    if (error == 0) {
      vio = std::max(act-relobj_, 0.0);
      if ((vio > solAbsTol_) &&
          (relobj_ == 0 || vio > fabs(relobj_)*solRelTol_)) {
#if SPEW
        logger_->msgStream(LogDebug) << me_ << " objective violated at LP "
          << " solution with violation = " << vio << std::endl;
#endif
        act = o->eval(nlpx, &error);
        if (error == 0) {
          f = o->getFunction();
          LinearFunctionPtr lf = LinearFunctionPtr();
          linearAt_(f, act, nlpx, &c, &lf, &error);
          if (error == 0) {
            vio = std::max(c+lf->eval(lpx)-relobj_, 0.0);
            if ((vio > solAbsTol_) && ((relobj_-c)==0
                                     || vio > fabs(relobj_-c)*solRelTol_)) {
#if SPEW
              logger_->msgStream(LogDebug) << me_ << "linearization of "
                "objective violated at LP solution with violation = " <<
                vio << ". OA cut added." << std::endl;
#endif
              ++(stats_->cuts);
              sstm << "_OAObjcut_";
              sstm << stats_->cuts;
              lf->addTerm(objVar_, -1.0);
              *status = SepaResolve;
              f = (FunctionPtr) new Function(lf);
              newcon = rel_->newConstraint(f, -INFINITY, -1.0*c, sstm.str());
            }
          }
        }
      }  else {
#if SPEW
        logger_->msgStream(LogDebug) << me_ << " objective feasible at LP "
          << " solution. No OA cut to be added." << std::endl;
#endif
      }
    }	else {
      logger_->msgStream(LogError) << me_
        << " objective not defined at this solution point." << std::endl;
#if SPEW
      logger_->msgStream(LogDebug) << me_ << " objective not defined at this "
        << " point." << std::endl;
#endif
    }
  }
  return;
}

void ParQGHandler::addInitLinearX_(const double *x)
{
  int error=0;
  FunctionPtr f;
  double c, act, cUb;
  std::stringstream sstm;
  ConstraintPtr con, newcon;
  LinearFunctionPtr lf = LinearFunctionPtr();

  for (CCIter it=nlCons_.begin(); it!=nlCons_.end(); ++it) {
    con = *it;
    act = con->getActivity(x, &error);
    if (error == 0) {
      f = con->getFunction();
      linearAt_(f, act, x, &c, &lf, &error);
      if (error == 0) {
        cUb = con->getUb();
#if USE_OPENMP
        sstm << "Thr_" << omp_get_thread_num();
#endif
        ++(stats_->cuts);
        sstm << "_OAcut_";
        sstm << stats_->cuts;
        sstm << "_AtRoot";
        f = (FunctionPtr) new Function(lf);
        newcon = rel_->newConstraint(f, -INFINITY, cUb-c, sstm.str());
        sstm.str("");
      }
    }	else {
      logger_->msgStream(LogError) << me_ << "Constraint" <<  con->getName() <<
        " is not defined at this point." << std::endl;
#if SPEW
      logger_->msgStream(LogDebug) << me_ << "constraint " <<
        con->getName() << " is not defined at this point." << std::endl;
#endif
    }
  }

  if (oNl_) {
    error = 0;
    ObjectivePtr o = minlp_->getObjective();
    act = o->eval(x, &error);
    if (error==0) {
      ++(stats_->cuts);
      sstm << "_OAObjcut_";
      sstm << stats_->cuts;
      sstm << "_AtRoot";
      f = o->getFunction();
      linearAt_(f, act, x, &c, &lf, &error);
      if (error == 0) {
        lf->addTerm(objVar_, -1.0);
        f = (FunctionPtr) new Function(lf);
        newcon = rel_->newConstraint(f, -INFINITY, -1.0*c, sstm.str());
      }
    }	else {
      logger_->msgStream(LogError) << me_ <<
        "Objective not defined at this point." << std::endl;
#if SPEW
      logger_->msgStream(LogDebug) << me_ <<
        "Objective not defined at this point." << std::endl;
#endif
    }
  }
  return;
}

BranchPtr CxUnivarHandler::doBranch_(BranchDirection UpOrDown,
                                     ConstVariablePtr v, double bvalue)
{
  BranchPtr branch;
  LinModsPtr linmods;

#if defined(DEBUG_CXUNIVARHANDLER)
  std::cout << "CxUnivarHandler, Branching: " << (UpOrDown == DownBranch ? "Down" : "Up")
            << " at value: " << bvalue << " on: " << std::endl;
  v->write(std::cout);
#endif

  // Zero out tmpX and grad each time, or else bad things happen
  for (UInt j = 0; j < tmpX_.size(); ++j) {
    tmpX_[j] = 0.0;
    grad_[j] = 0.0;
  }
    
  branch = (BranchPtr) new Branch();
  linmods = (LinModsPtr) new LinMods();

  // Change bounds on the x var (called v here)
  if (UpOrDown == DownBranch) {

    VarBoundModPtr mod = (VarBoundModPtr) new VarBoundMod(v, Upper, bvalue);
    linmods->insert(mod);

    // Find *all* cons_data that has v as an input variable.
    CxUnivarConstraintIterator dit; 
  
    for (dit = cons_data_.begin(); dit != cons_data_.end(); ++dit) {
      if ((*dit)->getRInputVar() == v) {

	ConstVariablePtr rov = (*dit)->getROutVar();
	FunctionPtr fn = (*dit)->getOriginalCon()->getFunction();
	int error;

	// Change the secant constraint
	ConstraintPtr secCon = (*dit)->getSecantCon();

	LinearFunctionPtr lf; 
	FunctionPtr f;

	double xlb = v->getLb();
	double xub = bvalue;

	// TODO: Check the error value!
	tmpX_[v->getIndex()] = xlb;
	double fxlb =  fn->eval(tmpX_, &error);
	tmpX_[v->getIndex()] = xub;
	double fxub =  fn->eval(tmpX_, &error);
	tmpX_[v->getIndex()] = 0.0;

	// TODO: check/remedy numerical issues in this division
	double m = 0.0;
	if (xub - xlb > 10e-7) {
	  m = (fxub - fxlb)/(xub - xlb);
	}
	double intercept = fxlb - m*xlb;
	lf = (LinearFunctionPtr) new LinearFunction();
	lf->addTerm(rov, 1.0);
	lf->addTerm(v, -m);

        LinConModPtr lcmod = (LinConModPtr) new LinConMod(secCon, lf,
                                                          -INFINITY,
                                                          intercept);
	linmods->insert(lcmod);

	// Change all linearization constraints
	ConstraintVector::iterator lin_it;
        for(lin_it = (*dit)->linConsBegin(); lin_it != (*dit)->linConsEnd();
            ++lin_it) {
	  ConstraintPtr c = *lin_it;
	}
      }
    }
  } else {
    VarBoundModPtr mod = (VarBoundModPtr) new VarBoundMod(v, Lower, bvalue);
    linmods->insert(mod);
  }

  assert(!"add Mod correctly here.");
  branch->addPMod(linmods);
  return branch;

}