walksat.c

sp组合优化算法的一个改进算法

/* walksat version 35 */
/* version 1 by Bram Cohen 7/93 */
/* versions 2 - 35  by Henry Kautz */

/************************************/
/* Compilation flags                */
/************************************/

/* If the constant HUGE is set, then dynamic arrays are used instead of static arrays.
   This allows very large or very small problems to be handled, without
   having to adjust the constants MAXATOM and/or MAXCLAUSE.  However, on some
   architectures (e.g. SGI) the compiled program is about 25% slower, because not
   as many optimizations can be performed. */

#define HUGE 1

/* If the constant NT is set, then the program is modified to compile under Windows NT.
   Currently, this eliminates the timing functionality. */

/* #define NT 1 */

/************************************/
/* Standard includes                */
/************************************/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <sys/types.h>
#include <limits.h>
#include <signal.h>

#ifndef NT
#include <sys/times.h>
#include <sys/time.h>
#endif

#ifndef CLK_TCK
#define CLK_TCK 60
#endif

#ifdef NT
#define random() rand()
#define srandom(seed) srand(seed)
#endif

/************************************/
/* Constant parameters              */
/************************************/

#define MAXATOM 1000000		/* maximum possible number of atoms */

#ifdef Huge
#define STOREBLOCK 2000000	/* size of block to malloc each time */
#else
#define STOREBLOCK 2000000	/* size of block to malloc each time */
#define MAXCLAUSE 500000	/* maximum possible number of clauses */
#endif

#define TRUE 1
#define FALSE 0

#define MAXLENGTH 500           /* maximum number of literals which can be in any clause */


/************************************/
/* Internal constants               */
/************************************/

enum heuristics { RANDOM, BEST, TABU, NOVELTY, RNOVELTY };

#define NOVALUE -1
#define INIT_PARTIAL 1
#define HISTMAX 101		/* length of histogram of tail */

#define Var(CLAUSE, POSITION) (ABS(clause[CLAUSE][POSITION]))

static int scratch;
#define ABS(x) ((scratch=(x))>0?(scratch):(-scratch))

#define BIG 100000000



/************************************/
/* Main data structures             */
/************************************/

/* Atoms start at 1 */
/* Not a is recorded as -1 * a */
/* One dimensional arrays are statically allocated. */
/* Two dimensional arrays are dynamically allocated in */
/* the second dimension only.  */

int numatom;
int numclause;
int numliterals;

#ifdef Huge
int ** clause;			/* clauses to be satisfied */
				/* indexed as clause[clause_num][literal_num] */
int * size;			/* length of each clause */
int * false;			/* clauses which are false */
int * lowfalse;
int * wherefalse;		/* where each clause is listed in false */
int * numtruelit;		/* number of true literals in each clause */
#else
int * clause[MAXCLAUSE];	/* clauses to be satisfied */
				/* indexed as clause[clause_num][literal_num] */
int size[MAXCLAUSE];		/* length of each clause */
int false[MAXCLAUSE];		/* clauses which are false */
int lowfalse[MAXCLAUSE];
int wherefalse[MAXCLAUSE];	/* where each clause is listed in false */
int numtruelit[MAXCLAUSE];	/* number of true literals in each clause */
#endif

int *occurence[2*MAXATOM+1];	/* where each literal occurs */
				/* indexed as occurence[literal+MAXATOM][occurence_num] */

int numoccurence[2*MAXATOM+1];	/* number of times each literal occurs */


int atom[MAXATOM+1];		/* value of each atom */ 
int lowatom[MAXATOM+1];
int solution[MAXATOM+1];

int changed[MAXATOM+1];		/* step at which atom was last flipped */

int breakcount[MAXATOM+1];	/* number of clauses that become unsat if var if flipped */
int makecount[MAXATOM+1];	/* number of clauses that become sat if var if flipped */

int numfalse;			/* number of false clauses */

/************************************/
/* Global flags and parameters      */
/************************************/

int abort_flag;

