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

	AllocaInsn::AllocaInsn(const VariablePtr& variable, const ValuePtr& size, const ValueList& dimensions) :
		Insn(IC_Stack, IT_Alloca), size(size), variable(variable), dimensions(dimensions) {
		assert(variable && "variable must not be null");
		assert(variable->getValueType() == Value::VT_Memory && "variable must not be a temporary");
		assert(size && "size must not be null");
		assert(size->getType()->isInt() && "size must be of type int");

		if (analysis::types::isArray(variable->getType())) {
			auto type = cast<ArrayType>(variable->getType());
			assert(type->getNumOfDimensions() == dimensions.size() && "dimension information mismatch with type");
		}
	}

Branches MultilinearTermsHandler::getBranches(BrCandPtr cand, DoubleVector & x,
                                         RelaxationPtr, SolutionPoolPtr)
{
  BrVarCandPtr vcand = boost::dynamic_pointer_cast <BrVarCand> (cand);
  VariablePtr v = vcand->getVar();
  double value = x[v->getIndex()];

  // can't branch on something that is at its bounds.
  if (!(value > v->getLb()+1e-8 && value < v->getUb()-1e-8)) {
    cerr << "Warning!  Branching on variable with bounds/value: [" << 
      v->getLb() << " , " << value << "  " << v->getUb() << " ]" << endl;
    //assert(value > v->getLb()+1e-8 && value < v->getUb()-1e-8);
  }

  Branches branches = (Branches) new BranchPtrVector();
  BranchPtr branch;
  branch = doBranch_(DownBranch, v, value);
  branches->push_back(branch);
  branch = doBranch_(UpBranch, v, value);
  branches->push_back(branch);


  logger_->msgStream(LogDebug2) << "branching on " << v->getName();
  logger_->msgStream(LogDebug2) << " <= " << value << " or " 
    << " >= " << value << endl;

#if defined(DEBUG_MULTILINEARTERMS_HANDLER)  
  cout << "branching on " << v->getName();
  cout << " <= " << value << " or " << " >= " << value << endl;
#endif

  return branches;
}

SecantMod::SecantMod(ConstraintPtr con, LinearFunctionPtr new_lf,
                     double new_rhs, VariablePtr x, BoundType lu, double new_b,
                     VariablePtr y)
{
  double y_lb, y_ub, b2;
  if (lu==Lower) {
    b2 = x->getUb();
    BoundsOnSquare(new_b, b2, y_lb, y_ub);
  } else {
    b2 = x->getLb();
    BoundsOnSquare(b2, new_b, y_lb, y_ub);
  }
  ymod_ = (VarBoundMod2Ptr) new VarBoundMod2(y, y_lb, y_ub);
  xmod_ = (VarBoundModPtr) new VarBoundMod(x, lu, new_b);
  lmod_ = (LinConModPtr) new LinConMod(con, new_lf, -INFINITY, new_rhs);
}

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 Sgd<T>::update_impl(const string &key, VariablePtr param) {
  Size_t size = param->size();

  const T *grad = param->get_grad_pointer<T>(this->ctx_);
  T *data = param->cast_data_and_get_pointer<T>(this->ctx_);
  std::transform(grad, grad + size, data, data,
                 [this](T g, T x) { return x - lr_ * g; });
}

void ParQGHandler::fixInts_(const double *x)
{
  double xval;
  VariablePtr v;
  VarBoundMod2 *m = 0;
  for (VariableConstIterator vit=minlp_->varsBegin(); vit!=minlp_->varsEnd();
       ++vit) {
    v = *vit;
    if (v->getType()==Binary || v->getType()==Integer) {
      xval = x[v->getIndex()];
      xval = floor(xval + 0.5);
      m = new VarBoundMod2(v, xval, xval);
      m->applyToProblem(minlp_);
      nlpMods_.push(m);
    }
  }
  return;
}

