2013 UPDATE

island record.jpg

 

 

 

Robert Ziff at the University of Michigan  -

Ziff got interested in this idea that a randomly distributed 2 level system could retain more water than a randomly distributed 3 level system when the size of the square was greater than 51 x51.  Percolation theory nicely explains this counter intuitive finding.  There was a resulting Physical Review Letter on this in 2012

1.  a b c Knecht, Craig; Walter Trump, Daniel ben-Avraham, Robert M. Ziff (2012). "Retention capacity of random surfaces". Physical Review Letters 108 (4): 045703. arXiv:1110.6166.Bibcode:2012PhRvL.108d5703Kdoi:10.1103/PhysRevLett.108.045703.

 

Ziff started a Wikipedia page on this subject  ..

http://en.wikipedia.org/wiki/Water_retention_on_mathematical_surfaces

which inspired me to add content to the wiki below ….

http://en.wikipedia.org/wiki/Associative_magic_square

Ziff gave a lecture on this general topic which is provided on the link below.

http://www.compphys.uni-oldenburg.de/en/download/Alexander/dpg_school/ZiffTalkPart5.pdf

#include "stdlib.h"

#include "stdio.h"

#include "math.h"

 

void randinit(long seed);

 

#define DEPTH 2

#define HEIGHT 54   //four larger than the size of the lattice

#define WIDTH HEIGHT

#define RUNSMAX  50000

#define PRINTFREQ 1000

#define S 65535

#define SEED 120

#define DIRMAX 4

#define M  16383             

#define PutOnStack(X,Y) {xlist[pptr]=X; ylist[pptr]=Y; ++pptr;}

#define NewRandomInteger (++nd,ra[nd&M] = ra[(nd-471)&M]^ra[(nd-1586)&M]^ra[(nd-6988)&M]^ra[(nd-9689)&M])                                              

 

long ra[M+1], nd;

long  lat[WIDTH][HEIGHT], lat2[WIDTH][HEIGHT];

long xlist[S+1], ylist[S+1];

FILE    *fp1;

 

int main(void)