int heuristic = BEST;		/* heuristic to be used */
int numerator = NOVALUE;	/* make random flip with numerator/denominator frequency */
int denominator = 100;		
int tabu_length;		/* length of tabu list */

long int numflip;		/* number of changes so far */
long int numnullflip;		/*  number of times a clause was picked, but no  */
				/*  variable from it was flipped  */
int numrun = 10;
int cutoff = 100000;
int base_cutoff = 100000;
int target = 0;
int numtry = 0;			/* total attempts at solutions */
int numsol = NOVALUE;		/* stop after this many tries succeeds */
int superlinear = FALSE;

int makeflag = FALSE;		/* set to true by heuristics that require the make values to be calculated */

/* Histogram of tail */

long int tailhist[HISTMAX];	/* histogram of num unsat in tail of run */
long histtotal;
int tail = 3;
int tail_start_flip;

/* Printing options */

int printonlysol = FALSE;
int printsolcnf = FALSE;
int printfalse = FALSE;
int printlow = FALSE;
int printhist = FALSE;
int printtrace = FALSE;
int trace_assign = FALSE;

/* Initialization options */

char initfile[100] = { 0 };
int initoptions = FALSE;

/* Randomization */

int seed;			/* seed for random */
struct timeval tv;
struct timezone tzp;

/* Statistics */

double expertime;
long flips_this_solution;
long int lowbad;		/* lowest number of bad clauses during try */
long int totalflip = 0;		/* total number of flips in all tries so far */
long int totalsuccessflip = 0;	/* total number of flips in all tries which succeeded so far */
int numsuccesstry = 0;		/* total found solutions */

long x;
long integer_sum_x = 0;
double sum_x = 0.0;
double sum_x_squared = 0.0;
double mean_x;
double second_moment_x;
double variance_x;
double std_dev_x;
double std_error_mean_x;
double seconds_per_flip;
int r;
int sum_r = 0;
double sum_r_squared = 0.0;
double mean_r;
double variance_r;
double std_dev_r;
double std_error_mean_r;

double avgfalse;
double sumfalse;
double sumfalse_squared;
double second_moment_avgfalse, variance_avgfalse, std_dev_avgfalse, ratio_avgfalse;
double std_dev_avgfalse;
double f;
double sample_size;

double sum_avgfalse = 0.0;
double sum_std_dev_avgfalse = 0.0;
double mean_avgfalse;
double mean_std_dev_avgfalse;
int number_sampled_runs = 0;
double ratio_mean_avgfalse;

double suc_sum_avgfalse = 0.0;
double suc_sum_std_dev_avgfalse = 0.0;
double suc_mean_avgfalse;
double suc_mean_std_dev_avgfalse;
int suc_number_sampled_runs = 0;
double suc_ratio_mean_avgfalse;

double nonsuc_sum_avgfalse = 0.0;
double nonsuc_sum_std_dev_avgfalse = 0.0;
double nonsuc_mean_avgfalse;
double nonsuc_mean_std_dev_avgfalse;
int nonsuc_number_sampled_runs = 0;
double nonsuc_ratio_mean_avgfalse;

/* Hamming calcualations */

char hamming_target_file[512] = { 0 };
char hamming_data_file[512] = { 0 };
int hamming_sample_freq;
int hamming_flag = FALSE;
int hamming_distance;
int hamming_target[MAXATOM+1];
void read_hamming_file(char initfile[]);
void open_hamming_data(char initfile[]);
int calc_hamming_dist(int atom[], int hamming_target[], int numatom);
FILE * hamming_fp;

/* Noise level */
int samplefreq = 1;

/************************************/
/* Forward declarations             */
/************************************/

void parse_parameters(int argc,char *argv[]);
void print_parameters(int argc, char * argv[]);

int pickrandom(void);
int pickbest(void);
int picktabu(void);
int picknovelty(void);
int pickrnovelty(void);
char * heuristic_names[] = { "random", "best", "tabu", "novelty", "rnovelty" };
static int (*pickcode[])(void) =
      {pickrandom, pickbest, picktabu, picknovelty, pickrnovelty};