void LinFeasPump::constructObj_(ProblemPtr, ConstSolutionPtr)
{
  double value, lb, ub;
  VariablePtr variable;
  UInt i;
  //double obj_relaxation;
  double obj_weight;
  double constant;
  FunctionPtr f;
  LinearFunctionPtr olf_mod;

  //obj_relaxation = sol->getPrimal()[r_->getNumVars()-1];
  obj_weight     = 0; //std::min(1/(fabs(obj_relaxation)+1e-7), 0.1);
  obj_weight     = 0.0;
  i              = 0;
  constant       = 0;
  olf_mod        = olfClone_->clone();
  (*olf_mod)    *= obj_weight;

  for (VariableConstIterator v_iter=r_->varsBegin();
      v_iter!=r_->varsEnd(); ++v_iter, ++i) {
    variable = *v_iter;
    if (variable->getType() == Binary) {
      value  = roundedSol_[i];
      lb     = variable->getLb();
      ub     = variable->getUb();
      if (fabs(value - lb) > intTol_) { 
        olf_mod->addTerm(variable,-1.0);
        constant += ub;
      } else if (fabs(value - ub) > intTol_) {
        olf_mod->addTerm(variable,1.0);
        constant -= lb;
      } else {
        // add a new variable with coeff 1
        // add two constraints for absolute value
      }      
    }
  }
  f = (FunctionPtr) new Function(olf_mod);
  r_->changeObj(f, constant);
}

int VariableLuaMetaData::Index(lua_State * L, const VariablePtr var,
								const std::string & index)
{
	///*if (index == "Const")
	//{
	//	lua_pushboolean(L, var->GetType()->IsConst() ? 1 : 0);
	//	return 1;
	//}
	//else if (index == "ConstPointer")
	//{
	//	lua_pushboolean(L, var->GetType()->IsConstPointer() ? 1 : 0);
	//	return 1;
	//}
	//else if (index == "Pointer")
	//{
	//	lua_pushboolean(L, var->GetType()->IsPointer() ? 1 : 0);
	//	return 1;
	//}
	//else if (index == "Reference")
	//{
	//	lua_pushboolean(L, var->GetType()->IsReference() ? 1 : 0);
	//	return 1;
	//}
	//else if (index == "TypeNode")
	//{
	//	Model::NodeLuaMetaData::PutInstanceOnStack(L, var->GetType()->GetNode());
	//	return 1;
	//}*/
	if (index == "Initializer")
	{
		lua_pushstring(L, var->GetInitializer().c_str());
		return 1;
	}
	else if (index == "Type")
	{
		Model::TypeLuaMetaData::PutInstanceOnStack(L, var->GetType());
		return 1;
	}
	return 0;
}

// -----------------------------------------------------------------------------------
void CBCSolver::setUpVariablesAndObj(std::vector<VariablePtr>& variables, 
                                     ObjFuncPtr& obj) {
  for (int i = 0; i < variables.size(); i++) {
    VariablePtr v = variables.at(i);
    switch(v->type()) {
    case Variable::INT:
      builder_.setColumnBounds(i,intBound(v->lbound()),intBound(v->ubound()));
      builder_.setInteger(i);
      break;
    case Variable::LINEAR:
      builder_.setColumnBounds(i,doubleBound(v->lbound()),doubleBound(v->ubound()));
      break;
    }
    builder_.setObjective(i,obj->getModifier(v));
  }
}

// Implement Handler::getBranches().
Branches CxUnivarHandler::getBranches(BrCandPtr cand, DoubleVector &x,
                                      RelaxationPtr , SolutionPoolPtr )
{
  double minFromBds = 0.1;
  BrVarCandPtr vcand = boost::dynamic_pointer_cast <BrVarCand> (cand);
  VariablePtr v = vcand->getVar();

  double xval = x[v->getIndex()];
  double value = xval;  // Make sure branch value is not too close to an end point
  double len = v->getUb() - v->getLb();
  if (value < v->getLb() + minFromBds*len) {
    value = v->getLb() + minFromBds*len;
  } else if (value > v->getUb() - minFromBds*len) {
    value = v->getUb() - minFromBds*len; 
  }

  // can't branch on something that is at its bounds.
  if (!(value > v->getLb()+1e-8 && value < v->getUb()-1e-8)) {
    std::cerr << "Warning!  Branching on variable with bounds/value: [" << 
      v->getLb() << " , " << value << "  " << v->getUb() << " ]" << std::endl;
    //assert(value > v->getLb()+1e-8 && value < v->getUb()-1e-8);
  }

  Branches branches = (Branches) new BranchPtrVector();

  BranchPtr branch = (BranchPtr) new Branch();
  VarBoundModPtr mod = (VarBoundModPtr) new VarBoundMod(v, Upper, value);
  assert(!"add Mod correctly here.");
  branch->addPMod(mod);
  branch->setActivity((v->getUb()-value)/len);
  branches->push_back(branch);

  branch = (BranchPtr) new Branch();
  mod = (VarBoundModPtr) new VarBoundMod(v, Lower, value);
  assert(!"add Mod correctly here.");
  branch->addPMod(mod);
  branch->setActivity((value - v->getLb())/len);
  branches->push_back(branch);


  logger_->msgStream(LogDebug2) << "branching on " << v->getName();
  logger_->msgStream(LogDebug2) << " <= " << value << " or " 
    << " >= " << value << std::endl;

#if defined(DEBUG_CXUNIVARHANDLER)  
  std::cout << "branching on " << v->getName();
  std::cout << " <= " << value << " or " << " >= " << value << std::endl;
#endif
  
  return branches;

}

