sp.c

人工智能的一个重要方法

				if(1-l->eta>EPS) {
					v[var].pi.p*=1-l->eta;
				} else {
					v[var].pi.pzero++;
				}
#endif
			} else {
#ifdef FAST_ITERATION
				v[var].pi.m*=1-l->eta;
#else
				if(1-l->eta>EPS) {
					v[var].pi.m*=1-l->eta;
				} else {
					v[var].pi.mzero++;
				}
#endif				
			}
		}
	}
}

double update_eta(int cl)
// updates all eta's and pi's in clause cl 
{
	int i;
	struct clausestruct *c;
	struct literalstruct *l;
	struct pistruct *pi;
	double eps;
	double neweta;
	double allprod=1;
	double wt,wn;
	int zeroes=0;
#ifndef FAST_ITERATION
	double m,p;
#endif
	c=&(clause[cl]);
	for(i=0,l=c->literal; i<c->lits; i++,l++) if(v[l->var].spin==0) {
		pi=&(v[l->var].pi);
#ifdef FAST_ITERATION
		if(l->bar) {
			wt=pi->m;
			wn=pi->p/(1-l->eta)*(1-wt*norho);
		} else {
			wt=pi->p;
			wn=pi->m/(1-l->eta)*(1-wt*norho);
		}
#else
		if(l->bar) {
			m=pi->mzero?0:pi->m;
			if(pi->pzero==0) {
				p=pi->p/(1-l->eta);
			} else if (pi->pzero==1 && 1-l->eta<EPS) {
				p=pi->p;
			} else {
				p=0;
			}
			wn=p*(1-m*norho);
			wt=m;
		} else { 
			p=pi->pzero?0:pi->p;
			if(pi->mzero==0) {
				m=pi->m/(1-l->eta);
			} else if (pi->mzero==1 && 1-l->eta<EPS) {
				m=pi->m;
			} else {
				m=0;
			}
			wn=m*(1-p*norho);
			wt=p;
		}
#endif
		prod[i]=wn/(wn+wt);
		if(prod[i]<EPS) {
			if(++zeroes == 2)
				break;
		} else {
			allprod*=prod[i];
		}
	}
	eps=0;
	for(i=0,l=c->literal; i<c->lits; i++,l++) if(v[l->var].spin==0) {
		if(!zeroes){
			neweta=allprod/prod[i];
		} else if (zeroes==1 && prod[i]<EPS) {
			neweta=allprod;
		} else {
			neweta=0;
		}

		pi=&(v[l->var].pi);
#ifdef FAST_ITERATION
		if(l->bar) {
			pi->p*=(1-neweta)/(1-l->eta);
		} else {
			pi->m*=(1-neweta)/(1-l->eta);
		}
#else			
		if(l->bar) {
			if(1-l->eta>EPS) {
				if(1-neweta>EPS) {
					pi->p*=(1-neweta)/(1-l->eta);
				} else {
					pi->p/=(1-l->eta);
					pi->pzero++;
				} 
			} else {
				if(1-neweta>EPS) {
					pi->p*=(1-neweta);
					pi->pzero--;
				} 
			}
		} else {
			if(1-l->eta>EPS) {
				if(1-neweta>EPS) {
					pi->m*=(1-neweta)/(1-l->eta);
				} else {
					pi->m/=(1-l->eta);
					pi->mzero++;
				} 
			} else {
				if(1-neweta>EPS) {
					pi->m*=(1-neweta);
					pi->mzero--;
				} 
			}
		}
#endif
		if(eps<fabs(l->eta-neweta))
			eps=fabs(fabs(l->eta-neweta));
		l->eta=neweta;
	}
	return eps;
}

struct weightstruct compute_field(int var)
//compute H field of the variable var
{
	struct weightstruct H;
	double p,m;
#ifdef FAST_ITERATION
	p=v[var].pi.p;
	m=v[var].pi.m;
#else
	p=v[var].pi.pzero?0:v[var].pi.p;
	m=v[var].pi.mzero?0:v[var].pi.m;
#endif
	H.z=p*m;
	H.p=m-H.z;
	H.m=p-H.z;
	return normalize(H);
}

void print_fields()
//print all H (non-cavity) fields
{
	int var;
	struct weightstruct H;
	compute_pi();
	for(var=1; var<=N; var++) if(v[var].spin==0) {
		H=compute_field(var);
		printf("#H(%i)={%f,%f,%f}\n",var,H.p,H.z,H.m);
	}
}

void print_eta()
//print all etas
{
	int c,l;
	for(c=0; c<M; c++)  if(clause[c].type){
		for(l=0;l<clause[c].lits;l++) if(!v[clause[c].literal[l].var].spin){
			printf("#eta(%i,%i)=%f\n",c,l,clause[c].literal[l].eta);
		}
	}
}

int fix_balanced()
//fix some balanced spin
{
	int var,spin;
	int pap[MAXPAP];
	int npap,minvar=1,totsites=0;
	struct weightstruct H;
	double hphmmin=1000,maxmag=0,summag=0;
	npap=0;
	for(var=1; var<=N && npap<MAXPAP; var++) if(v[var].spin==0) {
		totsites++;
		H=compute_field(var);
		if(fabs(H.p-H.m)<MEBIL && H.z<MEZERO && fabs(H.p-H.m)>MENOBIL) {
			pap[npap++]=var;
			if(fabs(H.p-H.m)<hphmmin) {
				hphmmin=fabs(H.p-H.m);
				minvar=var;
			}
		}
		maxmag=H.p>H.m?H.p:H.m;
		summag+=maxmag;
	}
	if(npap==0)
		return 0;
	var=pap[randint(npap)];
	H=compute_field(var);
	spin=H.p<H.m?1:-1;
	printf("fixed balanced\n");
	return !fix(var,spin);
}

