C++ MIP类代码示例

void LinearConstraint::augmentStructure(MIP &mip, Cost &cost, int varindex, map<Value, Cost> &delta) {

	for (map<Value, Cost>::iterator i = delta.begin(); i != delta.end();i++) {
		
		int var1 = mapvar[varindex][i->first];
		mip.objCoeff(var1, mip.objCoeff(var1)-i->second); // update unary cost
		if (mip.sol(var1) == 1){ // using this value?
			cost -= i->second;
		}
	}

}

void LinearConstraint::checkRemoved(MIP &mip, Cost &cost, vector<int> &rmv) {


	pair<Cost, bool> result;
	vector<int> cDomain, cDomain2;
	//bool deleted = false;
	bool flag = false;
	for (int i=0;i<arity_;i++) {
		cDomain.clear();
		getDomainFromMIP(mip, i, cDomain); // getDomain
		sort(cDomain.begin(), cDomain.end());
		
		
		EnumeratedVariable* y = (EnumeratedVariable*)getVar(i);
		for (EnumeratedVariable::iterator v = y->begin(); v != y->end();++v) {
		
			vector<int>::iterator it = find(cDomain.begin(), cDomain.end(), *v);
			if (it == cDomain.end()) {
				cout << "non exist a value ?" << endl;
				for (vector<int>::iterator v=cDomain.begin();v != cDomain.end();v++) {
					cout << *v << " ";
				} cout << endl;
				for (EnumeratedVariable::iterator v = y->begin(); v != y->end();++v) {
					cout << *v << " ";
				} cout << endl;
				exit(0);
			}
			cDomain.erase(it);
			//deleted = true;
		}
		
		
		if (!cDomain.empty()) {
			cDomain2.clear();
			rmv.push_back(i);
			for (vector<int>::iterator v=cDomain.begin();v != cDomain.end();v++) {
				int var1 = mapvar[i][*v];
				if (mip.sol(var1) == 1){ // checking if this value is being used
					flag = true;
				}
				mip.colUpperBound(var1, 0); // removeDomain
			}
			//deleted = true;
		}
	}
	if (flag){
		cost = solveMIP(); // solve
	}
}

  static void GenerateMatrix (const FEL & bfel, const MIP & mip,
			      MAT & mat, LocalHeap & lh)
  {
    // must get the right elements, otherwise an exception is thrown.

    const CompoundFiniteElement & cfel = 
      dynamic_cast<const CompoundFiniteElement&> (bfel);

    // a scalar H1 element
    const ScalarFiniteElement<2> & fel_u = 
      dynamic_cast<const ScalarFiniteElement<2>&> (cfel[0]);
    const ScalarFiniteElement<2> & fel_p = 
      dynamic_cast<const ScalarFiniteElement<2>&> (cfel[2]);
    
    int nd_u = fel_u.GetNDof();
    int nd_p = fel_p.GetNDof();
    
    // transformation of derivatives from reference element to general element:
    FlatMatrixFixWidth<2> gradu(nd_u, lh);
    fel_u.CalcMappedDShape (mip, gradu);

    // the shape functions of the pressure
    FlatVector<> vecp(nd_p, lh);
    fel_p.CalcShape (mip.IP(), vecp);

    mat = 0;

    // the first nd_u shape functions belong to u_x, the next nd_u belong to u_y:
    mat.Rows(0,2).Cols(cfel.GetRange(0)) = Trans (gradu);
    mat.Rows(2,4).Cols(cfel.GetRange(1)) = Trans (gradu);

    // ... and finally nd_p shape functions for the pressure:
    mat.Row(4).Range(cfel.GetRange(2)) = vecp;
  }

    static void GenerateMatrix (const AFEL & bfel, const MIP & sip,
                                MAT & mat, LocalHeap & lh)
    {
      HeapReset hr(lh);
      const HDivDivFiniteElement<2> & fel = 
        dynamic_cast<const HDivDivFiniteElement<2>&> (bfel);
      
      int nd = fel.GetNDof();
      
      Mat<3,2> jac = sip.GetJacobian();
      double det = fabs (sip.GetJacobiDet());

      FlatMatrix<> shape(nd, 3, lh);
      fel.CalcShape (sip.IP(), shape);
      
      Mat<3,9> trans;
      for (int i = 0; i < 3; i++)
        {
          Mat<2> sigma_ref;
          sigma_ref = 0.0;
          switch (i)
            {
            case 0: sigma_ref(0,0) = 1.0; break;
            case 1: sigma_ref(1,1) = 1.0; break;
            case 2: sigma_ref(0,1) = sigma_ref(1,0) = 1.0; break;
            }
          auto hm = jac * sigma_ref;
          auto sigma = hm * Trans(jac);
          sigma *= (1.0 / sqr(det));
          
          trans ( i, 0 ) = sigma(0,0);
          trans ( i, 1 ) = sigma(0,1);
          trans ( i, 2 ) = sigma(0,2);
          trans ( i, 3 ) = sigma(1,0);
          trans ( i, 4 ) = sigma(1,1);
          trans ( i, 5 ) = sigma(1,2);
          trans ( i, 6 ) = sigma(2,0);
          trans ( i, 7 ) = sigma(2,1);
          trans ( i, 8 ) = sigma(2,2);
        }
      mat = Trans(trans) * Trans (shape);
    }

void LinearConstraint::findProjection(MIP &mip, Cost &cost, int varindex, map<Value, Cost> &delta) {

	pair<Cost, bool> result;
	delta.clear();
	EnumeratedVariable* x = (EnumeratedVariable*)getVar(varindex);
	for (EnumeratedVariable::iterator j = x->begin(); j != x->end(); ++j) {
		
		int var1 = mapvar[varindex][*j];
		int tmp = cost;
		cost = tmp = mip.augment(var1); // make sure this value is used...
		
		assert(tmp >= 0);
		delta[*j] = tmp;
	}

}

	virtual int solve() { //return the optimal value from the MIP
             if (solver) solver->solve();     
             return 0;
        }

	virtual int sol(int var1) {
            if (solver) return solver->sol(var1);     
            return 0;
        }

	virtual int solValue() {
            if (solver) return solver->solValue();     
            return 0;
        }

	virtual int restore() {
            if (solver) return solver->restore();                 
            return 0;            
        }

	virtual void objCoeff(int var1, int i) { // Set the coefficient of the variable 
             if (solver) solver->objCoeff(var1, i);     
        }

	virtual int augment(int var1) {
            if (solver) return solver->augment(var1);     
            return 0;
        }

	virtual void addBool(int n) {         
                if (solver) solver->addBool(n);         
        }

	virtual void addInt(int n) {
                if (solver) solver->addInt(n); 
        }

	virtual void addRows(int n) {
            if (solver) solver->addRows(n); 
        }

	virtual void end() {
            if (solver) solver->end(); 
        }