sudoku.cpp

// Copyright (C) 2005, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

#include <cassert>
#include <iomanip>

#include "CoinPragma.hpp"
// For Branch and bound
#include "OsiSolverInterface.hpp"
#include "CbcModel.hpp"
#include "CoinModel.hpp"
// For all different
#include "CbcBranchCut.hpp"
#include "CbcBranchActual.hpp"
#include "CbcBranchAllDifferent.hpp"
#include "CbcCutGenerator.hpp"
#include "CglAllDifferent.hpp"
#include "OsiClpSolverInterface.hpp"
#include "CglStored.hpp"

#include "CoinTime.hpp"

/************************************************************************

This shows how we can define a new branching method to solve problems with
all different constraints.

We are going to solve a sudoku problem such as

1, , ,4, , ,7, ,
 ,2, , ,5, , ,8,
8,7,3, , ,6, , ,9
4, , ,7, , ,1, ,
 ,5, , ,8, , ,2,
 , ,6, ,4,9, , ,3
7, , ,1, , ,4, ,
 ,8, ,6,2, , ,5,
 , ,9, ,7,3, ,1,6

The input should be exported from spreadsheet as a csv file where cells are
empty unless value is known.

We set up a fake objective and simple constraints to say sum must be 45

and then we add all different branching (CbcBranchAllDifferent) 
and all different cuts (to fix variables) (CglAllDifferent).

CbcBranchAllDifferent is really just an example of a cut branch.  If we wish to stop x==y
then we can have two branches - one x <= y-1 and the other x >= y+1.  It should be easy for the user to 
make up similar cut branches for other uses.

Note - this is all we need to solve most 9 x 9 puzzles because they seem to solve
at root node or easily.  To solve 16 x 16 puzzles we need more.  All different cuts
need general integer variables to be fixed while we can branch so they are just at bounds.  
To get round that we can introduce extra 0-1 variables such that general integer x = sum j * delta j
and then do N way branching on these (CbcNWay) so that we fix one delta j to 1.  At the same time we use the
new class CbcConsequence (used in N way branching) which when delta j goes to 1 fixes other variables.
So it will fix x to the correct value and while we are at it we can fix some delta variables in other
sets to zero (as per all different rules).  Finally as well as storing the instructions which say if 
delta 11 is 1 then delta 21 is 0 we can also add this in as a cut using new trivial cut class
CglStored.

************************************************************************/