// Implement Handler::getBrMod().
ModificationPtr CxUnivarHandler::getBrMod(BrCandPtr cand, DoubleVector &x, 
                                          RelaxationPtr ,
                                          BranchDirection brdir) 
{

  // This method is used in Reliability branching.
  //XXX If you want it to be more accurate, you should add the new secant and linearizations
  //   into lmods.

  LinModsPtr lmods;
  lmods = (LinModsPtr) new LinMods();

  double minFromBds = 0.1;
  BrVarCandPtr vcand = boost::dynamic_pointer_cast <BrVarCand> (cand);
  VariablePtr v = vcand->getVar();

#if defined(DEBUG_CXUNIVARHANDLER)  
  std::cout << "Branching mod candidate (working problem) ID: " << v->getId() 
            << " address: " << (  v.get() ) << std::endl;
#endif  

  // x is a *relaxation* solution, while we have put the *original* (or
  // working) problem variables into the BrCandPtr, so we need to
  // update our value appropriately...
  
  double xval = x[v->getIndex()];
  double value = xval;  // Make sure branch value is not too close to an end point
  double len = v->getUb() - v->getLb();
  if (value < v->getLb() + minFromBds*len) {
    value = v->getLb() + minFromBds*len;
  } else if (value > v->getUb() - minFromBds*len) {
    value = v->getUb() - minFromBds*len; 
  }

  // can't branch on something that is at its bounds.
  if (!(value > v->getLb()+1e-8 && value < v->getUb()-1e-8)) {
    std::cerr << "Warning!  Branching on variable with bounds/value: [" << 
      v->getLb() << " , " << value << "  " << v->getUb() << " ]" << std::endl;
    //assert(value > v->getLb()+1e-8 && value < v->getUb()-1e-8);
  }

  if (brdir == DownBranch) {
    // down branch
    VarBoundModPtr mod = (VarBoundModPtr) new VarBoundMod(v, Upper, value);
    lmods->insert(mod);
  } else if (brdir ==  UpBranch) {
  // up branch
    VarBoundModPtr mod    = (VarBoundModPtr) new VarBoundMod(v, Lower, value);
    lmods->insert(mod);
  }

  return lmods;

}

bool IpoptEngine::presolve_()
{
  bool should_stop = false;
  VariablePtr v;
  double diff;
  bool all_fixed = true;
  int e=0;

  status_ = EngineUnknownStatus;
  // visit each variable and check bounds. Stop if bad bounds are found or if
  // all variables are fixed.
  for (VariableConstIterator it=problem_->varsBegin(); it!=problem_->varsEnd(); 
      ++it) {
    v = *it;
    diff = v->getUb() - v->getLb();
    if (diff < -etol_) {
      status_ = ProvenLocalInfeasible;
      return true;
    }
    if (fabs(diff)>etol_) {
      all_fixed = false;
    }
  }

  if (all_fixed == true) {
    double obj_value;
    double *x;
    double vio = -INFINITY;
    double act;

    x = new double[problem_->getNumVars()];
    should_stop = true;
    for (VariableConstIterator it=problem_->varsBegin();
        it!=problem_->varsEnd(); ++it) {
      v = *it;
      x[v->getIndex()] = v->getUb();
    }
    obj_value = problem_->getObjValue(x, &e);

    // check if constraints are violated.
    for (ConstraintConstIterator it=problem_->consBegin(); 
        it!=problem_->consEnd(); ++it) {
      act = (*it)->getActivity(x, &e);
      vio = std::max(act-(*it)->getUb(), (*it)->getLb()-act);
      vio = std::max(vio, 0.);
      if (vio>etol_ || e!=0) {
        break;
      }
    }

   // sol_ must be set even if the problem is infeasible because QG handler
   // may need this point for cut generation etc.
    sol_->setPrimal(x); 
    if (vio>etol_) {
      status_ = ProvenLocalInfeasible;
    } else {
      sol_->setObjValue(obj_value);
      mynlp_->setSolution(sol_);
      status_ = ProvenLocalOptimal;
    }
    delete [] x;
  }
  return should_stop;
}

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