{

      long    dx[4] = {1, 0, -1, 0}, dy[4] = {0, 1, 0, -1};

      long x, y, xo, yo, xp, yp, dir, prob, gptr, pptr, nocc, runs, i, j, hit;

      long    nocctot, level;

      double  sum, sumtot;

      FILE    *fp1;

     

      randinit(SEED);

     

 

      for (x = 0; x < WIDTH; ++x) lat2[x][0] = lat2[x][HEIGHT-1] = -1;

      for (y = 0; y < HEIGHT; ++y) lat2[0][y] = lat2[WIDTH-1][y] = -1;

     

      for (x = 1; x < WIDTH-1; ++x) lat[x][0] = lat[x][HEIGHT-2] = 0;

      for (y = 1; y < HEIGHT-1; ++y) lat[0][y] = lat[WIDTH-2][y] = 0;

 

      for (runs = 1; runs <= RUNSMAX; ++runs)

      {

                 

            for (x = 2; x < WIDTH-2; ++x)

            for (y = 2; y < HEIGHT-2; ++y)                 

                  lat[x][y] = 1+(int)(DEPTH*(NewRandomInteger/2147483648.0));

                 

            //faster way for two levels lat[x][y] = 1+(NewRandomInteger & 1);

           

            for (x = 1; x < WIDTH-1; ++x)

            for (y = 1; y < HEIGHT-1; ++y)

                  lat2[x][y] = 0;  

                 

             

            /* to print out the lattice for program checking

            for (y = 0; y < HEIGHT; ++y)

            { 

                  for (x = 0; x < WIDTH; ++x)

            printf("%3ld",lat[x][y]);

                  printf("\n");

            }

            printf("\n");

            for (y = 0; y < HEIGHT; ++y)

            {

                  for (x = 0; x < WIDTH; ++x)

            printf("%3ld",lat2[x][y]);

                  printf("\n");

            }

            printf("\n");    

            */   

 

             

          pptr = 0;    

            level = 1;

            PutOnStack(1,1)        

            lat2[1][1] = level;

                  do

                  {    

                        y = 2147483647;

                        for (i = 0; i < pptr; ++i)

                              if (lat[xlist[i]][ylist[i]] < y) {y = lat[xlist[i]][ylist[i]]; j = i;}

                       

                        x = xlist[j];

                        y = ylist[j];

                        if (lat[x][y] > level) level = lat[x][y];

                        lat2[x][y] = -level;

                        xlist[j] = xlist[pptr-1];

                        ylist[j] = ylist[pptr-1];

                        --pptr;

                 

                        for (dir = 0; dir < DIRMAX; ++dir)

                        {

                              xp = x + dx[dir];

                              yp = y + dy[dir];

                       

                              if (lat2[xp][yp] >= 0)

                              {

                                    lat2[xp][yp] = -1;

                                    PutOnStack(xp,yp)

                              }

                        }

                  }

                  while (pptr > 0);

             

             

            /* to print out the lattice for program checking

            printf("\n");

            for (y = 0; y < HEIGHT; ++y)

            { 

                  for (x = 0; x < WIDTH; ++x)

            printf("%3ld",lat[x][y]);

                  printf("\n");

            }

            printf("\n");

            for (y = 0; y < HEIGHT; ++y)

            {

                  for (x = 0; x < WIDTH; ++x)

            printf("%3ld",-lat2[x][y]);

                  printf("\n");

            }

           

            printf("\n");          

            for (y = 0; y < HEIGHT; ++y)

            {

                  for (x = 0; x < WIDTH; ++x)

            printf("%3ld",-lat[x][y]-lat2[x][y]);

                  printf("\n");

            }

            */

           

           

            sum = 0;

            for (x = 2; x < WIDTH-2; ++x)

            for (y = 2; y < HEIGHT-2; ++y)     

            sum += lat[x][y]+lat2[x][y];

           

            //printf("sum %10ld\n", -sum);

            sumtot -= sum;

           

            if ((runs % PRINTFREQ) == 0)

            printf("runs sumtot %10ld%18.8e%18.8e\n", runs, sumtot, sumtot/runs);

            }

      }

 

void randinit(long seed)

{    

      double a, ee = -1 + 1/2147483648.0;

      long i;

      extern long nd, ra[M+1];

     

      a = seed/2147483648.0;

      for (nd = 0; nd <= M; nd++)

      {

            a *= 16807;

            a += ee * (long)(a);

            if (a >= 1) a += ee;

            ra[nd] = (long) (2147483648.0 * a);

      }

      nd = M;

      for(i = 0; i<100001; i++)

            NewRandomInteger;

}

 

 

Seung Ki Baek

Seung Ki worked independently on this model.  He wrote an additional paper on the subject.

Baek.jpg

He graciously contributed this code for finding the water retention of any size of square.

///////////////////////////////////////////

// Water retention model by Craig Knecht  

//

// three levels

// depth-first burning algorithm           

// tree structure to find avaiable sites

//

// 2011. 10.

// Seung Ki Baek

////////////////////////////////////////////

 

#include <stdlib.h>

#include <stdio.h>

#include <math.h>

 

#define DIM 2

#define KMAX (2*DIM)

 

#define UNEXAMINED 0

#define BLOCK 1

#define WATER 2

#define HIGHBLOCK 3

#define HIGHWATER 4

#define BOUNDARY 5

 

#define MALLOC(p,s) \

if(!((p)=malloc(s))) { \

fprintf(stderr,"Insufficient memory");\

exit(EXIT_FAILURE);\

}

 

int totnum(int r, int D)

{

  int i, n;

  n = 1;

  for(i=0; i<D; i++) n *= r;

  return n;

}

 

void linkage(int k, int l, int* KK, int** nn)

{

  nn[k][KK[k]] = l;  KK[k]++;

  nn[l][KK[l]] = k;  KK[l]++;

}

 

void tessellation(int r, int D, int *KK, int** nn)

