tpp.c

#line 2 "tpp.lw"
#include<LEDA/graph.h>
#include<scil/scil.h>
#include<scil/constraints/tour.h>

using namespace SCIL;
using namespace LEDA;
using namespace std;

//#include"cp_cuts.h"

/*
class rrt : public cp_sym_constraint {
  private:
  var**** V;
  int n;

  public:
  rrt(int n_, var**** V_) {
    V=V_;
    n=n_;
  }

  virtual
  void infer(infer_status& is) {
    // cout<<"trying to infer\n";

    list<var> dummy; // TODO use it!
    var one_fixed_i, one_fixed_o;
    
    int nv[n];
    list<var> nvn[n];

    for(int i=0; i<n; i++) { // Test for Tours
      for(int j=0; j<n; j++) nv[j]=0;
      for(int l=0; l<2*n-3; l++) { // For every time step but last
        int lb_tci=0, ub_tci=0, lb_tco=0, ub_tco=0;
        int lb_nci, ub_nci, lb_nco, ub_nco;
        for(int j=0; j<n; j++) {
          lb_nci=0; ub_nci=0; lb_nco=0; ub_nco=0;
          
          // count node degrees
          for(int k=0; k<n; k++) {
            lb_nci+=is.lower_bound(V[i][k][j][l]);
            ub_nci+=is.upper_bound(V[i][k][j][l]);
            if(lb_nci+lb_tco>1) {
              dummy.append(V[i][k][j][l]);
              dummy.append(one_fixed_i);
              is.unsolvable(dummy);
              dummy.clear();
              return;
            };
            if(is.one_fixed(V[i][k][j][l])) {
              one_fixed_i=V[i][k][j][l];
              nv[j]++;
              nvn[j].append(V[i][k][j][l]);
            };
            lb_nco+=is.lower_bound(V[i][j][k][l+1]);
            ub_nco+=is.upper_bound(V[i][j][k][l+1]);
            if(lb_nco+lb_tco>1) {
              dummy.append(V[i][j][k][l+1]);
              dummy.append(one_fixed_o);
              is.unsolvable(dummy);
              dummy.clear();
              return;
            };
            if(is.one_fixed(V[i][j][k][l+1])) one_fixed_o=V[i][j][k][l+1];
          }

          if(lb_nci==1) { 
            // if one incoming is fixed to one: 
            // fix all other incoming in column to 0
            // fix all outgoing edges of other nodes to 0
            dummy.append(one_fixed_i);
            for(int j1=0; j1<n; j1++)
              for(int k=0; k<n; k++) {
                if(!is.fixed(V[i][k][j1][l])) {
                  is.fix(V[i][k][j1][l], 0, dummy);
                };
                if((j1!=j) && (!is.fixed(V[i][j1][k][l+1]))) {
                  is.fix(V[i][j1][k][l+1], 0, dummy);
                };
              }
         
            if((ub_nco<=1) && (lb_nco==0)) {
              var unfixed;
              for(int k=0; k<n; k++) 
                if(is.fixed(V[i][j][k][l+1])) { 
                  dummy.append(V[i][j][k][l+1]);
                } else {
                  unfixed=V[i][j][k][l+1];
                }
              if(ub_nco==0) {
                is.unsolvable(dummy);
                return;
              }
              is.fix(unfixed, 1, dummy);
              ub_nco=1; lb_nco=1;
            };
           
            dummy.clear();
          };

          if(lb_nco==1) { 
            // analog above
            dummy.append(one_fixed_o);
            for(int j1=0; j1<n; j1++) 
              for(int k=0; k<n; k++) {
                if(!is.fixed(V[i][j1][k][l+1])) {
                  is.fix(V[i][j1][k][l+1], 0, dummy);
                };
                if((j1!=j) && (!is.fixed(V[i][k][j1][l]))) {
                  is.fix(V[i][k][j1][l],0, dummy);
                }
              }

            if((ub_nci<=1) && (lb_nci==0)) {
              var unfixed;
              for(int k=0; k<n; k++) 
                if(is.fixed(V[i][k][j][l])) { 
                  dummy.append(V[i][k][j][l]);
                } else {
                  unfixed=V[i][k][j][l];
                }
              if(ub_nci==0) {
                is.unsolvable(dummy);
                return;
              }
              is.fix(unfixed, 1, dummy);
              lb_nci=1; ub_nci=1;
            };
            dummy.clear();
          };
  
          if((lb_nci==1) && (ub_nco==1) && (lb_nco==0)) { 
            // if one incoming edge and only one outgoing edge left: fix it
            var unfixed;
            dummy.append(one_fixed_i);
            for(int k=0; k<n; k++) 
              if(!is.fixed(V[i][j][k][l+1])) {
                unfixed=V[i][j][k][l+1];
              } else dummy.append(V[i][j][k][l+1]);
                                  
            is.fix(unfixed, 1, dummy);
            dummy.clear();
          };
          if((lb_nco==1) && (ub_nci==1) && (lb_nci==0)) {
            // as above
            var unfixed;
            dummy.append(one_fixed_o);
            for(int k=0; k<n; k++) 
              if(!is.fixed(V[i][k][j][l])) {
                unfixed=V[i][k][j][l];
              } else dummy.append(V[i][k][j][l]);
            is.fix(unfixed, 1, dummy);
            dummy.clear();
          };

          if((ub_nci==0) && (ub_nco!=0)) {
            // if indegree is 0 outdegree is 0
            for(int k=0; k<n; k++) {
              dummy.append(V[i][k][j][l]);
            };
            for(int k=0; k<n; k++) if(!is.fixed(V[i][j][k][l+1])) {
              is.fix(V[i][j][k][l+1], 0, dummy);
            };
            dummy.clear();
          };

          if((ub_nco==0) && (ub_nci!=0)) {
            // if outdegree is 0 indegree is 0
            for(int k=0; k<n; k++) {
              dummy.append(V[i][j][k][l+1]);
            };
            for(int k=0; k<n; k++) if(!is.fixed(V[i][k][j][l])) {
              is.fix(V[i][k][j][l], 0, dummy);
            };
            dummy.clear();
          };



          // add to total degrees
          lb_tci+=lb_nci;
          ub_tci+=ub_nci;
          lb_tco+=lb_nco;
          ub_tco+=ub_nco;
        };

        // infer on total degrees
        if(ub_tci==0) {
          is.unsolvable(dummy);
          cout<<"not implemented jet 1\n";
          return;
        };
        if (ub_tco==0) {
          is.unsolvable(dummy);
          cout<<"not implemented jet 2\n";
          return;
        }
        if((ub_tci==1) && (lb_tci==0)) {
          // only one in column left
          cout<<"one in column left\n";
          var unfixed;
          for(int j=0; j<n; j++)
            for(int k=0; k<n; k++) 
              if(is.fixed(V[i][k][j][l])) {
                dummy.append(V[i][k][j][l]);
              } else unfixed=V[i][k][j][l];
          
          is.fix(unfixed, 1, dummy);
          dummy.clear();
        };

        if((ub_tco==1) && (lb_tco==0)) {
          // only one in column left
          var unfixed;
          for(int j=0; j<n; j++)
            for(int k=0; k<n; k++) 
              if(is.fixed(V[i][j][k][l+1])) {
                dummy.append(V[i][j][k][l+1]);
              } else unfixed=V[i][j][k][l+1];
          is.fix(unfixed, 1, dummy);
        };
      }
      for(int j=0; j<n; j++) {
        if((i!=j) && (nv[j]>1)) {
          var v;
          forall(v, nvn[j]) dummy.append(v);
          is.unsolvable(dummy);
          dummy.clear();
          return;
        }
        if((i==j) && (nv[j]>n-1)) {
          var v;
          forall(v, nvn[j]) dummy.append(v);
          is.unsolvable(dummy);
          return;
        }
      }
    };
  };
};
*/

