smknap.cc

#include <iostream>
#include <cmath>
#include<scil/scil.h>
#include<scil/constraints/submodular.h>

using namespace std;
using namespace SCIL;

class mysub : public submodular {
   public:
      std::map<var, double>& co;

      mysub(var& v_, var_map<int>& VM_, std::map<var, double>& co_, int n_) 
         : submodular(v_, VM_, n_), co(co_)
      {
      }

   private:
      double value(solution& S) {
         double w = 0;
         for( int i = 0; i < n; i++ )
            w += S.value(VM[i]) * co[VM[i]];
         return sqrt(w);
      }
};


int main()
{
   int n;
   double b, k;
   double lb = 0.;
   double ub = 0.;
   cout << "Type n: ";
   cin >> n;
   double* s = new double[n];
   double* t = new double[n];
   double* bs = new double[n];
   //compute initial lower and upper bound on k
   for( int i = 0; i < n; i++ ) {
      cin >> s[i];
      if( s[i] < 0. )
         lb += s[i];
      else
         ub += s[i];
   }
   for( int i = 0; i < n; i++ )
      cin >> t[i];
   cin >> b;
   
   int no = 0;
   int ns = 0;
   int nl = 0;
   while( ub-lb > 1 ) {
      no++;
      ILP_Problem IP(Optsense_Min);
      IP.primal_bound(b + .1);
      IP.set_silent();
      var_map<int> VM;
      std::map<var, double> co;
      row r;
      var v = IP.add_variable(1., -1e4, 1e25, Vartype_Float);
      for( int i = 0; i < n; i++ ) {
         VM[i] = IP.add_binary_variable(0.);
         r += VM[i] * s[i];
      }
      for( int i = 0; i < n; i++ ) {
         co[VM[i]] = t[i];
      }
      IP.add_sym_constraint(new mysub(v, VM, co, n));
      k = lb + (ub - lb)/2.;
      IP.add_basic_constraint(r>=k, Static);

      IP.optimize();
      solution sol = IP.get_solution();
      cout << "\n#" << no << "  lb: " << lb << "  ub: " << ub << "  k: " << k << endl;
      cout << "feasible: " << IP.feasibleFound() << "  status: " << IP.status();
      cout << "  opt value: " << IP.get_optimum() << endl;
      for( int i = 0; i < n; i++ )
         cout << sol.value(VM[i]) << "  ";
      cout << endl;
      ns += IP.nSub();
      nl += IP.nLp();

      if( IP.feasibleFound() && IP.status() == 0 && IP.get_optimum() <= b ) {
         lb = k;
         for( int i = 0; i < n; i++ )
            bs[i] = sol.value(VM[i]);
      }
      else
         ub = k;

   }

   double val = 0.;
   cout << endl << "Best solution: ";
   for( int i = 0; i < n; i++ ) {
      cout << bs[i] << " ";
      val += bs[i] * s[i];
   }
   cout << endl << "value: " << val << endl;
   cout << "#solved IPs: " << no << "   #subs: " << ns << "   #LPs: " << nl << endl;

   delete[]s;
   delete[]t;
   delete[]bs;
   return 0;
}

---