{

  int i, j, k, l, kd, kd2;

  int v[D];

  int N = totnum(r, D);

  for(i=0; i<N; i++) KK[i] = 0;

 

  for(k=0; k<N; k++){

    for(i=0; i<D; i++){

      kd = (int)(pow(r,i)+0.5);

      kd2 = (int)(pow(r,i+1)+0.5);

      l = k+kd;

      if(l>=N) continue;

      if((k%kd2)/kd==r-1) continue;

      linkage(k,l,KK,nn);

    }

  }

}

 

void grow(int p, double prob, int* KK, int** nn, int* state, int* size, int*

touch, int* tree, int* level, int height)

{

  int i, j, k;

 

  for(i=0; i<KK[p]; i++){

    j = nn[p][i];

    if(state[j]==BOUNDARY){ (*touch)++; continue; }

    if(state[j]==WATER || state[j]==BLOCK || state[j]==HIGHBLOCK) continue;

    state[j] = WATER;

    (*size)++;

    for(k=0; k<height; k++) tree[level[k] + (j >> (height-1-k))]--;

    grow(j, prob, KK, nn, state, size, touch, tree, level, height);

  }

}

 

void regrow(int p, double prob, int* KK, int** nn, int* state, int* size, int*

touch, int* tree, int* level, int height)

{

  int i, j, k;

 

  for(i=0; i<KK[p]; i++){

    j = nn[p][i];

    if(state[j]==BOUNDARY){ (*touch)++; continue; }

    if(state[j]==HIGHWATER || state[j]==HIGHBLOCK) continue;

    state[j] = HIGHWATER;

    (*size)++;

    for(k=0; k<height; k++) tree[level[k] + (j >> (height-1-k))]--;

    regrow(j, prob, KK, nn, state, size, touch, tree, level, height);

  }

}

 

void draw(int* state, int L, int avail, int size, int remain, int touch, int water)

{

  int i, j, k;

  printf("site %d: size %d\t remain %d\t touch %d\t water %d\n", avail, size, remain, touch, water);

  for(i=0; i<L; i++){

    for(j=0; j<L; j++){

      k = i*L + j;

      if(state[k]==BOUNDARY) printf(" #");

      else if(state[k]==UNEXAMINED) printf("  ");

      else if(state[k]==BLOCK) printf(" .");

      else if(state[k]==HIGHBLOCK) printf(" o");

      else if(state[k]==HIGHWATER) printf(" ~");

      else printf(" +");

    }

    printf("\n");

  }

}

 

main(int argc, char** argv)

