Solve http://q.hatena.ne.jp/1355418360
#include <stdio.h> #include <math.h> #include <stdlib.h> class WangLandau { public: // 1-dim constructor WangLandau(long long emin, long long emax,double pip = 0.2, int histmin=10000); // constructor for multi-dimensional energy optimization WangLandau(int dim, long long *emin, long long *emax,double pip = 0.2, int histmin=10000); ~WangLandau(); // 1-dim version bool Accept(long long enow, long long enew); // n-dim version bool Accept(long long *enow, long long *enew); // save/load histogram void DumpHist(char *fn); void LoadHist(char *fn); void SetQuench(bool q0); private: //Energy range long long EMin, EMax; // multi-dimensional histogram int Dim; #define WLMAXDIM 3 long long EMind[WLMAXDIM], EMaxd[WLMAXDIM]; int ELen[WLMAXDIM]; // LogW increase per hit double Pip; // Tabulise exp(-n*Pip) int *ThTab; void InitThTab(); // Max table size: if n exceeds dEmax, exp(-n*Pip) is approximated by 0 int dEmax; // Switch to quench mode bool Quench; // Energy histogram unsigned int *Hist; unsigned int *HistOld; int HistSize; void DumpHist1D(char *fn); void DumpHistnD(char *fn); void LoadHist1D(char *fn); void LoadHistnD(char *fn); // refresh histogram when enough hits are recorded int HistMin; int HistAdr(long long *ene); void ClrHist(); void RefreshHist(); // EMin boundary breached! void ExtendHist1D(); }; // 1D constructor WangLandau::WangLandau(long long emin, long long emax,double pip, int histmin) { Dim=1; EMin = emin; EMax = emax; HistSize = emax-emin+1; Hist = (unsigned int *)malloc(HistSize*sizeof(int)); HistOld = (unsigned int *)malloc(HistSize*sizeof(int)); if (Hist==NULL || HistOld==NULL) { fprintf(stderr, "#WangLandau:: Can't alloc Hist, emax-emin=%d too large\n", HistSize); exit(1); } ClrHist(); Pip = pip; HistMin = histmin; Quench = false; dEmax = (int)(5.0/pip); ThTab = (int *)malloc(dEmax * sizeof(int)); if (ThTab==NULL) { fprintf(stderr, "#WangLandau:: Can't alloc ThTab, size=%d pip %f too small\n", dEmax, pip); exit(1); } InitThTab(); }; WangLandau::~WangLandau() { free(Hist); free(HistOld); free(ThTab); }; // n-Dim constructor WangLandau::WangLandau(int dim, long long *emin, long long *emax, double pip , int histmin) { if (dim > WLMAXDIM) { printf("#WangLandau: dim=%d exceeds max limit and not practical\n", dim); exit(1); } Dim=dim; int i; HistSize = 1; for(i=0;i<dim;i++) { EMind[i] = emin[i]; EMaxd[i] = emax[i]; ELen[i] = emax[i]-emin[i]+1; HistSize *= ELen[i]; } Hist = (unsigned int *)malloc(HistSize*sizeof(int)); HistOld = (unsigned int *)malloc(HistSize*sizeof(int)); if (Hist==NULL || HistOld==NULL) { fprintf(stderr, "#WangLandau:: Can't alloc Hist, size=%d too large\n", HistSize); exit(1); } ClrHist(); Pip = pip; HistMin = histmin; Quench = false; dEmax = (int)(5.0/pip); ThTab = (int *)malloc(dEmax * sizeof(int)); if (ThTab==NULL) { fprintf(stderr, "#WangLandau:: Can't alloc ThTab, size=%d pip %f too small\n", dEmax, pip); exit(1); } InitThTab(); }; void WangLandau::SetQuench(bool q0) { Quench = q0; } void WangLandau::ClrHist() { int i; for(i=0;i<HistSize;i++) Hist[i]=HistOld[i]=0; } void WangLandau::InitThTab() { int i; ThTab[0] = RAND_MAX/2; for(i=1;i<dEmax;i++) ThTab[i] = (int)(RAND_MAX * exp(-Pip*i)); } void WangLandau::DumpHist(char *fn) { if(Dim==1) DumpHist1D(fn); else DumpHistnD(fn); } void WangLandau::DumpHist1D(char *fn) { FILE *fp=fopen(fn, "w"); int i; int di = 10;// energy interval to calc temperature for(i=0;i<HistSize;i++) { int h = Hist[i]; if(h==0)continue; double beta=0.0; if (i+10 <HistSize) beta = (Hist[i+di] - Hist[i])/(1.0*di); fprintf(fp, "%lld %d %.20e %e\n", EMin+i, h, -h*Pip, beta); } fclose(fp); } void WangLandau::DumpHistnD(char *fn) { FILE *fp=fopen(fn, "w"); char fn2[1024]; sprintf(fn2,"%s-c",fn); FILE *fp2=fopen(fn2, "w"); sprintf(fn2,"%s-n",fn); FILE *fp3=fopen(fn2, "w"); int i, d; int div[WLMAXDIM]; for(d=0;d<Dim;d++) div[d] = 1+ELen[d]/256; int laste1 = -1; for(i=0;i<HistSize;i++) { long long ee[WLMAXDIM]; int h = Hist[i]; int ho = HistOld[i]; int ii = i; bool do_dump2 = true; for(d=0;d<Dim;d++) { int e2 = ii%ELen[d]; long long ene = EMind[d] + e2; ee[d]= EMaxd[d]-e2; if ((e2 % div[d]) !=0) do_dump2=false; ii /= ELen[d]; if(h!=0)fprintf(fp, "%lld ", ene); if(h!=0 && h!= ho)fprintf(fp3, "%lld ", ene); } if (do_dump2) { if(ee[1]!=laste1) { fprintf(fp2, "\n"); laste1=ee[1]; } int hh=Hist[HistAdr(ee)]; //if (hh!=0) { for(d=0;d<Dim;d++) fprintf(fp2, "%lld ", ee[d]); // temperature ee[0]++; int dh1 = 0; int ad1=HistAdr(ee); if (ad1!=-1) dh1 = Hist[ad1]-hh; int dh2 = 0; ee[0]--; ee[1]++; ad1=HistAdr(ee); if (ad1!=-1) dh2 = Hist[ad1]-hh; fprintf(fp2, "%d %e %e", hh, dh1*Pip, dh2*Pip); } } if(h!=0) fprintf(fp, "%d %.20e\n", h, -h*Pip); if(h!=0 && h!= ho)fprintf(fp3, " %d\n", h-ho); HistOld[i]=Hist[i]; } fclose(fp); fclose(fp2); fclose(fp3); } void WangLandau::LoadHist(char *fn) { if(Dim==1) LoadHist1D(fn); else LoadHistnD(fn); } void WangLandau::LoadHist1D(char *fn) { ClrHist(); FILE *fp=fopen(fn, "r"); char buf[256]; while(1) { fgets(buf,255,fp); if (feof(fp))break; long long i; int h; double dd; sscanf(buf, "%lld %d %le", &i, &h, &dd); if (i>=EMin && i<=EMax) { //Hist[i-EMin]= (int)(-dd/Pip); Hist[i-EMin]= h; } } fclose(fp); int i; for(i=1;i<EMax-EMin+1;i++) { unsigned int h = Hist[i]; if (h<Hist[i-1]) Hist[i]=Hist[i-1]; } } int WangLandau::HistAdr(long long *ene) { int d; int adrs=0, inc=1; for(d=0;d<Dim;d++) { int ee = ene[d] - EMind[d]; if (ee<0 || ee >=ELen[d]) return -1; adrs += inc * ee; inc *= ELen[d]; } return adrs; } void WangLandau::LoadHistnD(char *fn) { ClrHist(); FILE *fp=fopen(fn, "r"); char buf[256]; int nh=0, nh2=0; while(1) { fgets(buf,255,fp); if (feof(fp))break; long long ene[WLMAXDIM]; int h; double dd; //Hack! if (Dim==2) sscanf(buf, "%lld %lld %d %le", &ene[0], &ene[1], &h, &dd); if (Dim==3) sscanf(buf, "%lld %lld %lld %d %le", &ene[0], &ene[1], &ene[2], &h, &dd); int adrs= HistAdr(ene); if (adrs== -1)continue; if (adrs>=HistSize) { printf("#Its a bug!\n");exit(1); } Hist[adrs]= h; nh++; nh2+=h; } printf("#WangLandau: Hist %s Loaded %d entries %d hits\n",fn, nh, nh2); // Fill In fclose(fp); } // 1D bool WangLandau::Accept(long long enow, long long enew) { if (Quench) { if (enew <= enow) return true; if (enow == enew && rand()<RAND_MAX/2) return true; return false; } if (enow > EMax || enew > EMax) { if (enew < enow) { if (enew <=EMax) Hist[enew-EMin]++; return true; } else { if (enow <=EMax) Hist[enow-EMin]++; return false; } } if (enew<EMin) ExtendHist1D(); int en1 = enow-EMin; int en2 = enew-EMin; int de = Hist[en2] - Hist[en1]; if (de<0) { Hist[en2]++; return true; } else if (de >= dEmax) { Hist[en1]++; return false; } int th = ThTab[de]; if (rand() < th) { Hist[en2]++; return true; } Hist[en1]++; return false; }; // n-D bool WangLandau::Accept(long long *enow, long long *enew) { if (Quench) { // use ene[0] if (enew[0] <= enow[0]) return true; if (enow[0] == enew[0] && rand()<RAND_MAX/2) return true; return false; } int adNow = HistAdr(enow); int adNew = HistAdr(enew); //if(enew[1]!=enow[1]) printf("dE=%lld %lld %c %c\n", enew[0]-enow[0], enew[1]-enow[1], (adNow==-1)?'O':'I', (adNew==-1)?'O':'I'); if (adNew == -1 || adNow == -1) { if (adNow != -1) { Hist[adNow]++; return false; } int es1=0, es2=0; int d; for(d=0;d<Dim;d++) { es1+=enow[d]; es2+=enew[d]; if(enow[d]<enew[d]) return false; } return true; //if (es2 <= es1) return true; //return false; } int de = Hist[adNew] - Hist[adNow]; if (de<0) { Hist[adNew]++; return true; } else if (de >= dEmax) { Hist[adNow]++; return false; } int th = ThTab[de]; if (rand() < th) { Hist[adNew]++; return true; } Hist[adNow]++; return false; }; void WangLandau::ExtendHist1D() { long long EMinOld = EMin; int elen = (EMax - EMin)*3/2; EMin = EMax +1 - elen; printf("#WangLandau: EMin reached, lowering EMin %lld to %lld and extending buffer %f Mb\n", EMinOld, EMin,elen*4.0/1e6); unsigned int *newHist = (unsigned int *)malloc(elen*sizeof(int)); int i; for(i=0;i<elen;i++) newHist[i]=0; // migrate data for(i=EMinOld;i<=EMax;i++) newHist[i-EMin] = Hist[i-EMinOld]; free(Hist); Hist = newHist; } // default group size #define GS 5 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 int eofn(int n) { return n*(n-1)/2; } 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]; // energy //res += (n*(n-1))/2; res += eofn(n); }} 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++) { int n= mcount[p1][p1o[k]]; de += eofn(n-1)-eofn(n); //de -= mcount[p1][p1o[k]] -1; } //remove p2 for(k=0;k<nmember[g2]-1;k++) { //de -= mcount[p2][p2o[k]] -1; int n= mcount[p2][p2o[k]]; de += eofn(n-1)-eofn(n); } //add p1 for(k=0;k<nmember[g2]-1;k++) { int n=mcount[p1][p2o[k]]; //de += mcount[p1][p2o[k]]; de += eofn(n+1)-eofn(n); } //add p2 for(k=0;k<nmember[g1]-1;k++) { int n=mcount[p2][p1o[k]]; de += eofn(n+1)-eofn(n); //de += mcount[p2][p1o[k]]; } return de; } int eval_personal(int m1) { int res=0; int i; for(i=0;i<np;i++)if(i!=m1)res+=eofn(mcount[m1][i]); return res; } int max_personal() { int max=eval_personal(0); int min=max; int sum=max*max; int i; for(i=1;i<np;i++) { int d=eval_personal(i); if(d>max)max=d; if(d<min)min=d; sum+=d*d; } //return max-min; return sum; } 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_exchange(s, g1,g2,i1,i2); if(de<0) { accept_exchange(s,g1,g2,i1,i2); enow+=de; } } } int ebest=enow; int e2best = max_personal(); 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; int maxp; if(enow==ebest) maxp=max_personal(); if(enow<ebest || (enow==ebest && maxp<e2best)) { ebest=enow; e2best = maxp; 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 Pmax %d\n", enow, nmin, mmax, e2best); printf( "# %d = %d", np, nmember[0]); for(i=1;i<ngroup;i++) { printf( "+%d", nmember[i]); } printf( "\n"); // dump for(i=0;i<np;i++) { printf( "%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]; printf( "%c",'A'+g0); } int mm; for(mm=mmax;mm>=2;mm--) { for(i2=0;i2<np;i2++) { if(mcount[i][i2]!=mm)continue; if(mcount[i][i2]>1) printf( " %02d", i2+1); if(mcount[i][i2]>1) printf("x%d", mcount[i][i2]); } } printf( "\n"); } //if(nmin==15)exit(1); }//nmin }//accep }//stage ss++; }//step }