int fix_chunk(int quant)
//fix quant spins at a time, taken from list[]
{
	int i,totsites;
	struct weightstruct H;
	totsites=0;
	i=0;
	while (freespin && quant--) {
		while(v[list[i]].spin)
			i++;
		H=compute_field(list[i]);
		if(verbose>2)
			print_stats();
		if(verbose>1)
			printf("H(%i)={%f,%f,%f} ---> %s\n",list[i],H.p,H.z,H.m,H.p>H.m?"-":"+");
		if(fix(list[i],H.p>H.m?-1:1))
			exit(0);
	}
	return quant;
}

double iterate()
//update etas of clauses in a random permutation order
{
	int cl,vi=0,quant,i;

	double eps,maxeps;
	eps=0;
	maxeps=0;
	for(quant=M-ncl[0];quant;--quant) {
		cl=perm[i=randint(quant)];
		perm[i]=perm[quant-1];
		perm[quant-1]=cl;
		eps=update_eta(cl);
		if(eps>epsilon) {
			vi++;
		}
		if(eps>maxeps) {
			maxeps=eps;
		}
	}
	return maxeps;
}


double compute_sigma()
//compute complexity
{
	int var,cl,i,n;
	double allprod,allprod1,wt,wn;
	double sigmabonds=0, sigmasites=0, sigma=0;
	struct literalstruct *l;
	struct clausestruct *c;
	struct pistruct *pi;
	struct clauselist *cli;
	sigmabonds=0;
	for(cl=0, c=clause; cl<M; c++, cl++) if(c->type) {
		allprod=1;
		allprod1=1;
		for(i=0,l=c->literal; i<c->lits; i++,l++) if(v[l->var].spin==0) {
			pi=&(v[l->var].pi);
#ifdef FAST_ITERATION
			if(l->bar) {
				wn=pi->p/(1-l->eta);
				wt=pi->m;
			} else {
				wn=pi->m/(1-l->eta);
				wt=pi->p;
			}
#else
			if(l->bar) {
				if(1-l->eta>EPS) {
					wn=pi->p/(1-l->eta);
				} else if (pi->pzero==1) {
					wn=pi->p;
				} else {
					wn=0;
				}
				wt=pi->mzero?0:pi->m;
			} else { 
				if(1-l->eta>EPS) {
					wn=pi->m/(1-l->eta);
				} else if (pi->mzero==1) {
					wn=pi->m;
				} else {
					wn=0;
				}
				wt=pi->pzero?0:pi->p;
			}
#endif
			allprod*=wn*(1-wt);
			allprod1*=wt+wn-wt*wn;
		}
		sigmabonds+=log(allprod1-allprod);
	}
	sigmasites=0;
	for(var=1;var<=N;var++) if(v[var].spin==0){
		n=0;
		for(cli=v[var].clauselist,cl=0; cl<v[var].clauses; cli++,cl++) {
			if(cli->clause->type)
				n++;
		}
#ifdef FAST_ITERATION
		wt=v[var].pi.p;
		wn=v[var].pi.m;
#else
		wt=v[var].pi.pzero?0:v[var].pi.p;
		wn=v[var].pi.mzero?0:v[var].pi.m;
#endif
		sigmasites+=log(wt+wn-wt*wn)*(n-1);
	}
	sigma=(sigmabonds-sigmasites)/freespin;
	printf("bonds=%f sites=%f sigma=%f\n", sigmabonds,sigmasites,sigma);
	return sigma;
}

int sequential_converge()
//call iterate() until the maximum difference of eta with previous 
//step is small.
{
	double eps=0;
	int iter=0,cl,quant;
	compute_pi();
	for(cl=0,quant=0; cl<M; cl++) if(clause[cl].type) {
		perm[quant++]=cl;
	}
	do {
		eps=iterate();
		fprintf(stderr,".");
		fflush(stderr);
	} while (eps>epsilon && iter++<iterations);
	if(eps<=epsilon) {
		fprintf(stderr,":-)\n");
		return 1;
	} else {
		fprintf(stderr,"[%f]:-(\n",eps);
		writeformula(fopen(NOCONVERGENCE,"w+"));
		return 0;
	}
	
}

#ifdef QUEUE

extern int Qlength;


//only neighbors to significatively changed messages are updated
int lazy_converge()
{
	int cl,j,k,steps=0,tot=M,aux,cl2,max=0;
	compute_pi();
//	queue_reset();
//	for(cl=0; cl<M; cl++) if(clause[cl].type) {
//		queue_add(cl);
//		tot++;
//	}
	while((cl=queue_get())!=-1) if(clause[cl].type) {
		if(max<Qlength)
			max=Qlength;
		if(steps++%tot==0) {
			fprintf(stderr, ".");
			fflush(stderr);
		}
		if(steps/tot>iterations) {
			fprintf(stderr, ")^:\n");
			return 0;
		}
		
		if(update_eta(cl)>epsilon) {
			for(j=0; j<clause[cl].lits; j++) { 
				aux=clause[cl].literal[j].var;
				if(!v[aux].spin) {
					for(k=0; k<v[aux].clauses; k++) if(clause[cl2=v[aux].clauselist[k].clause-clause].type) {
						queue_add(cl2);
						
					}
				}
			}
		}
	}
	fprintf(stderr,"(^:\n");
	return 1;
}

#endif //QUEUE

int converge()
{
#ifdef QUEUE
	if(uselazy)
		return lazy_converge();
	else
#endif	//QUEUE	
		return sequential_converge();
}