bool OBBT::tightenVariableBound(ConstraintPtr cs, VariablePtr variable)
{
    assert(cs->hasVariable(variable));

    auto vars = cs->getVariables();

    // Store and set objective costs to zero
    std::vector<double> costs;

    for (auto var : vars)
    {
        costs.push_back(var->getCost());
    }

    // Set costs for lower bound problem
    for (auto var : vars)
    {
        if (var == variable)
            var->setCost(1); // Min. variable
        else
            var->setCost(0);
    }

    //SolverIpopt solver_min(cs);
    SolverGurobi solver_min(cs);
    SolverResult result_min = solver_min.optimize();

    // Set costs for upper bound problem
    for (auto var : vars)
    {
        if (var == variable)
            var->setCost(-1); // Max. variable
        else
            var->setCost(0);
    }

    //SolverIpopt solver_max(cs);
    SolverGurobi solver_max(cs);
    SolverResult result_max = solver_max.optimize();

    // Reset costs
    int counter = 0;
    for (auto var : vars)
    {
        var->setCost(costs.at(counter));
        counter++;
    }

    // Check for infeasibility
    if (result_min.status == SolverStatus::INFEASIBLE
        || result_max.status == SolverStatus::INFEASIBLE)
        return false;

    // Update lower bound
    if (result_min.status == SolverStatus::OPTIMAL)
    {
        if (!variable->updateLowerBound(result_min.objectiveValue))
        {
            // This should never happen!
            cout << "Min bound" << endl;
            cout << *variable << endl;
            cout << result_min << endl;
            return false;
        }
    }

    // Update upper bound
    if (result_max.status == SolverStatus::OPTIMAL)
    {
        if (!variable->updateUpperBound(-result_max.objectiveValue))
        {
            cout << std::setprecision(10) << -result_max.objectiveValue << endl;
            cout << std::setprecision(10) << variable->getLowerBound() << endl;
            cout << std::setprecision(10) << variable->getLowerBound() + result_max.objectiveValue << endl;
            // This should never happen!
            cout << "Max bound" << endl;
            cout << *variable << endl;
            cout << result_max << endl;
            return false;
        }
    }

    // No update
    return true;
}

