sp.c

人工智能的一个重要方法

/*
    Copyright 2002 Alfredo Braunstein, Michele Leone, Marc M閦ard, 
                   Martin Weigt and Riccardo Zecchina

    This file is part of SP (Survey Propagation).

    SP is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or 
    (at your option) any later version.

    SP is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with SP; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <unistd.h>
#include "random.h"
#include "formula.h"
#include "sp.h"

#ifdef QUEUE

#include "queue.h"

#endif //QUEUE

//global system vars & structures
struct vstruct *v=NULL; //all per var information
struct clausestruct *clause=NULL; //all per clause information
double *magn=NULL; //spin magnetization list (for fixing ordering)
int *perm=NULL; //permutation, for update ordering
int M=0; //number of clauses
int N=0; //number of variables
int *ncl=NULL; //clause size histogram
int maxlit=0; //maximum clause size
int freespin; //number of unfixed variables
double epsilon=EPSILONDEF; //convergence criterion
int maxconn=0; //maximum connectivity

//auxiliary vars
int *list=NULL;
double *prod=NULL;

//flags & options
double percent=0;
int fixperstep=1;
int ndanger=0;
int generate=0;
int verbose=0;
int fields=0;
int stop=0;
int complexity=0;
FILE *paraplot=NULL;
FILE *replay=NULL;
char *load=NULL;
int iterations=ITERATIONS;
int K=3;
int uselazy=0;
double rho=0;
double norho=1;

#define EPS (1.0e-16)

int main(int argc, char ** argv)
{
	int oldfreespin;
//parse command line
	parsecommandline(argc,argv);
//read or generate a formula
	if(generate) {
		randomformula(K);
	} else if (load) {
		readvarformula(load);
	} else {
		fprintf(stderr, "error: you must specify some formula (-l or -n & -a)\n");
		return 0;
	}
//define fixperstep if based on percent
	if(percent) {
		fixperstep=N*percent/100.0;
	}
//allocate mem for dynamics
	initmem();
//possibly simplify the formula
	oneclauses();
//write the formula
	writeformula(fopen(FORMULA,"w+"));
//pick an initial starting point for the dynamics
	randomize_eta();
//a first normal convergence helps
	if(uselazy)
		sequential_converge();
//fix ndanger balanced spins (useful for restarts)
	while(ndanger--){
		if(!converge())
			return -1;
		build_list(&index_pap);
		oldfreespin=freespin;
		fix_chunk(1);
		printf("fixed %i PAP var (+%i ucp)\n",1,oldfreespin-freespin-1);
	}
//converge and fix fixperstep spins at a time
	while(converge()) {
		if(complexity)
			compute_sigma();
		oldfreespin=freespin;
		build_list(&index_biased);
		if(fields)
			print_fields();
		if(stop)
			return 0;
		fix_chunk(fixperstep);
		if(ncl[0]==M) {
			writesolution(fopen(SOLUTION, "w+"));
			printf("FOUND!\n");
			return 0;
		}
		printf("fixed %i biased var%s (+%i ucp)\n",fixperstep,fixperstep>1?"s":"",oldfreespin-freespin-fixperstep);
		if(verbose)
			print_stats();
	}
	return -1;
}
void oneclauses()
//initially simplify the formula, eliminating all 1-clauses 
{
	int c,var;
	for(c=0; c<M; c++) { 
		if(clause[c].lits==1) {
			var=clause[c].literal[0].var;
			if(v[var].spin==0) {
				printf("eliminating var %i (in 1-clause)\n",var);
				fix(var,clause[c].literal[0].bar?-1:1);
			}
		}
	}
}

int parsecommandline(int argc, char **argv)
//get command line parameters
{
	double alpha=0;
	char c;
	generate=0;
	usrand(1);
	while((c=getopt(argc, argv,
#ifdef QUEUE
			"z"
#endif
			"R:k:cN:M:r:n:m:a:s:d:hf:v%:e:p:l:F/i:"))!=-1) {
		switch (c) {
		case 'R':
			rho=atof(optarg);
			norho=1-rho;
			break;
		case 'z':
			uselazy=1;
			break;
		case 'l':
			load=optarg;
			break;
		case 'p':
			paraplot=fopen(optarg,"w+");
			break;
		case 'r':
			replay=fopen(optarg,"w+");
			break;
		case 'N':
		case 'n':
			N=atoi(optarg);
			break;
		case 'M':
		case 'm':
			M=atoi(optarg);
			break;
		case 'a':
			alpha=atof(optarg);
			break;
		case 'e':
			epsilon=atof(optarg);
			break;
		case 's':
			usrand(atoi(optarg));
			srand(atoi(optarg));
			break;
		case '/':
			stop=1;
			break;
		case 'k':
		case 'K':
			K=atoi(optarg);
			break;
		case 'v':
			verbose++;
			break;
		case 'F':
			fields=1;
			break;
		case 'f':
			fixperstep=atoi(optarg);
			break;
		case '%':
			percent=atof(optarg);
			break;
		case 'd':
			ndanger=atoi(optarg);
			break;
		case 'i':
			iterations=atoi(optarg);
			break;
		case 'c':
			complexity=1;
			break;
		case 'h':
		default:
			fprintf(stderr, "%s [options]\n"
				"\n"
				"  formula\n"
				"\t-n <numvars>\n"
				"\t-m <numclauses>\n"
				"\t-a <alpha>\n"
				"\t-R <rho>\t modified dynamics (0=sp, 1=bp)\n"
				"\t\t\t (real values inbetween may make sense)\n"
				"\t-l <filename>\t reads formula from file\n"
				"  solving\n"
				"\t-d <danger>\t fix this much PAP spins (experimental)\n"
				"\t-f <fixes>\t per step\n"
				"\t-%% <fixes>\t per step (%%)\n"
				"\t-e <error>\t criterion for convergence\n"
				"\t-z \t\t use lazy convergence instead of sequential\n"
#ifndef QUEUE
				"\t\t\t (unavaliable, recompile with -DQUEUE)\n"
#endif
				"\t-i <iter>\t maximum number of iterations until convergence\n"
				"  stats\n"
				"\t-p <filename>\t output a magneticity plot\n"
				"\t-r <filename>\t replay file\n"
				"\t-c \t\t computes complexity\n"
				"\t-F \t\t print fields\n"
				"\t-v \t\t increase verbosity\n"
				"  misc\n"
				"\t-s <seed>\t (0=use time, default=1)\n"
				"\t-/\t\t stop after first convergence\n"
				"\t-h\t\t this help\n",argv[0]);
			exit(-1);
		}
	}
	if(load && !N && !M && !alpha) {
		generate=0;
	} else if(N && alpha && !M) {
		M=N*alpha;
		generate=1;
	} else if(M && alpha && !N) {
		N=M/alpha;
		generate=1;
	} else if(M && N && alpha==0) {
		generate=1;
	} else {
		fprintf(stderr, "error: you have to specify exactly TWO of -n,-m and -a, or -l FILE (and then a formula is read from FILE)\n");
		exit(-1);
	}
	return 0;
}

inline struct weightstruct normalize(struct weightstruct H)
//normalize a triplet
{
	double norm;
	norm=H.m+H.z+H.p;
	H.m/=norm;
	H.p/=norm;
	H.z/=norm;
	return H;
}

int order(void const *a, void const *b)
//order relation for qsort, uses ranking in magn[]
{
	double aux;
	aux=magn[*((int *)b)]-magn[*((int *)a)];
	return aux<0?-1:(aux>0?1:0);
}

double index_biased(struct weightstruct H)
//most biased ranking
{
	return fabs(H.p-H.m);
}
double index_pap(struct weightstruct H)
//most biased with some fuss
{
	return fabs(H.p-H.m)+randreal()*0.1;
}

double index_para(struct weightstruct H)
//least paramagnetic ranking
{
	return H.z;
}

double index_frozen(struct weightstruct H)
//most paramagnetic ranking
{
	return -H.z;
}

double index_best(struct weightstruct H)
//min(H.p,H.m) ranking
{
	return -(H.p>H.m?H.m:H.p);
}

int build_list(double (* index)(struct weightstruct))
//build an ordered list with order *index which is one of index_?
{
	int var,totsites;
	struct weightstruct H;
	double summag;
	double maxmag;
	summag=0;
	totsites=0;
	for(var=1; var<=N; var++) if(v[var].spin==0) {
		H=compute_field(var);
		list[totsites++]=var;
		magn[var]=(*index)(H);
		maxmag=H.p>H.m?H.p:H.m;
		summag+=maxmag;
	}
	qsort(list, totsites, sizeof(int), &order);
	printf("<bias>:%f\n",summag/totsites);
	if(paraplot) {
		fprintf(paraplot,"%f %f %i %i %i\n", (N-freespin)*1.0/N, summag/totsites, freespin, ncl[2], ncl[3]);
		fflush(paraplot);
	}
	if(summag/totsites<PARAMAGNET) {
		printf("paramagnetic state\n");
		printf("sub-formula has:\n");
		print_stats();
		writesolution(fopen(SPSOL,"w+"));
		solvesubsystem();
		exit(0);
	}
	return 0;
}

void randomize_eta()
//pick initial random values
{
	int i,j;
	for(i=0; i<M; i++) {
		for(j=0; j<clause[i].lits;j++) {
			clause[i].literal[j].eta=randreal();
		}
	}
}

void initmem()
//allocate mem (can be called more than once)
{
	free(perm);
	free(list);
	free(magn);
	free(prod);
	perm=calloc(M,sizeof(int));
	list=calloc(N+1,sizeof(int));
	magn=calloc(N+1,sizeof(double));
	prod=calloc(maxlit,sizeof(double));
#ifdef QUEUE
	queue_init(M);
#endif //QUEUE
	if(!perm||!list||!magn ||!div||!prod){
		fprintf(stderr, "Not enough memory for internal structures\n");
		exit(-1);
	}
}

void compute_pi()
//compute pi products of all vars from scratch
{
	int i,var;
	struct clauselist *cl;
	struct literalstruct *l;
	for(var=1; var<=N; ++var) if(v[var].spin==0) {
		v[var].pi.p=1;
		v[var].pi.m=1;
#ifndef FAST_ITERATION
		v[var].pi.pzero=0;
		v[var].pi.mzero=0;
#endif
		for(i=0,cl=v[var].clauselist;i<v[var].clauses; i++,cl++) if(cl->clause->type) {
			l=cl->clause->literal+cl->lit;
			if(l->bar) {
#ifdef FAST_ITERATION
				v[var].pi.p*=1-l->eta;
#else