disjunctive.cc

#include<scil/constraints/disjunctive.h>
#include<scil/cons_map.h>

#define infinity 1e+30

class disjunctive::disjunctive_cut : public cons_obj {
 private:
  cons_obj* c;
  h_array<var,var>& VM;
  var v;
  bool de;

 public:
  disjunctive_cut(cons_obj* c_, h_array<var,var>& VM_, var v_, bool de_=false) 
    : cons_obj(c_->sense(), 0), VM(VM_), de(de_) {
    c=c_;
    v=v_;
  };

  void non_zero_entries(row& r) {
    row r1;
    c->non_zero_entries(r1);

    row_entry re;
    forall(re, r1) {
      if(VM[re.Var]!=nil)
        r+=re.coeff*VM[re.Var];
    };
    r-=c->rhs()*v;
    //cout<<r1.size()<<" size of row\n";
  };

  ~disjunctive_cut() {
    if(de) delete c;
  };
};

#ifdef ABACUS_COMPILER_GCC
class my_subproblem : public subproblem {
#else
class my_subproblem : public Std_subproblem_ {
#endif
private:
  h_array<var,var>& V;
  var z;
  subproblem& S;
  subproblem& S1;
  double d;

public:
  my_subproblem(h_array<var,var>& V_, var z_, subproblem& S_,
                subproblem& S1_) :
#ifdef ABACUS_COMPILER_GCC
    subproblem(&S_.GM), 
#endif
    V(V_), z(z_), S(S_), S1(S1_) {
    d=S1.value(z);
  };

  double
  value(var v) {
    if(V[v]==nil) return 0;
    return S1.value(V[v])/d;
  };

  cons add_basic_constraint(cons_obj* c, Activation a = Static) {
    //cout<<"add modified constraint\n";
    return S1.add_basic_constraint(new disjunctive::disjunctive_cut(c, V, z, true));
  };

  bool my_separate() {
    sym_constraint* sc;
    bool b=false;
    forall_sym_constraints(sc, S) {
      disjunctive* d=dynamic_cast<disjunctive*>(sc);
      if(d==nil) {
        sym_constraint::status s=sc->separate(*this);
        if(s==sym_constraint::constraint_found) b=true;
      };
    };
    return b;
  };

  virtual
  ~my_subproblem() {
  };

  void before_delete() {
#ifndef ABACUS_COMPILER_GCC
    fatherp=this;
#endif
    //sym_constraints.clear();
  };
};

class separate_disj : public sym_constraint {
private:
  h_array<var,var>& v1;
  h_array<var,var>& v2;
  var z1, z2;
  bool b1;
  //cons_map<var>& VM1;
  //cons_map<var>& VM2;
  cons b;
  list<cons>& A;

  h_array<var, double> D1;
  h_array<var, double> D2;
  cons_map<cons>& C1;
  cons_map<cons>& C2;

  subproblem& OS;
  h_array<var,bool>& IN;
  h_array<cons, double> Vv;
  int count;

public:

  separate_disj(h_array<var,var>& v1_, h_array<var,var>& v2_, 
                var z1_, var z2_,
                //cons_map<var>& VM1_, cons_map<var>& VM2_, 
                cons b_, 
                list<cons>& A_, cons_map<cons>& C1_, cons_map<cons>& C2_,
                subproblem& OS_, h_array<var, bool>& IN_)
    : v1(v1_), v2(v2_), z1(z1_), z2(z2_), // VM1(VM1_), VM2(VM2_), 
      b(b_), A(A_),
      C1(C1_), C2(C2_), OS(OS_), IN(IN_) {
    count=0;
  };

  status feasible(solution& S) {
    if(count<1000) return infeasible_solution;
    return feasible_solution;
  };

  status separate(subproblem& s) {
    count++;
    cout<<count<<" counter\n";
    if(count<=1) {
      solution sol;
      my_subproblem* SP=new my_subproblem(v1, z1, OS, s);
      b1=SP->my_separate();
      SP->before_delete();
      delete SP;
      SP=new my_subproblem(v2, z2, OS, s);
      bool b2=SP->my_separate();
      SP->before_delete();
      delete SP;
      b1=b1 || b2;
      
      if(b1)
        return constraint_found;
    } else b1=false;

    cout<<"SAVE SOLUTOIN\n";
    //V1v.clear();
    //V2v.clear();
    cons c;
    forall_basic_constraints(c, s) { // TODO improve: only some values need to be saved
      Vv[c]=s.value(c);
    };

    var v;
    forall_variables(v, OS) if(IN[v]) {
      D1[v]=s.red_cost(v1[v]);
      D2[v]=s.red_cost(v2[v]);
    };

    double d;
    forall(c, A) {
      d=s.value(C1[c]);
      if(d!=0) {
        row r1;
        c.cons_pointer()->non_zero_entries(r1);
        row_entry re;
        forall(re, r1) {
          D1[re.Var]+=re.coeff*d;
        };
      };
      d=s.value(C2[c]);
      if(d!=0) {
        row r1;
        c.cons_pointer()->non_zero_entries(r1);
        row_entry re;
        forall(re, r1) {
          D2[re.Var]+=re.coeff*d;
        };
      };
    };

    //forall_defined(c, VM2) {
    //  Vv[c]=s.value(c);
    //};
    //Vv[b]=s.value(b);

    count=1000;
    return fathom;
  };