{

  int N, N1, Ne;

  int i, j, k, l;

  int *KK, **nn;

  int L;

  double r, prob1, prob2;

  int ensemble;

  int DRAW;

  int *state;

  int remain, avail;

  int size, touch;

  int water;

  double frac;

  double avg, stdev;

 

  int height;

  int length;

  int *level;

  int *tree;

 

  if(argc!=6){ printf("perc L prob1 prob2 ensemble DRAW\n"); exit(1);}

 

  L = atoi(argv[1]);

  prob1 = atof(argv[2]);

  prob2 = atof(argv[3]);

  ensemble = atoi(argv[4]);

  DRAW = atoi(argv[5]);

 

  if(prob1 + prob2>1.0){

    printf("%d %f %f %e %e\n", L, prob1, prob2, 2.0, 0.0);

    return;

  }

 

  N = totnum(L, DIM);

  N1 = N - 1;

  Ne = totnum(L-2,DIM); // effective number of sites excluding boundary points

  MALLOC(KK, N*sizeof(int));

  MALLOC(nn, N*sizeof(int*));

  for(i=0; i<N; i++) MALLOC(nn[i], KMAX*sizeof(int));

  MALLOC(state, N*sizeof(int));

 

  height = ceil(log(N)/log(2));

  MALLOC(level, height*sizeof(int));

  length = 0;

  level[0] = length;

  for(i=1; i<=height; i++){

    length += 1 << i;

    level[i] = length;

  }

  MALLOC(tree, length*sizeof(int));

  for(i=0; i<length; i++) tree[i] = 0;

 

  srand48(time(0));

  tessellation(L, DIM, KK, nn);

 

  avg = stdev = 0.0;

  for(l=0; l<ensemble; l++){

 

    remain = 0;

    for(i=0; i<N; i++){

      r = drand48();

      if(KK[i]<KMAX) state[i] = BOUNDARY;

      else if(r < prob1) state[i] = BLOCK;

      else if(r >= prob1 && r < prob1+prob2) state[i] = HIGHBLOCK;

      else{

        state[i] = UNEXAMINED;

        remain++;

        for(j=0; j<height; j++) tree[level[j] + (i >> (height-1-j))]++;

      }

    }

 

    if(DRAW) draw(state, L, 0, 0, remain, 0, 0);

 

    avail = 0;

    water = 0;

    while(remain){

      k = 0;

      for(i=0; i<height; i++){

        for(j=k; j<k+2; j++){

          if(tree[level[i] + j]){

            k = 2*j;

            break;

          }

        }

      }

      avail = j;

 

      state[avail] = WATER;

      size = 1;

      for(j=0; j<height; j++) tree[level[j] + (avail >> (height-1-j))]--;

 

      touch = 0;

      grow(avail, 1.0, KK, nn, state, &size, &touch, tree, level, height);

      remain -= size;

      if(touch==0) water += size;

 

      if(DRAW) draw(state, L, avail, size, remain, touch, water);

    }

 

 

    // second... height 2

    remain = 0;

    for(i=0; i<length; i++) tree[i] = 0;

    for(i=0; i<N; i++){

      if(state[i] == WATER || state[i] == BLOCK){

        remain++;

        for(j=0; j<height; j++) tree[level[j] + (i >> (height-1-j))]++;

      }

    }

 

    avail = 0;

    while(remain){

      k = 0;

      for(i=0; i<height; i++){

        for(j=k; j<k+2; j++){

          if(tree[level[i] + j]){

            k = 2*j;

            break;

          }

        }

      }

      avail = j;

 

      state[avail] = HIGHWATER;

      size = 1;

      for(j=0; j<height; j++) tree[level[j] + (avail >> (height-1-j))]--;

 

      touch = 0;

      regrow(avail, 1.0, KK, nn, state, &size, &touch, tree, level, height);

      remain -= size;

      if(touch==0) water += size;

 

      if(DRAW) draw(state, L, avail, size, remain, touch, water);

    }

 

    frac = (double)water/Ne;

    avg += frac;

    stdev += frac*frac;

  }

 

  avg /= ensemble;

  stdev /= ensemble;

  stdev = sqrt((stdev - avg*avg)/(ensemble));

  printf("%d %f %f %e %e\n", L, prob1, prob2, avg, stdev);

 

  free(KK);

  for(i=0; i<N; i++) free(nn[i]);

  free(nn);

  free(state);

  free(level);

  free(tree);

 

  return;

}

 

 

 

 

 

 

 

 

 

 

 

 

Johan Ofverstedt

Johhan Ofverstedt at Uppsala University did a student thesis on this subject.

http://uu.diva-portal.org/smash/record.jsf?pid=diva2:534020

He applied the computer science departments expertise in constraint programming to find solutions for maximum water retention on magic squares.  He posted the program on Sourceforge

http://sourceforge.net/projects/wrmssolver/

 Johan’s program can be easily modified to find solutions for all varieties of magic squares … ie semi-magic, associative magic squares etc.

 

Harry White

http://users.eastlink.ca/~sharrywhite/Download.html

Magic_Square_Graphic.png

25 x 25 magic square –

Harry White has been a force in the past few years  providing useful utilities.   He modified Johan Ofverstedt’s program and produced a “complete” utility that makes this terraformed graphic possible.

 

Albrecht Durer

1514.jpg

Albrecht Durer created this famous 4x4 associative magic square in 1514.  He placed the two numbers in the lower row 15  14  to mark the date of its creation.  The 14 x 14 magic square is a tribute to the upcoming 500th year anniversary.  This square contains 2014 units of water.   The numbers in red mark the date of Durer’s birth and death.

http://www.futilitycloset.com/2013/03/30/stormy-weather-3/

Greg Ross has a nice blog on this  ( link above)

 

 

 

Latest Information :

http://en.wikipedia.org/wiki/User_talk:Knecht03

 

history.jpg