http://q.hatena.ne.jp/1317373518
#include <stdio.h> #include <stdlib.h> #include <strings.h> #include "wl.h" // default group size #define GS 4 int np;// number of attendants int ngroup; // number of groups int ns; // number of stages int *groupid, *nmember, *g_ofs; int *nmeet, *n_smallg, min_nmeet; //np int **seating; //ns*np int **mcount; // np*np void alloc(int np0, int ns0) { int i; np=np0; ns=ns0; int nd=np % GS; ngroup=np/GS; if(nd!=0)ngroup++; #define MYALLOC(x) (int *)malloc(sizeof(int)*(x)) #define MYALLOC2(x) (int **)malloc(sizeof(int *)*(x)) groupid = MYALLOC(np); nmember = MYALLOC(ngroup); g_ofs = MYALLOC(ngroup); nmeet = MYALLOC(np); n_smallg = MYALLOC(np); seating = MYALLOC2(ns); mcount = MYALLOC2(np); for(i=0;i<ns;i++) seating[i] = MYALLOC(np); for(i=0;i<np;i++) mcount[i] = MYALLOC(np); int n_small = GS*ngroup-np; int tot=0; for(i=0;i<ngroup;i++) { if(i<n_small)nmember[i]=GS-1; else nmember[i]=GS; int j; for(j=0;j<nmember[i];j++) groupid[tot+j]=i; g_ofs[i]=tot; tot += nmember[i]; } #if 0 printf("#"); for(i=0;i<np;i++) { printf("%d", groupid[i]); } printf("\n"); #endif } int init_mcount(int **st, int **t) { int i1,i2,s,g; for(i1=0;i1<np;i1++) for(i2=0;i2<np;i2++) t[i1][i2]=0; for(s=0;s<ns;s++) { for(g=0;g<ngroup;g++) { for(i1=0;i1<nmember[g];i1++){ for(i2=i1+1;i2<nmember[g];i2++){ int p1=st[s][g_ofs[g]+i1]; int p2=st[s][g_ofs[g]+i2]; t[p1][p2]++; t[p2][p1]++; } } } } int res=0; for(i1=0;i1<np;i1++){ for(i2=i1+1;i2<np;i2++){ int n=t[i1][i2]; res += (n*(n-1))/2; }} for(i1=0;i1<np;i1++)nmeet[i1]=0; for(i1=0;i1<np;i1++) for(i2=0;i2<np;i2++)if(t[i1][i2]+t[i2][i1]!=0)nmeet[i1]++; return res; } int p1o[GS], p2o[GS]; int eval_exchange(int s, int g1, int g2, int i1, int i2) { int o1=g_ofs[g1]; int o2=g_ofs[g2]; int p1 = seating[s][o1+i1]; int p2 = seating[s][o2+i2]; int k; for(k=1;k<nmember[g1];k++) p1o[k-1]=seating[s][o1+(i1+k)%nmember[g1]]; for(k=1;k<nmember[g2];k++) p2o[k-1]=seating[s][o2+(i2+k)%nmember[g2]]; int de=0; int n; //remove p1 for(k=0;k<nmember[g1]-1;k++) de -= mcount[p1][p1o[k]] -1; //remove p2 for(k=0;k<nmember[g2]-1;k++) de -= mcount[p2][p2o[k]] -1; //add p1 for(k=0;k<nmember[g2]-1;k++) de += mcount[p1][p2o[k]]; //add p2 for(k=0;k<nmember[g1]-1;k++) de += mcount[p2][p1o[k]]; return de; } void accept_exchange(int s, int g1, int g2, int i1, int i2) { int o1=g_ofs[g1]; int o2=g_ofs[g2]; int p1 = seating[s][o1+i1]; int p2 = seating[s][o2+i2]; int k; //remove p1 for(k=0;k<nmember[g1]-1;k++) { mcount[p1][p1o[k]]--; mcount[p1o[k]][p1]--; if(mcount[p1][p1o[k]]==0) { nmeet[p1]--; nmeet[p1o[k]]--; } } //remove p2 for(k=0;k<nmember[g2]-1;k++) { mcount[p2][p2o[k]]--; mcount[p2o[k]][p2]--; if(mcount[p2][p2o[k]]==0) { nmeet[p2]--; nmeet[p2o[k]]--; } } #if 0 int k1,k2; for(k1=0;k1<np;k1++) for(k2=0;k2<np;k2++) if(mcount[k1][k2]<0) { printf("#err G %d %d I %d %d\n",g1,g2,i1,i2); exit(1); } #endif //add p1 for(k=0;k<nmember[g1]-1;k++) { mcount[p2][p1o[k]]++; mcount[p1o[k]][p2]++; if(mcount[p2][p1o[k]]==1) { nmeet[p2]++; nmeet[p1o[k]]++; } } //add p2 for(k=0;k<nmember[g2]-1;k++) { mcount[p1][p2o[k]]++; mcount[p2o[k]][p1]++; if(mcount[p1][p2o[k]]==1) { nmeet[p1]++; nmeet[p2o[k]]++; } } seating[s][o1+i1]=p2; seating[s][o2+i2]=p1; } // Debug void dump_mcount() { int i1,i2; for(i1=0;i1<np;i1++){ for(i2=0;i2<np;i2++){ printf("%02d ", mcount[i1][i2]); } printf("\n"); } printf("\n"); } void get_shuffle(int *g1, int *g2, int *i1, int *i2) { *g1=(rand()>>8)%ngroup; *g2=(*g1+1+((rand()>>8)%(ngroup-1)))%ngroup; *i1=(rand()>>8)%(nmember[*g1]); *i2=(rand()>>8)%(nmember[*g2]); } int get_nmin() { int ii1; int nmin0=999; for(ii1=0;ii1<np;ii1++) { if(nmeet[ii1]<nmin0)nmin0=nmeet[ii1]; } return nmin0; } int main(int argc, char **argv) { int s,ss; int i1,i2,g1,g2; int mmax=999; if(argc!=3) { printf("%s n.of.stages n.of.people\n",argv[0]); exit(1); } int ns0 = atoi(argv[1]); int np0 = atoi(argv[2]); alloc(np0, ns0); //init for(s=0;s<ns;s++) for(i1=0;i1<np;i1++)seating[s][i1]=i1; int enow=init_mcount(seating, mcount); /////////////// int emax=enow; printf("#Emax=%d\n",emax); // Init energy range WangLandau wl(0, emax/2); // Find Local min for(ss=0;ss<10000;ss++) { for(s=1;s<ns;s++) { get_shuffle(&g1, &g2, &i1,&i2); int de=eval_excchange(s, g1,g2,i1,i2); if(de<0) { accept_exchange(s,g1,g2,i1,i2); enow+=de; } } } int ebest=enow; printf("# Enow=%d\n",enow); int nmin=0; ss=0; while(1) { for(s=1;s<ns;s++) { get_shuffle(&g1, &g2, &i1,&i2); int de=eval_exchange(s, g1,g2,i1,i2); // Wang-Landau update if (wl.Accept(enow, enow+de)) { accept_exchange(s,g1,g2,i1,i2); enow+=de; int nmin0=get_nmin(); //printf("#nmin %d %d\n", nmin0, nmin); int newbest=0; if(enow<ebest) { ebest=enow; newbest=1; //nmin=0; } else if(enow==ebest) newbest=2; int ii1,ii2; int mmax0=0; for(ii1=0;ii1<np;ii1++) for(ii2=ii1+1;ii2<np;ii2++) if(mmax0<mcount[ii1][ii2])mmax0=mcount[ii1][ii2]; // New record if((nmin0>nmin &&newbest!=0)||(nmin0==nmin && newbest==1)||(nmin0==nmin && newbest==2 && mmax0<mmax)) { int i; nmin=nmin0; mmax=mmax0; printf("# ebest %d Nmin %d Mmax %d\n", enow, nmin, mmax); char fn[256]; sprintf(fn, "best-%02d-%02d.log", ns, np); FILE *fp=fopen(fn,"w"); fprintf(fp, "# %d\n", np); // dump for(i=0;i<np;i++) { fprintf(fp, "%02d %d ",i+1, nmeet[i]); int s; for(s=0;s<ns;s++) { int i2; int g0=-1; for(i2=0;i2<np;i2++)if(seating[s][i2]==i)g0=groupid[i2]; fprintf(fp, "%c",'A'+g0); } for(i2=0;i2<np;i2++) { if(mcount[i][i2]>1) fprintf(fp, " %02d", i2+1); if(mcount[i][i2]>2) fprintf(fp,"x%d", mcount[i][i2]); } fprintf(fp, "\n"); } fclose(fp); //if(nmin==15)exit(1); }//nmin }//accep }//stage ss++; if(ss% 10000000==0) wl.DumpHist("h"); //if(ss==10000000) exit(1); }//step }
