SuperCon 2006

#include    <stdio.h>
#include    <stdlib.h>
#include    <math.h>
#include    <string.h>

#define N 3000
typedef unsigned char mtype;

mtype mij0[N][N];
// for compressed row matrix storage
int nnz[N];
unsigned short *jnz[N];
mtype *mij[N];

// Number of switch
#define MAXQ 6
unsigned char spin[N];
unsigned char spin_bak[N];
int hist[MAXQ];

int e_total, e_best, e_bound;

// For exp(-beta e) table
#define EMAX 1024
int etab[EMAX];
// Probability table precision
#define PMAX (1<<16)

void init_etab(double beta)
{
    int i;
    for(i=0;i<EMAX;i++)
	etab[i]=(int)(PMAX*exp(-beta*i));
}

// Cal. effective field
void cal_h(int i, int h[MAXQ])
{
    int j,q;

    unsigned short *jptr = jnz[i];
    mtype *mptr = mij[i];

    for(q=0;q<MAXQ;q++) h[q]= hist[q];
    h[spin[i]]--;

#if 0
    for(j=0;j<N;j++)h[spin[j]] -= mij0[i][j];
#else
    // Hot spot
    while(*mptr != 0)
    {
	h[ spin[*jptr] ] -= *mptr;
	jptr++;
	mptr++;
    }
#endif
}

// Heatbath update
void do_hb(int i, int zero)
{
    int j, q, q0, si_new;
    int h[MAXQ], hmin;
    int p[MAXQ], psum;

    cal_h(i, h);

    hmin=h[0];
    q0=0;
    for(q=1;q<MAXQ;q++)
    {
	if(hmin>h[q])
	{
	    hmin=h[q];
	    q0=q;
	}
    }

    if(zero==1)
    {
	// Quench
	si_new=q0;
    }
    else
    {
        psum=0;
        for(q=0;q<MAXQ;q++)
        {
	    int e = h[q]-hmin;
	    p[q] = (e>=EMAX) ? 0:etab[e];
	    psum += p[q];
        }
        psum = (int)(psum*1.0*rand()/RAND_MAX);

        q=0;
        while(psum >p[q])
        {
	    psum -= p[q];
	    q++;
        }
        si_new = q;
    }

    hist[spin[i]]--;
    hist[si_new]++;

    e_total += -h[spin[i]]+h[si_new];

    spin[i]=si_new;

    if(e_total < e_best)
    {
	e_best = e_total;
	memcpy(spin_bak,spin,N);
    }
}

void read_file(char *fn)
{
    char buf[N*2+100];
    FILE *fp=fopen(fn,"r");
    int i,j;

    fgets(buf,999,fp);

    e_bound = 0;

    for(i=0;i<N;i++)
    {
	int nz=0;
        fgets(buf,N+10,fp);
	for(j=0;j<N;j++)
	{
	    mij0[i][j] = buf[j]-'0';
	    if(mij0[i][j]!=0 && i!=j) nz++;
	}
	nnz[i]=nz;
	mij0[i][i]=0;
	e_bound += nz;
    }
    fclose(fp);

    e_bound/=2;

    // Prepare compressed row data
    for(i=0;i<N;i++)
    {
	int nz=nnz[i];
	jnz[i] = (unsigned short *)malloc(sizeof(unsigned short)*nz);
	mij[i] = (mtype *)malloc(sizeof(mtype)*(nz+1));

	nz=0;
	for(j=0;j<N;j++)
	{
	    if(mij0[i][j]!=0)
	    {
		jnz[i][nz] = j;
		mij[i][nz] = mij0[i][j];
		nz++;
	    }
	}
	mij[i][nz] =0;
    }
}

int ene()
{
    int c1=0;
    int c2=0;
    int i,j;

    for(i=0;i<N;i++) {
    for(j=0;j<i;j++) {
	c1 += (spin[i]==spin[j])? 0:mij0[i][j];
	c2 += (spin[i]==spin[j])? 1:0;
    }}
    return c1+c2;
}


// Mean field critical beta
double beta_mf(double m_ave)
{
    int i,j;

    return MAXQ*1.0/N/m_ave;
}

void print_m()
{
    int q;
    for(q=0;q<MAXQ;q++) printf(" %d",hist[q]);
}

int main(int argc, char **argv)
{
    int i,j,q,step, maxstep;
    double beta, beta_min, beta_max, beta_inc;
    double beta_c;

    read_file(argv[1]);
    maxstep = atoi(argv[2]);

    beta_c = beta_mf(9.0);
    printf("#MF Beta_c = %f\n", beta_c);

    beta=beta_min=0.04;
    beta_max = 2.0;
    beta_inc = exp(log(beta_max/beta_min)/maxstep);

    for(q=0;q<MAXQ;q++) hist[q]=0;

    for(i=0;i<N;i++)
    {
	spin[i] = (rand()>>8)%MAXQ;
	hist[spin[i]]++;
    }

    e_total = e_best = ene();

    // Main loop
    for(step=0;step<maxstep;step++)
    {
	init_etab(beta);
	for(i=0;i<N;i++)
	{
	    do_hb(i,0);
	}
	if(step%10==0)
	{
	    int k;
	    printf("%e %d  %d\n",beta,e_best, e_best-e_bound);
	}
	if (e_total == e_bound) break;
	beta *= beta_inc;
    }

    // Quench from the best config
    memcpy(spin,spin_bak,N);
    for(q=0;q<MAXQ;q++) hist[q]=0;
    for(i=0;i<N;i++) hist[spin[i]]++;

    for(step=0;step<40;step++)
    {
	for(i=0;i<N;i++) do_hb(i,1);
    }

    printf("#%d steps %d\n",maxstep,ene());

    printf("# CSize ");
    print_m();
    printf("\n");

    //put_image("matrix.ppm", spin, mij0);
}