int main() {

  int n, U;
  
  ifstream is("data.txt");

  is>>n;
  is>>U;
  cout<<n<<" cities\n";
  cout<<U<<" succesive\n";

  int D[n][n];

  for(int i=0; i<n; i++) {
    for(int j=0; j<n; j++) {
      is>>D[i][j];
      cout<<D[i][j]<<" ";
    }
    cout<<endl;
  }

  ILP_Problem IP(Optsense_Min);

  //cp_cuts* CPC=new cp_cuts;
  cons c;

  var**** V;
  V=new var***[n];
  for(int i=0; i<n; i++) {
    V[i]=new var**[n];
    for(int j=0; j<n; j++) {
      V[i][j]=new var*[n];
      for(int k=0; k<n; k++) {
        V[i][j][k]=new var[2*n-2];
      }
    }
  };

  int jj, kk;
  for(int i=0; i<n; i++) {
    for(int j=0; j<n; j++) {
      for(int k=0; k<n; k++) {
        for(int l=0; l<2*n-2; l++) {
          jj=1; if((l==0) && (i!=j)) jj=0; // If l=0, start from home
          if((k==j) && (i!=j)) jj=0;       // No edges from j to j for i
          kk=0; if(l==2*n-3) { 
            kk=D[k][i];
          }
          V[i][j][k][l]=IP.add_variable(D[j][k]+kk,0,jj,Vartype_Integer);
        }
      }
    }
  }

  // is home of i = sum j drives to i
  for(int i=0; i<n; i++) {
    for(int l=0; l<2*n-2; l++) {
      row r;
      for(int j=0; j<n; j++) {
        for(int k=0; k<n; k++) if(k!=i) {
          r+=(row) V[k][j][i][l];
        }
        r-=(row) V[i][j][i][l];
      }
      c=IP.add_basic_constraint(r==0);
      //CPC->add_sym_constraint(new cp_basic_constraint(c));
    }
  };

  // i visits k exactly once
  for(int i=0; i<n; i++) {
    for(int k=0; k<n; k++) if(k!=i) {
      row r;
      for(int j=0; j<n; j++) {
        for(int l=0; l<2*n-2; l++) {
          r+=(row) V[i][j][k][l];
        }
      }
      c=IP.add_basic_constraint(r==1);
      //CPC->add_sym_constraint(new cp_basic_constraint(c));
    }
  }

  // inflow = outflow
  for(int i=0; i<n; i++) {
    for(int j=0; j<n; j++) {
      for(int l=1; l<2*n-2; l++) {
        row r;
        for(int k=0; k<n; k++) {
          r+=(row) V[i][j][k][l];
        };
        for(int k=0; k<n; k++) {
          r-=(row) V[i][k][j][l-1];
        };
        c=IP.add_basic_constraint(r==0);
        //CPC->add_sym_constraint(new cp_basic_constraint(c));
      }
    }
  }

  // every flow is one
  for(int i=0; i<n; i++) {
    row r;
    for(int j=0; j<n; j++) {
      for(int k=0; k<n; k++) {
        r+=(row) V[i][j][k][0];
      }
    }
    c=IP.add_basic_constraint(r==1);
    //CPC->add_sym_constraint(new cp_basic_constraint(c));
  };
  // in the first U+1 steps everyone leaves home
  for(int i=0; i<n; i++) {
    row r;
    for(int j=0; j<n; j++) if(j!=i) {
      for(int l=0; l<U+1; l++) {
        r+=(row) V[i][i][j][l];
      }
    };
    c=IP.add_basic_constraint(r>=1);
    //CPC->add_sym_constraint(new cp_basic_constraint(c));
  };
  /*
  // in every interval of length 2U everyone leaves home and everyone comes home
  for(int l=1; l<2*n-2-2*U+1; l++) {
    for(int i=0; i<n; i++) {
      row R, R1;
      for(int l1=0; l1<2*U; l1++) {
        for(int k=0; k<n; k++) if(k!=i) {
          R+=(row) V[i][i][k][l+l1];
          R1+=(row) V[i][k][i][l+l1];
        }
      }
      c=IP.add_basic__constraint(R>=1);
      //CPC->add_sym_constraint(new cp_basic_constraint(c));
      c=IP.add_basic__constraint(R1>=1);
      //CPC->add_sym_constraint(new cp_basic_constraint(c));
    };
  }
  // in every interval of length U+1 at least n/2 leave home
  for(int l=0; l<2*n-2-U; l++) {
    row r;
    for(int i=0; i<n; i++) {
      for(int j=0; j<n; j++) if(j!=l) {
        for(int l1=0; l1<U+1; l1++) {
          r+=(row) V[i][i][j][l+l1];
        }
      }
    }
    IP.add_basic__constraint(r>=n/2);
    //CPC->add_sym_constraint(new cp_basic_constraint(c));
  };
  */
  // in every interval of length U+1, at least one edge home and one edge away.
  // TODO does not work
  /*
  for(int i=0; i<n; i++) {
    for(int l=0; l<2*n-2-U-1; l++) {
      row r1, r2;
      for(int j=0; j<n; j++) {
        for(int l1=l; l1<l+U+1; l1++) {
          if(l1!=2*n-2)
            r1+=(row) V[i][i][j][l1];
          if(j!=i) {
            for(int k=0; k<n; k++) {
              r2+=(row) V[i][j][k][l1];
            }
          }
        }
      }
      c=IP.add_basic__constraint(r1>=1);
      //CPC->add_sym_constraint(new cp_basic_constraint(c));
      c=IP.add_basic__constraint(r2>=1);
      //CPC->add_sym_constraint(new cp_basic_constraint(c));
    }
  };
  */
  /*
  for(int i=0; i<n; i++) {
    for(int l=0; l<2*n-2-U; l++) {
      row r;
      for(int j=0; j<n; j++) if(j!=i) {
        for(int l1=l; l1<l+U+1; l1++) {
          if(l1!=2*n-2) {
            r+=(row) V[i][i][j][l1];
            r+=(row) V[i][j][i][l1];
          }
        }
      }
      c=IP.add_basic__constraint(r>=1);
      //CPC->add_sym_constraint(new cp_basic_constraint(c));
    }
  };
  */

  // some clieqes
  for(int i=0; i<n; i++) {
    for(int l=0; l<2*n-4; l++) {
      for(int j=0; j<n; j++) if(j!=i) {
        row r;
        r+=(row) V[i][i][i][l+1];
        for(int k=0; k<n; k++) if(k!=j) {
          r+=(row) V[i][k][j][l];
          if(k!=i) r+=(row) V[i][j][k][l];
          if(k!=i) r+=(row) V[i][k][j][l+2];
          r+=(row) V[i][j][k][l+2];
        };
        c=IP.add_basic_constraint(r<=1);
        //CPC->add_sym_constraint(new cp_basic_constraint(c));
      };
    };
  };

  /*
  for(int i=0; i<n; i++) {
    for(int l=0; l<2*n-3; l++) {
      for(int j=0; j<n; j++) if(j!=i) {
        for(int k=0; k<n; k++) if((k!=i) && (k!=j)) {
          row r;
          r+=V[i][j][k][l];
          r+=V[i][j][k][l+1];
          r+=V[i][k][j][l];
          r+=V[i][k][j][l+1];
          for(int k1=0; k1<n; k1++) if((k1!=k) && (k1!=j)) {
            for(int j1=0; j1<n; j1++) if((j1!=k) && (j1!=j)) {
              r+=V[i][j1][k1][l];
            }
          }
          c=IP.add_basic__constraint(r<=1);
          //CPC->add_sym_constraint(new cp_basic_constraint(c));
        }
      }
    }
  };
  */
  /*
  for(int i=0; i<n; i++) {
    for(int l=1; l<2*n-2-U-2; l++) {
      row r;
      for(int l1=0; l1<U+1; l1++) {
        r+=V[i][i][i][l+l1];
        for(int j=0; j<n; j++) if(j!=i) {
          for(int k=0; k<n; k++) if(k!=i) {
            r+V[i][j][k][l+l1];
          }
        }
      }
      for(int j=0; j<n; j++) if(j!=i) {
        r-=V[i][j][i][l+U+1];
        r-=V[i][i][j][l+U+1];
        r-=V[i][j][i][l-1];
        r-=V[i][i][j][l-1];
      };
      IP.add_basic__constraint(r<=U-2);
    }
  };
  */
  // Symmetry elimination.
  c=IP.add_basic_constraint(((row) V[0][0][0][0])-V[0][0][0][2*n-3]<=0);
  //CPC->add_sym_constraint(new cp_basic_constraint(c));
  /*
  for(int i=0; i<n; i++) {
    for(int l=0; l<2*n-3; l++) {
      for(int j=0; j<n; j++) if(j!=i) {
        for(int k=0; k<n; k++) if((k!=i) && (k!=j)) {
          for(int m=0; m<n; m++) if((m!=i) && (m!=j) && (m!=k)) {
            row r;
            r+=V[i][j][k][l];
            r+=V[i][j][k][l+1];
            r+=V[i][k][j][l];
            r+=V[i][k][j][l+1];
            r+=V[i][m][i][l];
            r+=V[i][i][m][l+1];
            c=IP.add_basic__constraint(r<=1);
            //CPC->add_sym_constraint(new cp_basic_constraint(c));
          }
        }
      }
    }
  };
  */

  // Entweger von j nach k oder von k nach j, nicht beides
  /*
  for(int i=0; i<n; i++) {
    for(int j=0; j<n; j++) if(j!=i) {
      for(int k=0; k<n; k++) if((k!=i) && (k!=j)) {
        row r;
        for(int l=0; l<2*n-2; l++) {
          r+=V[i][j][k][l];
          r+=V[i][k][j][l];
        }
        c=IP.add_basic__constraint(r<=1);
        //CPC->add_sym_constraint(new cp_basic_constraint(c));
      }
    }
  }
  */
  /* TODO durcheinander
  for(int l=1; l<2*n-1-(n/2-1); l++) {
    row r;
    for(int i=0; i<n; i++) {
      for(int j=0; j<n; j++) if(j!=i) {
        for(int l1=l; l1<l+n/2-1; l1++) {
          r+=V[i][i][j][l1];
        }
      }
    }
    IP.add_basic__constraint(r>=1);
  }
  */
  /*
  graph G;
  node N[n];
  edge E[n][n];
  for(int i=0; i<n; i++) {
    N[i]=G.new_node();
  }
  for(int i=0; i<n; i++) for(int j=i+1; j<n; j++) {
    E[i][j]=G.new_edge(N[i],N[j]);
  };

  var_map<edge> VM[n];
  for(int i=0; i<n; i++) {
    for(int j=0; j<n; j++) {
      for(int k=j+1; k<n; k++) {
        var v=IP.add_variable(0,0,2*n,Vartype_Integer);
        VM[i][E[j][k]]=v;
        row r;
        for(int l=0; l<2*n-2; l++) {
          r+=V[i][j][k][l];
          r+=V[i][k][j][l];
        }
        if(j==i) {
          for(int m=0; m<n; m++) {
            r+=V[i][m][k][2*n-3];
          }
        }
        if(k==i) {
          for(int m=0; m<n; m++) {
            r+=V[i][m][j][2*n-3];
          };
        };
        IP.add_basic__constraint(r==v);
      };
    };
    IP.add_sym_constraint(new TOUR(G, VM[i]));
  };
  */

  // at least one over every cut

  int pw2[n];
  pw2[0]=1;
  for(int i=1; i<n; i++) pw2[i]=2*pw2[i-1];
  for(int i=0; i<n; i++) {
    cout<<i<<" Start\n";
    for(int j=1; j<pw2[n-1]; j++) {
      int k=0;
      for(int l=0; l<n; l++) if((j & pw2[l])==pw2[l]) k++;
      if((k>1) && (k<n-1)) {
        row r;
        for(int j1=0; j1<n; j1++) if((j & pw2[j1])==0) {
          for(int k1=0; k1<n; k1++) if((j & pw2[k1])==pw2[k1]) {
            for(int l1=0; l1<2*n-2; l1++) {
              r+=(row) V[i][j1][k1][l1];
            }
          }
        }
        c=IP.add_basic_constraint(r >=1);
        //CPC->add_sym_constraint(new cp_basic_constraint(c));
      }
    }
  }

  //CPC->add_sym_constraint(new rrt(n, V));

  //IP.add_sym_constraint(CPC);

  IP.optimize();

  cout<<"Binaries\n";
  for(int i=0; i<n*n*n*(2*n-2); i++) {
    cout<<" x"<<i;
    if(i%10==0) cout<<endl;
  }
  cout<<"\nEnd\n";


  for(int i=0; i<n; i++) {
    cout<<endl<<i<<": ";
    for(int l=0; l<2*n-2; l++) {
      for(int j=0; j<n; j++) {
        double dd=0;
        for(int k=0; k<n; k++) {
          
          if(IP.get_solution(V[i][k][j][l])>0.01) 
            dd+=IP.get_solution(V[i][k][j][l]);
            //cout<<l<<" "<<j<<" "<<k<<" "<<IP.get_solution(V[i][j][k][l])<<endl;
        }
        if(dd>0.01) cout<<l<<" "<<j<<" "<<dd<<endl;
      }
    }
  }

  cout<<endl;
}

---