double elapsed_seconds(void);
int countunsat(void);
void scanone(int argc, char *argv[], int i, int *varptr);
void init(char initfile[], int initoptions);
void initprob(void);                 /* create a new problem */
void flipatom(int toflip);           /* changes the assignment of the given literal */

void print_false_clauses(long int lowbad);
void save_false_clauses(long int lowbad);
void print_low_assign(long int lowbad);
void save_low_assign(void);
void save_solution(void);

void print_current_assign(void);
void handle_interrupt(int sig);

long super(int i);


void print_statistics_header(void);
void initialize_statistics(void);
void update_statistics_start_try(void);
void print_statistics_start_flip(void);
void update_and_print_statistics_end_try(void);
void update_statistics_end_flip(void);
void print_statistics_final(void);
void print_sol_cnf(void);

/************************************/
/* Main                             */
/************************************/


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

    gettimeofday(&tv,&tzp);
    seed = (( tv.tv_sec & 0177 ) * 1000000) + tv.tv_usec;
    parse_parameters(argc, argv);
    srandom(seed);
    print_parameters(argc, argv);
    initprob();
    initialize_statistics();
    print_statistics_header();
    signal(SIGINT, (void *) handle_interrupt);
    abort_flag = FALSE;
    (void) elapsed_seconds();

    while (! abort_flag && numsuccesstry < numsol && numtry < numrun) {
	numtry++;
	init(initfile, initoptions);
	update_statistics_start_try();
	numflip = 0;
	  
	if (superlinear) cutoff = base_cutoff * super(numtry);

	while((numfalse > target) && (numflip < cutoff)) {
	    print_statistics_start_flip();
	    numflip++;
	    flipatom((pickcode[heuristic])());
	    update_statistics_end_flip();
	}
	update_and_print_statistics_end_try();
    }
    expertime = elapsed_seconds();
    print_statistics_final();
    return 0;
}



