C++ LinearFunctionPtr类代码示例

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;
}

VariablePtr Transformer::newVar_(LinearFunctionPtr lf, double d,
                                 ProblemPtr newp)
{
  VariablePtr ov;
  FunctionPtr f;
  ConstraintPtr cnew;

  ov = yLfs_->findY(lf, d);
  if (!ov) {
    LinearFunctionPtr lf2 = lf->clone();
    ov = newp->newVariable(VarTran);
    yLfs_->insert(ov, lf2, d);
    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;
}

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;
  }
}

VariablePtr Transformer::newVar_(CGraphPtr cg, ProblemPtr newp)
{
  VariablePtr iv;
  VariablePtr ov = VariablePtr(); // NULL
  FunctionPtr f;
  LinearFunctionPtr lf;
  ConstraintPtr cnew;
  assert(cg);

  if (OpSumList!=cg->getOut()->getOp()) {
    ov = yUniExprs_->findY(cg);
  } 

  if (!ov) {
    ov = newp->newVariable(VarTran);
    lf = (LinearFunctionPtr) new LinearFunction();
    lf->addTerm(ov, -1.0);
    f = (FunctionPtr) new Function(lf, cg);
    cnew = newp->newConstraint(f, 0.0, 0.0);
#if SPEW
      logger_->msgStream(LogDebug) << me_ << "added new constraint "
                                   << std::endl;
      cnew->write(logger_->msgStream(LogDebug));
#endif 
    assignHandler_(cg, cnew);

    yUniExprs_->insert(ov, cg);
  }
  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 
Hypergraph::create(std::map<ConstVariablePtr, SetOfVars > const &terms)
{

  // Add the vertices
  std::map<ConstVariablePtr, SetOfVars >::const_iterator terms_it;
  for(terms_it = terms.begin(); terms_it != terms.end(); ++terms_it) {
    SetOfVars const &jt = terms_it->second;
    for(SetOfVars::iterator it = jt.begin(); it != jt.end(); ++it) {
      V_.insert(*it);
    }
  }

  // Add empty lists for the hyperedges
  for(SetOfVars::iterator v_it = V_.begin(); v_it != V_.end(); ++v_it) {
    std::list<SetOfVars> e_list;
    adjList_.insert(make_pair(*v_it, e_list));
  }

  // Now add the edges
  for(terms_it = terms.begin(); terms_it != terms.end(); ++terms_it) {
    SetOfVars const &jt = terms_it->second;
    E_.insert(jt);
    // Put it in the adjacency lists
    for (SetOfVars::const_iterator it = jt.begin(); it != jt.end(); ++it) {
      AdjListType::iterator pos = adjList_.find(*it);
      assert(pos != adjList_.end());
      std::list<SetOfVars> &e_list = pos->second;
      e_list.push_back(jt);
    }
    
    // Determine weight
    ConstVariablePtr zvar = terms_it->first;
    double zweight = 0.0;
    for (ConstrSet::iterator it2 = zvar->consBegin(); it2 != zvar->consEnd(); ++it2) {
      const LinearFunctionPtr lf = (*it2)->getLinearFunction();
      double w = lf->getWeight(zvar);
#if defined(DEBUG_MULTILINEARTERMS_HANDLER)
      zvar->write(std::cout);
      std::cout << "  has weight: " << w << " in constraint: ";
      (*it2)->write(std::cout);
#endif
      zweight += fabs(w);
    }

    // Add objective weight
    const LinearFunctionPtr obj = problem_->getObjective()->getLinearFunction();
    if (obj != 0) {
      double w = obj->getWeight(zvar);
      zweight += fabs(w);
    }

    weights_.insert(make_pair(jt, zweight));
    originalWeights_.insert(make_pair(jt, zweight));
  }
}

void ParQGHandler::addCut_(const double *nlpx, const double *lpx,
                        ConstraintPtr con, CutManager *cutman,
                        SeparationStatus *status)
{
  int error=0;
  ConstraintPtr newcon;
  std::stringstream sstm;
  double c, lpvio, act, cUb;
  FunctionPtr f = con->getFunction();
  LinearFunctionPtr lf = LinearFunctionPtr();

  act = con->getActivity(nlpx, &error);
  if (error == 0) {
    linearAt_(f, act, nlpx, &c, &lf, &error);
    if (error==0) {
      cUb = con->getUb();
      lpvio = std::max(lf->eval(lpx)-cUb+c, 0.0);
      if ((lpvio>solAbsTol_) &&
          ((cUb-c)==0 || (lpvio>fabs(cUb-c)*solRelTol_))) {
#if SPEW
        logger_->msgStream(LogDebug) << me_ << " linearization of constraint "
          << con->getName() << " violated at LP solution with violation = " <<
          lpvio << ". OA cut added." << std::endl;
#endif
        ++(stats_->cuts);
        sstm << "_OAcut_";
        sstm << stats_->cuts;
        *status = SepaResolve;
        f = (FunctionPtr) new Function(lf);
        newcon = rel_->newConstraint(f, -INFINITY, cUb-c, sstm.str());
        CutPtr cut = (CutPtr) new Cut(minlp_->getNumVars(),f, -INFINITY,
                                      cUb-c, false,false);
        cut->setCons(newcon);
        cutman->addCutToPool(cut);
        return;
      }
    }
  }	else {
    logger_->msgStream(LogError) << me_ << "Constraint not defined at"
      << " this point. "<<  std::endl;
#if SPEW
          logger_->msgStream(LogDebug) << me_ << "constraint " <<
            con->getName() << " not defined at this point." << std::endl;
#endif
  }
  return;
}

bool SOS1Handler::isGUB(SOS *sos)
{
  VariablePtr var;
  int nz = sos->getNz();
  ConstConstraintPtr c;
  FunctionPtr f;
  LinearFunctionPtr lf;
  bool isgub = false;


  for (VariableConstIterator viter = sos->varsBegin(); viter!=sos->varsEnd();
       ++viter) {
    var = *viter;
    if (Binary!=var->getType() || ImplBin!=var->getType()) {
      return false;
    }
  }

  // If we are here, then it means all variables are binary. Check if their sum
  // is <= 1.
  var = *(sos->varsBegin());
  for (ConstrSet::iterator it=var->consBegin();
       it!=var->consEnd() && false==isgub; ++it) {
    c = *it;
    f = c->getFunction();
    isgub = true;
    if (Linear!=f->getType()) {
      isgub = false;
    } else if ((UInt) nz==f->getNumVars()) {
      isgub = false;
    } else if (fabs(c->getUb()-1.0)>zTol_) {
      isgub = false;
    } else {
      lf = f->getLinearFunction();
      for (VariableConstIterator viter = sos->varsBegin();
           viter!=sos->varsEnd(); ++viter) {
        if ((lf->getWeight(*viter)-1.0)>zTol_) {
          isgub = false;
          break;
        }
      }
    }
  }

  return isgub;
}

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
MultilinearTermsHandler::setupWeights_()
{
  for(ConstTermIterator it = termsR_.begin(); it != termsR_.end(); ++it) {
    ConstVariablePtr zvar = it->first;
    double zweight = 0.0;
    for (ConstrSet::iterator it2 = zvar->consBegin(); it2 != zvar->consEnd(); ++it2) {
      const LinearFunctionPtr lf = (*it2)->getLinearFunction();
      double w = lf->getWeight(zvar);
#if defined(DEBUG_MULTILINEARTERMS_HANDLER)
      zvar->write(std::cout);
      std::cout << "  has weight: " << w << " in constraint: ";
      (*it2)->write(std::cout);
#endif
      zweight += fabs(w);
    }
    weights_.insert(std::make_pair(zvar, zweight));
  }
  
}

void Transformer::assignHandler_(CGraphPtr cg, ConstraintPtr c)
{
  switch (cg->getOut()->getOp()) {
  case OpMult:
  case OpSqr:
    qHandler_->addConstraint(c);
    break;
  default:
    {
    VariablePtr iv;
    VariablePtr ov = VariablePtr();
    LinearFunctionPtr lf = c->getFunction()->getLinearFunction();
    if (lf) {
      assert(lf->getNumTerms()==1);
      ov = lf->termsBegin()->first;
    }
    iv = *(c->getFunction()->getNonlinearFunction()->varsBegin());
    uHandler_->addConstraint(c, iv, ov, 'E');
    }
  }
}

void PolynomialFunction::removeLinear(LinearFunctionPtr lf)
{
  for (MonomialIter it = terms_.begin(); it!=terms_.end(); ) {
    if ((*it)->getDegree() == 1) {
      VarIntMapConstIterator it2 = (*it)->termsBegin();
      lf->incTerm(it2->first, (*it)->getCoeff());
      it = terms_.erase(it);
    } else {
      ++it;
    }
  }
}

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;

}

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
//.........这里部分代码省略.........

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);
//.........这里部分代码省略.........