Seating Optimization

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

}