  double value1(var v) {
    return D1[v];
  };

  double value2(var v) {
    return D2[v];
  };

  double value(cons c) {
    return Vv[c];
  };
};
    
disjunctive::status disjunctive::separate(subproblem& s) {

  // Compute a cut?

  if (s.found_constraint()) {
    return no_constraint_found;
  };

  cout<<"DISJUNCTIVE SEPARATION\n";
  
  status sta=no_constraint_found;
  list<cons_obj* > L;

  // Compute fractional variables

  list<var> FRAC;
  var v, w;
  double d;
  h_array<var,bool> IN;
  forall_variables(v, s) {
    d=s.value(v);
    if ((d>0.001) && (d<0.999)) 
      FRAC.append(v);
    if(d>0.001) IN[v]=true; else IN[v]=false;
  }

  int nk=0;     
  forall(v, FRAC) if((s.value(v)>0.2) && (s.value(v)<0.8) && (nk<10)) {
    nk++;
    cout<<nk<<"th round\n";
    // Create Initial Cut generating LP

    ILP_Problem P(Optsense_Min);
#ifdef ABACUS_COMPILER_GCC
    P.set_silent();
#endif
    P.add_variable(0,0,1,Vartype_Integer);
    h_array<var,var> v1;
    h_array<var,var> v2;
    var z1, z2, wo;
    list<cons> A;
    cons_map<cons> C1;
    cons_map<cons> C2;

    // Create Variables;
    z1=P.add_variable(0,0,1,Vartype_Float);
    z2=P.add_variable(0,0,1,Vartype_Float);

    forall_variables(w, s) {
      d=s.value(w);
      if (IN[w]) {
         v1[w]=P.add_variable(0,0,1,Vartype_Float);
         v2[w]=P.add_variable(0,0,1,Vartype_Float);
      } else {
        if (d<=0.01) {
          v1[w]=nil;
          v2[w]=nil;
        } else { // DOES NOT HAPPEN
          cout<<"HEIRERERTAERGARTHJFTAJSDFG\n";
          v1[w]=z1; //row(1);
          v2[w]=z2; //row(1);
        };
      };
    };
    wo=P.add_variable(1,0,100,Vartype_Float);

    // Create generic constraints
    cons b=P.add_basic_constraint(z1+z2==1);
    cons_map<var> VM1;
    cons_map<var> VM2;

    //forall(w, FRAC) {
    forall_variables(w, s) if(IN[w]) {
      VM1[w]=P.add_basic_constraint(v1[w]+v2[w]+wo>=s.value(w));
      VM2[w]=P.add_basic_constraint(v1[w]+v2[w]-wo<=s.value(w));
    };

    // Create constraints for constraints
    cons c;
    forall_basic_constraints(c, s) 
      //if((fabs(c.cons_pointer()->slack(s))<0.01) || 
      //         (c.cons_pointer()->Act==Static)) {
      if(true) {
      A.append(c);
      C1[c]=P.add_basic_constraint(new disjunctive_cut(c.cons_pointer(), v1, z1));
      C2[c]=P.add_basic_constraint(new disjunctive_cut(c.cons_pointer(), v2, z2));
    };

    cons u1=P.add_basic_constraint((row) v1[v]<=0);
    cons u2=P.add_basic_constraint((row) v2[v]>=z2);
    //forall(w, FRAC) {
    forall_variables(w, s) if(IN[w]) {
      P.add_basic_constraint((row) v2[w]<=z2);
      P.add_basic_constraint((row) v1[w]<=z1);
    };

    // Separate to strengthen
    separate_disj* sd=new separate_disj(v1,v2, z1, z2, //VM1, VM1, 
                                        b, A, C1, C2, s, IN);
    P.add_sym_constraint(sd);

    // Compute cut
    P.optimize();

    // Retrieve Cut
    row r; double d, d1, vi=0;
    //forall(w, FRAC) {
    forall_variables(w, s) {
      if(IN[w]) {
        d1=-sd->value(VM1[w])-sd->value(VM2[w]);
        /*
        d1=leda_min(sd->value1(w),sd->value2(w));
        if(w.var_pointer()==v.var_pointer()) 
          d1=leda_min(sd->value1(w)+sd->value(u1), sd->value2(w)+sd->value(u2));
        */
      } else {
        d1=leda_max(sd->value1(w),sd->value2(w));
      }
      /*
      if(fabs(d-d1)<0.001) { if(fabs(d)>0.001) cout<<"j"; } else {
        cout<<endl<<d<<" "<<sd->value(VM1[w])
            <<" "<<sd->value(VM2[w])<<"n"
            <<(w.var_pointer()==v.var_pointer())<<"\n";
        cout<<d1<<" "<<sd->value1(w)<<" "<<sd->value2(w)<<endl;
      }
      */
      if(d1!=0) {
        r+=d1*w;
        vi+=d1*s.value(w);
      }
    };
    if(true) { // TODO better criterion
      L.append(r >= sd->value(b));
      cout<<vi-sd->value(b)<<"Violation \n";
      if(vi-sd->value(b)>0) for(;;);
      sta=constraint_found;
    };

    delete sd;
  };

  cons_obj* c;
  forall(c, L) s.add_basic_constraint(c);

  return sta;
};

---