void parse_parameters(int argc,char *argv[])
{
  int i;
  int temp;

  for (i=1;i < argc;i++)
    {
      if (strcmp(argv[i],"-seed") == 0)
	scanone(argc,argv,++i,&seed);
      else if (strcmp(argv[i],"-hist") == 0)
	printhist = TRUE;
      else if (strcmp(argv[i],"-cutoff") == 0)
	scanone(argc,argv,++i,&cutoff);
      else if (strcmp(argv[i],"-random") == 0)
	heuristic = RANDOM;
      else if (strcmp(argv[i],"-novelty") == 0){
	heuristic = NOVELTY;
	makeflag = TRUE;
      }
      else if (strcmp(argv[i],"-rnovelty") == 0){
	heuristic = RNOVELTY;
	makeflag = TRUE;
      }
      else if (strcmp(argv[i],"-best") == 0)
	heuristic = BEST;
      else if (strcmp(argv[i],"-noise") == 0){
	scanone(argc,argv,++i,&numerator);
	if (i < argc-1 && sscanf(argv[i+1],"%i",&temp)==1){
	  denominator = temp;
	  i++;
	}
      }
      else if (strcmp(argv[i],"-init") == 0  && i < argc-1)
	sscanf(argv[++i], " %s", initfile);
      else if (strcmp(argv[i],"-hamming") == 0  && i < argc-3){
	sscanf(argv[++i], " %s", hamming_target_file);
	sscanf(argv[++i], " %s", hamming_data_file);
	sscanf(argv[++i], " %i", &hamming_sample_freq);
	hamming_flag = TRUE;
	numrun = 1;
      }
      else if (strcmp(argv[i],"-partial") == 0)
	initoptions = INIT_PARTIAL;
      else if (strcmp(argv[i],"-super") == 0)
	superlinear = TRUE;
      else if (strcmp(argv[i],"-tries") == 0)
	scanone(argc,argv,++i,&numrun);
      else if (strcmp(argv[i],"-target") == 0)
	scanone(argc,argv,++i,&target);
      else if (strcmp(argv[i],"-tail") == 0)
	scanone(argc,argv,++i,&tail);
      else if (strcmp(argv[i],"-sample") == 0)
	scanone(argc,argv,++i,&samplefreq);
      else if (strcmp(argv[i],"-tabu") == 0)
	{
	  scanone(argc,argv,++i,&tabu_length);
	  heuristic = TABU;
	}
      else if (strcmp(argv[i],"-low") == 0)
	printlow = TRUE;
      else if (strcmp(argv[i],"-sol") == 0)
	{
	  printonlysol = TRUE;
	  printlow = TRUE;
	}
      else if (strcmp(argv[i],"-solcnf") == 0)
	{
	  printsolcnf = TRUE;
	  if (numsol == NOVALUE) numsol = 1;
	}
      else if (strcmp(argv[i],"-bad") == 0)
	printfalse = TRUE;
      else if (strcmp(argv[i],"-numsol") == 0)
	scanone(argc,argv,++i,&numsol);
      else if (strcmp(argv[i],"-trace") == 0)
	scanone(argc,argv,++i,&printtrace);
      else if (strcmp(argv[i],"-assign") == 0){
	scanone(argc,argv,++i,&printtrace);
	trace_assign = TRUE;
      }
      else 
	{
	  if (strcmp(argv[i],"-help")!=0 && strcmp(argv[i],"-h")!=0 )
	    fprintf(stderr, "Bad argument %s\n", argv[i]);
	  fprintf(stderr, "General parameters:\n");
	  fprintf(stderr, "  -seed N -cutoff N -tries N -target N\n");
	  fprintf(stderr, "  -numsol N = stop after finding N solutions\n");
	  fprintf(stderr, "  -init FILE = set vars not included in FILE to false\n");
	  fprintf(stderr, "  -partial FILE = set vars not included in FILE randomly\n");
	  fprintf(stderr, "Heuristics:\n");
	  fprintf(stderr, "  -random -best -tabu N -novelty -rnovelty\n");
	  fprintf(stderr, "  -noise N or -noise N M (default M = 100)\n");
	  fprintf(stderr, "Printing:\n");
	  fprintf(stderr, "  -hamming TARGET_FILE DATA_FILE SAMPLE_FREQUENCY\n");
	  fprintf(stderr, "  -trace N = print statistics every N flips\n");
	  fprintf(stderr, "  -assign N = print assignments at flip N, 2N, ...\n");
	  fprintf(stderr, "  -sol = print satisfying assignments\n");
	  fprintf(stderr, "  -low = print lowest assignment each try\n");
	  fprintf(stderr, "  -bad = print unsat clauses each try\n");
	  fprintf(stderr, "  -tail N = assume tail begins at nvars*N\n");
	  fprintf(stderr, "  -solcnf = print sat assign in cnf format, and exit\n");
	  fprintf(stderr, "  -sample N = sample noise level every N flips\n");
	  exit(-1);
	}
    }
  base_cutoff = cutoff;
  if (numsol==NOVALUE || numsol>numrun) numsol = numrun;
  if (numerator==NOVALUE){
    switch(heuristic) {
    case BEST:
    case NOVELTY:
    case RNOVELTY:
      numerator = 50;
      break;
    default:
      numerator = 0;
      break;
    }
  }
}


void print_parameters(int argc, char * argv[])
{
  int i;

#ifdef Huge
  printf("walksat version 35 (Huge)\n");
#else
  printf("walksat version 35\n");
#endif
  printf("command line =");
  for (i=0;i < argc;i++){
    printf(" %s", argv[i]);
  }
  printf("\n");
  printf("seed = %i\n",seed);
  printf("cutoff = %i\n",cutoff);
  printf("tries = %i\n",numrun);

  printf("heuristic = ");
  switch(heuristic)
    {
    case TABU:
      printf("tabu %d", tabu_length);
      break;