サルでも解ける最適化問題

最適化問題用のC++のクラスを作りました。

離散的な状態に、整数のスコアが対応していて、そのスコアを最小化したいという問題。

#include "wl.h"

//スコアの関数を定義
int eval(状態);

int main(int argc, char **argv)
{

    状態＝ランダムな初期値
    int val = eval(状態);
    int vmin = val; // current best

    int s_min = スコアの下限、低めにざっくり見積り
    int s_max = val; //スコアの上限をランダムな状態のスコアに設定
    // クラス初期化
    WangLandau wl(s_min, s_max);

    while(1)
    {
        状態' = 状態とちょっと違う状態
        int vnew = eval(状態');

	if (wl.Accept(val, vnew))
	{
	    // update
	    val = vnew;
	    状態 = 状態';

	    // best score?
	    if (val < vmin)
	    {
		vmin = val;
		printf("#New record %d\n", val);
	    }
	}//accept
    }//loop
}

以下クラスのソース

wl.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();
};

wl.cc

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "wl.h"

// 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;
}