ExpressionPtr KernelReplacer::createKernelCallLambda(const ExpressionAddress localKernel,
		const ExpressionPtr work_dim, const ExpressionPtr local_work_size, const ExpressionPtr global_work_size) {
	NodeManager& mgr = prog->getNodeManager();
	IRBuilder builder(mgr);

	ExpressionPtr k = utils::getRootVariable(localKernel).as<ExpressionPtr>();
/*
std::cout << "searching " << *k << " from " << *localKernel << std::endl;
if(kernelFunctions.empty()) std::cout << "\tnothing\n";

for_each(kernelFunctions, [](std::pair<core::ExpressionPtr, core::LambdaExprPtr> kernel) {
	std::cout << "in " << *kernel.first << std::endl;
});
*/
	// try to find coresponding kernel function
	assert(kernelFunctions.find(k) != kernelFunctions.end() && "Cannot find OpenCL Kernel");
	const ExpressionPtr local = anythingToVec3(work_dim, local_work_size);
	const ExpressionPtr global = anythingToVec3(work_dim, global_work_size);

	LambdaExprPtr lambda = kernelFunctions[k].as<LambdaExprPtr>();

//dumpPretty(lambda);
//for_each(kernelTypes[k], [&](TypePtr ty) {
//	std::cout << "->\t" << *ty << std::endl;
//});

	/*    assert(kernelArgs.find(k) != kernelArgs.end() && "No arguments for call to kernel function found");
	const VariablePtr& args = kernelArgs[k];
	const TupleTypePtr& argTypes = dynamic_pointer_cast<const TupleType>(args->getType());*/
	const VariableList& interface = lambda->getParameterList()->getElements();

	vector<ExpressionPtr> innerArgs;
	const core::lang::BasicGenerator& gen = builder.getLangBasic();

	// construct call to kernel function
//		if(localMemDecls.find(k) == localMemDecls.end() || localMemDecls[k].size() == 0) {
//std::cout << "lmd " << localMemDecls[k] << std::endl;
	TypeList kernelType = kernelTypes[k];

	/* body of a newly created function which replaces clNDRangeKernel. It contains
	 *  the kernel call
	 *  return 0;
	 */
	StatementList body;

	// Kernel variable to be used inside the newly created function
	VariablePtr innerKernel = builder.variable(builder.tupleType(kernelType));
	for(size_t i = 0; i < interface.size() -2 /*argTypes->getElementTypes().size()*/; ++i) {
//??			TypePtr argTy = utils::vectorArrayTypeToScalarArrayType(interface.at(i)->getType(), builder);
		TypePtr argTy = interface.at(i)->getType();
		TypePtr memberTy = kernelType.at(i);
		ExpressionPtr tupleMemberAccess = builder.callExpr(memberTy, gen.getTupleMemberAccess(), utils::removeDoubleRef(innerKernel),
				builder.literal(gen.getUInt8(), toString(i)), builder.getTypeLiteral(memberTy));
		ExpressionPtr argument = handleArgument(argTy, memberTy, tupleMemberAccess, body);

		innerArgs.push_back(argument);
	}

	const TypePtr vecTy = builder.vectorType(gen.getUInt8(), builder.concreteIntTypeParam(static_cast<size_t>(3)));

	// local and global size to be used inside the newly created function
	VariablePtr innerGlobal = builder.variable(vecTy);
	VariablePtr innerLocal = builder.variable(vecTy);

	// add global and local size to arguments
	innerArgs.push_back(innerGlobal);
	innerArgs.push_back(innerLocal);

	ExpressionPtr kernelCall = builder.callExpr(gen.getInt4(), lambda, innerArgs);
	body.push_back(kernelCall);							   // calling the kernel function
	body.push_back(builder.returnStmt(builder.intLit(0))); // return CL_SUCCESS

	// create function type for inner function: kernel tuple, global size, local size
	TypeList innerFctInterface;
	innerFctInterface.push_back(innerKernel->getType());
	innerFctInterface.push_back(vecTy);
	innerFctInterface.push_back(vecTy);

	FunctionTypePtr innerFctTy = builder.functionType(innerFctInterface, gen.getInt4());

	// collect inner function parameters
	VariableList innerFctParams;
	innerFctParams.push_back(innerKernel);
	innerFctParams.push_back(innerGlobal);
	innerFctParams.push_back(innerLocal);

	// create lambda for inner function
	LambdaExprPtr innerLambda = builder.lambdaExpr(innerFctTy, builder.parameters(innerFctParams), builder.compoundStmt(body));

	return builder.callExpr(gen.getInt4(), innerLambda, builder.deref(localKernel), global, local);
}

	void ArrayType::setSize(const VariablePtr& size) {
		if(size) assert_pred1(size->getNodeManager().getLangBasic().isUIntInf, size->getType());
		this->size = size;
	}

Bool logobj(LinearFunctionPtr lf, CGraphPtr nlf)
{
  UInt numvars = lf->getNumTerms();
  // VariableSet vars;
  // lf->getVars(&vars);
  VariableGroup::const_iterator it;
  VariableGroup::const_iterator begin = lf->termsBegin();
  VariableGroup::const_iterator end   = lf->termsEnd();
  // Variable nodes.
  std::map<UInt,CNode*> varnodes;
  std::map<CNode*, CNode*> lognodes;
  std::map<UInt, CNode*> logidnodes;
  std::map<CNode*, CNode*> numnodes;
  std::map<UInt, CNode*> numidnodes;
  std::map<CNode*, CNode*> mulnodes;
  std::map<UInt, CNode*> mulidnodes;
  std::set<CNode*> multnodes;
  std::map<UInt, CNode*> multidnodes;
  CNode** multnodesvec = new CNode*[numvars];
  // Iterate through variables.
  CNode* tempnode;
  CNode* lognode;
  CNode* numnode;
  CNode* mulnode;
  CNode* onenode;
  CNode* sumnode;
  // Create numnode for one.
  UInt index = 0;
  for(it=begin; it!=end; ++it) {
    VariablePtr var = it->first;
    Double varlb = var->getLb();
    if (varlb >= 0) {
      UInt varid = var->getId();
      Double coeff = it->second;
      // Create variable node and add it to a set.
      tempnode = nlf->newNode(var);
      varnodes.insert(std::pair<UInt,CNode*>(varid,tempnode));
      onenode = nlf->newNode(1.0);
      // Create sumnode.      
      sumnode = nlf->newNode(OpPlus, onenode, tempnode);
      // Create log nodes.
      lognode = nlf->newNode(OpLog, &sumnode, 1);
      lognodes.insert(std::pair<CNode* , CNode*>(tempnode, lognode));
      logidnodes.insert(std::pair<UInt, CNode*>(varid, lognode));
      // Create num nodes.
      numnode = nlf->newNode(coeff);
      numnodes.insert(std::pair<CNode*, CNode*>(tempnode, numnode));
      numidnodes.insert(std::pair<UInt, CNode*>(varid, numnode));
      // Create mul nodes.
      mulnode = nlf->newNode(OpMult, numnode, lognode);
      multnodes.insert(mulnode);
      multidnodes.insert(std::pair<UInt, CNode*>(varid,mulnode));
      multnodesvec[index] = mulnode;
      index++;
    }
  }
  
  // This is the root node.
  CNode* sumlistnode = nlf->newNode(OpSumList, multnodesvec, index);  
  nlf->setOut(sumlistnode);
  nlf->finalize();
  
  return true;
}

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

