QuadFit.cc

#include<scil/scil.h>
#include<scil/core/monomial.h>
#include<list>

#include<iostream>
#include<stdio.h>
#include<vector>

#define SHIFT 0.5
// when deleting a monomial X*Y with coefficient CF, 
// increase coefficients of X and Y by SHIFT*CF each

#define SMALL(A,B) ((A) < (B) && (A) > -(B))
// is |A| < B ?

using namespace SCIL;

int main(int argc, char **argv) {

  int num_vars;       // number of variables (positions)
  int num_edges;      // number of edges (pairs with interaction)
  int num_sams;       // number of scenarios

  int *head,*tail;    // end nodes of edges
  double **coeff;     // coefficients
  double *m;          // target values

  double cuteps = 0;  // delete monomials with smaller coefficients
  double cst;         // constant part of the objective function

  // read data

  if(argc > 1) cuteps = atof(argv[1]);

  cin>>num_vars>>num_edges>>num_sams;

  head = new int[num_edges];
  tail = new int[num_edges];
  coeff = new double*[num_edges];
  for(int e = 0; e < num_edges; e++) coeff[e] = new double[num_sams];
  for(int e = 0; e < num_edges; e++) {
    cin>>head[e]>>tail[e];
    head[e]--;
    tail[e]--;
    for(int s = 0; s < num_sams; s++) cin>>coeff[e][s];
  }

  cst = 0;
  m = new double[num_sams];
  for(int s = 0; s < num_sams; s++) {
    cin>>m[s];
    cst += m[s]*m[s];
  }

  // determine coefficients for linear and quadratic monomials

  double *lobj = new double[num_vars];
  double *qobj = new double[num_edges];

  double cf;
  for(int v = 0; v < num_vars; v++) lobj[v] = 0;
  for(int e = 0; e < num_edges; e++) {
    cf = 0;
    for(int s = 0; s < num_sams; s++) {
      cf -= 2.0*m[s]*coeff[e][s];
      cf += coeff[e][s]*coeff[e][s];
    }
    qobj[e] = cf;
  }
  for(int e1 = 0; e1 < num_edges-1; e1++) {
    for(int e2 = e1+1; e2 < num_edges; e2++) {
      cf = 0;
      for(int s = 0; s < num_sams; s++) cf += 2.0*coeff[e1][s]*coeff[e2][s];
      if(SMALL(cf,cuteps)) {
        qobj[e1] += SHIFT*cf;
        qobj[e2] += SHIFT*cf;
      }
    }
  }
  for(int e = 0; e < num_edges; e++) {
    if(SMALL(qobj[e],cuteps)) {
      lobj[head[e]] += SHIFT*qobj[e];
      lobj[tail[e]] += SHIFT*qobj[e];
    }
  }

  // set up polynomial problem
  
  ILP_Problem IP(Optsense_Min);

  var *lvar = new var[num_vars];
  //var *qvar = new var[num_edges];
  std::vector<var> mon(4);

  for(int v = 0; v < num_vars; v++) 
    lvar[v] = IP.add_binary_variable(lobj[v]);
  for(int e = 0; e < num_edges; e++) {
    if(!SMALL(qobj[e],cuteps))
      /*qvar[e] = */ IP.add_polynomial(lvar[head[e]]*lvar[tail[e]]*qobj[e]);
  }
  monomial *nm;
  std::list<monomial*> lm;
  for(int e1 = 0; e1 < num_edges-1; e1++) {
    for(int e2 = e1+1; e2 < num_edges; e2++) {
      cf = 0;
      for(int s = 0; s < num_sams; s++) cf += 2.0*coeff[e1][s]*coeff[e2][s];
      if(!SMALL(cf,cuteps)) {
        mon[0] = lvar[head[e1]];
        mon[1] = lvar[tail[e1]];
        mon[2] = lvar[head[e2]];
        mon[3] = lvar[tail[e2]];
        nm = new monomial(cf,mon);
        IP.add_polynomial(*nm);
        lm.push_back(nm);
      }
    }
  }

  delete[] lobj;
  delete[] qobj;

  // optimize

  IP.optimize();

  // print solution

  bool *opt = new bool[num_vars];
  for(int v = 0; v < num_vars; v++)
    opt[v] = (IP.get_solution(lvar[v]) > 0.5);
  double robj = 0;
  for(int s = 0; s < num_sams; s++) {
    cf = m[s];
    for(int e = 0; e < num_edges; e++) 
      if(opt[head[e]] && opt[tail[e]]) cf -= coeff[e][s];
    robj += cf*cf;
  }
  cout<<"solution:";
  for(int v = 0; v < num_vars; v++) if(opt[v]) cout<<" "<<v;
  cout<<endl;
  cout<<"approximate objective value: "<<IP.get_optimum()+cst<<endl; 
  cout<<"real objective value:        "<<robj<<endl;
 
  // free
  list<monomial*>::iterator it;
  for( it = lm.begin(); it != lm.end(); it++ )
  delete (*it);

  delete[] opt;
         
  delete[] lvar;
  //delete[] qvar;

  delete[] head;
  delete[] tail;
  for(int e = 0; e < num_edges; e++) delete[] coeff[e];
  delete[] coeff;
  delete[] m;

  // return

  return 0;
}

---