int main(int argc, const char *argv[])
{
  // Get data
  std::string fileName = "./sudoku_sample.csv";
  if (argc >= 2)
    fileName = argv[1];
  FILE *fp = fopen(fileName.c_str(), "r");
  if (!fp) {
    printf("Unable to open file %s\n", fileName.c_str());
    exit(0);
  }
#define MAX_SIZE 16
  int valueOffset = 1;
  double lo[MAX_SIZE * MAX_SIZE], up[MAX_SIZE * MAX_SIZE];
  char line[80];
  int row, column;
  /*************************************** 
     Read .csv file and see if 9 or 16 Su Doku
  ***************************************/
  int size = 9;
  for (row = 0; row < size; row++) {
    fgets(line, 80, fp);
    // Get size of sudoku puzzle (9 or 16)
    if (!row) {
      int get = 0;
      size = 1;
      while (line[get] >= 32) {
        if (line[get] == ',')
          size++;
        get++;
      }
      assert(size == 9 || size == 16);
      printf("Solving Su Doku of size %d\n", size);
      if (size == 16)
        valueOffset = 0;
    }
    int get = 0;
    for (column = 0; column < size; column++) {
      lo[size * row + column] = valueOffset;
      up[size * row + column] = valueOffset - 1 + size;
      if (line[get] != ',' && line[get] >= 32) {
        // skip blanks
        if (line[get] == ' ') {
          get++;
          continue;
        }
        int value = line[get] - '0';
        if (size == 9) {
          assert(value >= 1 && value <= 9);
        } else {
          assert(size == 16);
          if (value < 0 || value > 9) {
            if (line[get] == '"') {
              get++;
              value = 10 + line[get] - 'A';
              if (value < 10 || value > 15) {
                value = 10 + line[get] - 'a';
              }
              get++;
            } else {
              value = 10 + line[get] - 'A';
              if (value < 10 || value > 15) {
                value = 10 + line[get] - 'a';
              }
            }
          }
          assert(value >= 0 && value <= 15);
        }
        lo[size * row + column] = value;
        up[size * row + column] = value;
        get++;
      }
      get++;
    }
  }
  int block_size = (int)sqrt((double)size);
  /*************************************** 
     Now build rules for all different
     3*9 or 3*16 sets of variables
     Number variables by row*size+column
  ***************************************/
  int starts[3 * MAX_SIZE + 1];
  int which[3 * MAX_SIZE * MAX_SIZE];
  int put = 0;
  int set = 0;
  starts[0] = 0;
  // By row
  for (row = 0; row < size; row++) {
    for (column = 0; column < size; column++)
      which[put++] = row * size + column;
    starts[set + 1] = put;
    set++;
  }
  // By column
  for (column = 0; column < size; column++) {
    for (row = 0; row < size; row++)
      which[put++] = row * size + column;
    starts[set + 1] = put;
    set++;
  }
  // By box
  for (row = 0; row < size; row += block_size) {
    for (column = 0; column < size; column += block_size) {
      for (int row2 = row; row2 < row + block_size; row2++) {
        for (int column2 = column; column2 < column + block_size; column2++)
          which[put++] = row2 * size + column2;
      }
      starts[set + 1] = put;
      set++;
    }
  }
  OsiClpSolverInterface solver1;

  /*************************************** 
     Create model
     Set variables to be general integer variables although
     priorities probably mean that it won't matter
  ***************************************/
  CoinModel build;
  // Columns
  char name[4];
  for (row = 0; row < size; row++) {
    for (column = 0; column < size; column++) {
      if (row < 10) {
        if (column < 10)
          sprintf(name, "X%d%d", row, column);
        else
          sprintf(name, "X%d%c", row, 'A' + (column - 10));
      } else {
        if (column < 10)
          sprintf(name, "X%c%d", 'A' + (row - 10), column);
        else
          sprintf(name, "X%c%c", 'A' + (row - 10), 'A' + (column - 10));
      }
      double value = CoinDrand48() * 100.0;
      build.addColumn(0, NULL, NULL, lo[size * row + column],
        up[size * row + column], value, name, true);
    }
  }
  /*************************************** 
     Now add in extra variables for N way branching
  ***************************************/
  for (row = 0; row < size; row++) {
    for (column = 0; column < size; column++) {
      int iColumn = size * row + column;
      double value = lo[iColumn];
      if (value < up[iColumn]) {
        for (int i = 0; i < size; i++)
          build.addColumn(0, NULL, NULL, 0.0, 1.0, 0.0);
      } else {
        // fixed
        // obviously could do better  if we missed out variables
        int which = ((int)value) - valueOffset;
        for (int i = 0; i < size; i++) {
          if (i != which)
            build.addColumn(0, NULL, NULL, 0.0, 0.0, 0.0);
          else
            build.addColumn(0, NULL, NULL, 0.0, 1.0, 0.0);
        }
      }
    }
  }

  /*************************************** 
     Now rows
  ***************************************/
  double values[] = { 1.0, 1.0, 1.0, 1.0,
    1.0, 1.0, 1.0, 1.0,
    1.0, 1.0, 1.0, 1.0,
    1.0, 1.0, 1.0, 1.0 };
  int indices[MAX_SIZE + 1];
  double rhs = size == 9 ? 45.0 : 120.0;
  for (row = 0; row < 3 * size; row++) {
    int iStart = starts[row];
    for (column = 0; column < size; column++)
      indices[column] = which[column + iStart];
    build.addRow(size, indices, values, rhs, rhs);
  }
  double values2[MAX_SIZE + 1];
  values2[0] = -1.0;
  for (row = 0; row < size; row++)
    values2[row + 1] = row + valueOffset;
  // Now add rows for extra variables
  for (row = 0; row < size; row++) {
    for (column = 0; column < size; column++) {
      int iColumn = row * size + column;
      int base = size * size + iColumn * size;
      indices[0] = iColumn;
      for (int i = 0; i < size; i++)
        indices[i + 1] = base + i;
      build.addRow(size + 1, indices, values2, 0.0, 0.0);
    }
  }
  solver1.loadFromCoinModel(build);
  build.writeMps("xx.mps");

  double time1 = CoinCpuTime();
  solver1.initialSolve();
  CbcModel model(solver1);
  model.solver()->setHintParam(OsiDoReducePrint, true, OsiHintTry);
  model.solver()->setHintParam(OsiDoScale, false, OsiHintTry);
  /*************************************** 
    Add in All different cut generator and All different branching
    So we will have integers then cut branching then N way branching
    in reverse priority order
  ***************************************/

  // Cut generator
  CglAllDifferent allDifferent(3 * size, starts, which);
  model.addCutGenerator(&allDifferent, -99, "allDifferent");
  model.cutGenerator(0)->setWhatDepth(5);
  CbcObject **objects = new CbcObject *[4 * size * size];
  int nObj = 0;
  for (row = 0; row < 3 * size; row++) {
    int iStart = starts[row];
    objects[row] = new CbcBranchAllDifferent(&model, size, which + iStart);
    objects[row]->setPriority(2000 + nObj); // do after rest satisfied
    nObj++;
  }
  /*************************************** 
     Add in N way branching and while we are at it add in cuts
  ***************************************/
  CglStored stored;
  for (row = 0; row < size; row++) {
    for (column = 0; column < size; column++) {
      int iColumn = row * size + column;
      int base = size * size + iColumn * size;
      int i;
      for (i = 0; i < size; i++)
        indices[i] = base + i;
      CbcNWay *obj = new CbcNWay(&model, size, indices, nObj);
      int seq[200];
      int newUpper[200];
      memset(newUpper, 0, sizeof(newUpper));
      for (i = 0; i < size; i++) {
        int state = 9999;
        int one = 1;
        int nFix = 1;
        // Fix real variable
        seq[0] = iColumn;
        newUpper[0] = valueOffset + i;
        int kColumn = base + i;
        int j;
        // same row
        for (j = 0; j < size; j++) {
          int jColumn = row * size + j;
          int jjColumn = size * size + jColumn * size + i;
          if (jjColumn != kColumn) {
            seq[nFix++] = jjColumn; // must be zero
          }
        }
        // same column
        for (j = 0; j < size; j++) {
          int jColumn = j * size + column;
          int jjColumn = size * size + jColumn * size + i;
          if (jjColumn != kColumn) {
            seq[nFix++] = jjColumn; // must be zero
          }
        }
        // same block
        int kRow = row / block_size;
        kRow *= block_size;
        int kCol = column / block_size;
        kCol *= block_size;
        for (j = kRow; j < kRow + block_size; j++) {
          for (int jc = kCol; jc < kCol + block_size; jc++) {
            int jColumn = j * size + jc;
            int jjColumn = size * size + jColumn * size + i;
            if (jjColumn != kColumn) {
              seq[nFix++] = jjColumn; // must be zero
            }
          }
        }
        // seem to need following?
        const int *upperAddress = newUpper;
        const int *seqAddress = seq;
        CbcFixVariable fix(1, &state, &one, &upperAddress, &seqAddress, &nFix, &upperAddress, &seqAddress);
        obj->setConsequence(indices[i], fix);
        // Now do as cuts
        for (int kk = 1; kk < nFix; kk++) {
          int jColumn = seq[kk];
          int cutInd[2];
          cutInd[0] = kColumn;
          if (jColumn > kColumn) {
            cutInd[1] = jColumn;
            stored.addCut(-COIN_DBL_MAX, 1.0, 2, cutInd, values);
          }
        }
      }
      objects[nObj] = obj;
      objects[nObj]->setPriority(nObj);
      nObj++;
    }
  }
  model.addObjects(nObj, objects);
  for (row = 0; row < nObj; row++)
    delete objects[row];
  delete[] objects;
  model.messageHandler()->setLogLevel(1);
  model.addCutGenerator(&stored, 1, "stored");
  // Say we want timings
  int numberGenerators = model.numberCutGenerators();
  int iGenerator;
  for (iGenerator = 0; iGenerator < numberGenerators; iGenerator++) {
    CbcCutGenerator *generator = model.cutGenerator(iGenerator);
    generator->setTiming(true);
  }
  // Set this to get all solutions (all ones in newspapers should only have one)
  //model.setCutoffIncrement(-1.0e6);
  /*************************************** 
     Do branch and bound
  ***************************************/
  // Do complete search
  model.branchAndBound();
  std::cout << "took " << CoinCpuTime() - time1 << " seconds, "
            << model.getNodeCount() << " nodes with objective "
            << model.getObjValue()
            << (!model.status() ? " Finished" : " Not finished")
            << std::endl;

  /*************************************** 
     Print solution and check it is feasible
     We could modify output so could be imported by spreadsheet
  ***************************************/
  if (model.getMinimizationObjValue() < 1.0e50) {

    const double *solution = model.bestSolution();
    int put = 0;
    for (row = 0; row < size; row++) {
      for (column = 0; column < size; column++) {
        int value = (int)floor(solution[row * size + column] + 0.5);
        assert(value >= lo[put] && value <= up[put]);
        // save for later test
        lo[put++] = value;
        printf("%d ", value);
      }
      printf("\n");
    }
    // check valid
    bool valid = true;
    // By row
    for (row = 0; row < size; row++) {
      put = 0;
      for (column = 0; column < size; column++)
        which[put++] = row * size + column;
      assert(put == size);
      int i;
      for (i = 0; i < put; i++)
        which[i] = (int)lo[which[i]];
      std::sort(which, which + put);
      int last = valueOffset - 1;
      for (i = 0; i < put; i++) {
        int value = which[i];
        if (value != last + 1)
          valid = false;
        last = value;
      }
    }
    // By column
    for (column = 0; column < size; column++) {
      put = 0;
      for (row = 0; row < size; row++)
        which[put++] = row * size + column;
      assert(put == size);
      int i;
      for (i = 0; i < put; i++)
        which[i] = (int)lo[which[i]];
      std::sort(which, which + put);
      int last = valueOffset - 1;
      for (i = 0; i < put; i++) {
        int value = which[i];
        if (value != last + 1)
          valid = false;
        last = value;
      }
    }
    // By box
    for (row = 0; row < size; row += block_size) {
      for (column = 0; column < size; column += block_size) {
        put = 0;
        for (int row2 = row; row2 < row + block_size; row2++) {
          for (int column2 = column; column2 < column + block_size; column2++)
            which[put++] = row2 * size + column2;
        }
        assert(put == size);
        int i;
        for (i = 0; i < put; i++)
          which[i] = (int)lo[which[i]];
        std::sort(which, which + put);
        int last = valueOffset - 1;
        for (i = 0; i < put; i++) {
          int value = which[i];
          if (value != last + 1)
            valid = false;
          last = value;
        }
      }
    }
    if (valid) {
      printf("solution is valid\n");
    } else {
      printf("solution is not valid\n");
      abort();
    }
  }
  return 0;
}