void POP03::runProblem()
{
    // Problem 3 (Nataraj)
    // cout << "\n\nSolving problem P03..." << endl;

    int dim = 4;

    // x1,x2,l1,l2
    std::vector<double> costs = {1, 0, 0, 0};
    std::vector<double> lb = {-10,-10,0,-1000};
    std::vector<double> ub = {10,10,100,1000};
    std::vector<double> z0(dim,0);

    std::vector<VariablePtr> vars;
    for (int i = 0; i < dim; i++)
    {
        auto var = std::make_shared<Variable>(costs.at(i), lb.at(i), ub.at(i));
        vars.push_back(var);
    }

    ConstraintSetPtr cs = std::make_shared<ConstraintSet>();

    { // x1^2 = l1
        VariablePtr var = vars.at(0);
        std::vector<double> thislb = {var->getLowerBound()};
        std::vector<double> thisub = {var->getUpperBound()};

        std::vector<unsigned int> deg = {2};

        DenseVector c(3);
        c.setZero();
        c(0) = 0;
        c(1) = 0;
        c(2) = 1;

        DenseMatrix T = getTransformationMatrix(deg, thislb, thisub);
        DenseMatrix coeffs = T*c;

        std::vector< std::vector<double> > knots = getRegularKnotVectors(deg, thislb, thisub);

        BSpline bs(coeffs.transpose(), knots, deg);

        std::vector<VariablePtr> cvars = {var, vars.at(2)};
        ConstraintPtr cbs = std::make_shared<ConstraintBSpline>(cvars, bs, true);
        cs->add(cbs);
    }

    { // x1^3 = l1
        VariablePtr var = vars.at(0);
        std::vector<double> thislb = {var->getLowerBound()};
        std::vector<double> thisub = {var->getUpperBound()};

        std::vector<unsigned int> deg = {3};

        DenseVector c(4); c.setZero();
        c(0) = 0;
        c(1) = 0;
        c(2) = 0;
        c(3) = 1;

        DenseMatrix T = getTransformationMatrix(deg, thislb, thisub);
        DenseMatrix coeffs = T*c;

        std::vector< std::vector<double> > knots = getRegularKnotVectors(deg, thislb, thisub);

        BSpline bs(coeffs.transpose(), knots, deg);

        std::vector<VariablePtr> cvars = {var, vars.at(3)};
        ConstraintPtr cbs = std::make_shared<ConstraintBSpline>(cvars, bs, true);
        cs->add(cbs);
    }

    { // x2 - l1 <= 0, x2 - l2 <= 0

        DenseMatrix A(2,3);
        A.setZero();
        A(0,0) = -1; A(0,1) = 1;
        A(1,0) = 1; A(1,1) = 2; A(1,2) = -1;

        DenseVector b; b.setZero(2); b(1) = -1e-5;

        std::vector<VariablePtr> cvars = {
            vars.at(1),
            vars.at(2),
            vars.at(3)
        };
        ConstraintPtr lincon = std::make_shared<ConstraintLinear>(cvars, A, b, false);
        cs->add(lincon);
    }

    BB::BranchAndBound bnb(cs);
    Timer timer;
    timer.start();
    SolverResult res = bnb.optimize();
    timer.stop();
    cout << "Time: " << timer.getMilliSeconds() << " (ms)" << endl;

    if (res.status == SolverStatus::OPTIMAL)
        fopt_found = res.objectiveValue;
}