ca.c

细胞自动机的一个源代码

#include "ca.h"

int (*ca_update_func)(ca_t *ca, int s, int e);


/* get flat position in memory array from point
   vector p */
int ca_point(ca_t *ca, int *p) {
  int i;
  int run=1, c, n=0;

  c = ca->cur;

  for (i=0; i<(ca->dim); i++) {
    n += run*p[i];
    run *= ca->d[i];
  }
  return(n);
}

/* get point vector from flat memroy point p
   and put into dest */
int ca_vec(ca_t *ca, int *dest, int p) {
  static int *d = NULL, *run = NULL, dim = 0;
  int i, j, k;

  if (dim!=(ca->dim)) {

    if (d) free(d);
    d = (int *)malloc(sizeof(int)*(ca->dim));
    if (run) free(run);
    run = (int *)malloc(sizeof(int)*(ca->dim));
    dim = ca->dim;
  }

  run[0] = 1;
  for (i=0; i<(dim-1); i++) {
    run[i+1] = ca->d[i]*run[i];
  }

  for (i=(dim-1); i>=0; i--) {
    dest[i] = p/run[i];
    p -= (p/run[i])*run[i];
  }
}  

int *ca_nei_addp(ca_t *ca, int *p, int *a) {
  int i, j, k;
  int carry=0;

  for (i=0; i<(ca->dim); i++) {
    p[i] += a[i] + ca->d[i];
    p[i] %= ca->d[i];

    /*
    carry = ( (p[i]>ca->d[i]) ? 1 : 0 );

    if (p[i]>=ca->d[i]) {
      carry = p[i] - ca->d[i] + 1;
    } else if (p[i]<0) {
      carry = ca->d[i] + p[i];
    } else carry=0;

    p[i] += ca->d[i];
    p[i] %= ca->d[i];
    */
  }
  return(p);

}

/* get sum of neighbors for outer totalistic
   CA */
int ca_tot_sum(ca_t *ca, int p) {
  int i, j, k=0;
  int c;
  int tfp;
  static int *po=NULL, *t=NULL, dim=0;

  if (dim!=ca->dim) {
    dim = ca->dim;
    if (po) {
      free(po);
      free(t);
    }
    po = (int *)malloc(sizeof(int)*(ca->dim));
    t = (int *)malloc(sizeof(int)*(ca->dim));
  }

  c = ca->cur;
  ca_vec(ca, t, p);

  for (i=0; i<(ca->nei_num); i++) {
    memcpy(po, t, sizeof(int)*dim);
    ca_nei_addp(ca, po, ca->nei[i]);
    tfp = ca_point(ca, po);

    k += (ca->ca[c][tfp])*(ca->w[i]);

    /* whoopsie, doesn't take into acount edges for > 1 dim */
    /*
    k += (ca->ca[c][(p + ca->fnei[i] + ca->fd)%ca->fd])*(ca->w[i]);
    */

  }

  /*
  if (ca->ca[c][p]) {
    printf("k: %i (%i)\n", k, ca->wadd);
  }
  */

  /* k is neighbor sum.  wadd is sum of weights.
     since its outer totalistic, add the middle cell
     value * sum of weights to get rule */
  return(k + (ca->ca[c][p]*(ca->wadd)));

}

/* get next state for cell.
   rule position = p_n p_(n-1) ... p_0 p */
int ca_flat_state(ca_t *ca, int p) {
  int i, j, k, c, n, *d;
  int dig=0, digp=1, tfp;
  static int *po=NULL, *t=NULL, dim=0;

  if (dim!=ca->dim) {
    dim = ca->dim;
    if (po) {
      free(po);
      free(t);
    }
    po = (int *)malloc(sizeof(int)*(ca->dim));
    t = (int *)malloc(sizeof(int)*(ca->dim));
  }

  c = ca->cur;
  ca_vec(ca, t, p);

  for (i=(ca->nei_num-1); i>=0; i--) {

    memcpy(po, t, sizeof(int)*dim);
    ca_nei_addp(ca, po, ca->nei[i]);
    tfp = ca_point(ca, po);

    dig += ca->ca[c][tfp]*digp;
    digp *= ca->k;

    /* whoops, doesn't handle edges in > 1 dimension */
    /*
    dig += ca->ca[c][ (ca->fnei[i] + p + ca->fd) % ca->fd ]*digp;
    digp *= ca->k;
    */

  }

  return(ca->r[dig]);
}

/* get rule position from flat position in memory */
int ca_flat_rule(ca_t *ca, int p) {
  int i, j, k, c, n, *d;
  int dig, digp=1;

  c = ca->cur;

  dig = 0;

  for (i=0; i<(ca->nei_num); i++) {
    dig += ca->ca[c][ (ca->fnei[i] + p + ca->fd) % ca->fd ]*digp;
    digp *= ca->dim;
  }

  /*
  dig += ca->ca[c][p];
  digp *= ca->dim;

  for (i = ((ca->nei_num)/2); i<(ca->nei_num); i++) {
    dig += ca->ca[c][ca->fnei[i]+p]*digp;
    digp *= ca->dim;
  }
  */

  return(dig);
}

int ca_tot_update(ca_t *ca, int s, int e) {
  int i, j, c, n;

  c = ca->cur;
  n = (ca->cur + 1)%(ca->wlen);

  for (i=s; i<=e; i++) {
    ca->ca[n][i] = ca->tr[ca_tot_sum(ca, i)];
  }

  ca->cur++;
  ca->cur%=ca->wlen;

  return(0);
}

/*the way i'm doing it, there really isn't 
  any use of having s and e...*/
int ca_flat_update(ca_t *ca, int s, int e) {
  int i, j, c, n;

  c = ca->cur;
  n = (ca->cur+1)%(ca->wlen);

  /* apply the rule to every cell */
  for (i=s; i<=e; i++) {
    ca->ca[n][i] = ca_flat_state(ca,i);
  }

  ca->cur++;
  ca->cur%=ca->wlen;

  return(0);

}

/* a more user friendly version of update.
   ps start point, pe end point
   DON'T USE (its not right)
*/
int ca_update(ca_t *ca, int *ps, int *pe) {

  return(ca_flat_update(ca, ca_point(ca, ps), ca_point(ca, pe)));
}


/*********************************
 memory functions (alloc, init, free etc.)
**********************************/

/* allocate the cell aut for the ca
   assumed ca->dim & ca->d allocated and initialized */
int ca_cell_alloc(ca_t *ca) {
  int i, n=1;

  /* allocate XxYxZ... size ca */
  for (i=0; i<(ca->dim); i++) {
    n*=ca->d[i];
  }


  //printf("ca: allocating ca (%i)\n", sizeof(int *)*(ca->wlen));

  if (!(ca->ca = (int **)malloc(sizeof(int *)*(ca->wlen))))
    return(-1);
  for (i=0; i<(ca->wlen); i++) {


    //printf("\t%i (%i)\n", i, sizeof(int)*n);

    if (!(ca->ca[i] = (int *)malloc(sizeof(int)*n)))
      return(-1);
  }

  ca->cur = 0;
  ca->next = 1;

  return(0);
}



ca_t *ca_alloc(int k, int dim, int wlen,
	       int *d, int nei_num, int **nei,
	       int *r, float **b, float *wb) {
  ca_t *ca;
  ca = (ca_t *)malloc(sizeof(ca_t));
  return(ca_init(ca, k, dim, wlen, d, nei_num, nei, r, b, wb));
}


/* non totalistic init function */
ca_t *ca_init(ca_t *ca, int k, int dim, int wlen,
	      int *d, int nei_num, int **nei,
	      int *r, float **b, float *wb) {
  int i, j;
  /*
    ca_t *ca;
    ca = (ca_t *)malloc(sizeof(ca_t));
  */

  ca->k = k;
  ca->dim = dim;
  ca->wlen = wlen;
  ca->slen = wlen;

  ca->go = 0;

  /* dimensions of ca */

  //printf("ca: allocating d (%i)\n", sizeof(int)*dim);

  if (!(ca->d = (int *)malloc(sizeof(int)*dim)))
    return(NULL);

  ca->fd=1;
  for (i=0; i<dim; i++) {
    ca->d[i] = d[i];
    ca->fd *= d[i];
  }


  ca->nei_num = nei_num;

  /* neigbors */


  //printf("ca: allocating nei & fnei (%i)\n", sizeof(int *)*nei_num);

  if (!(ca->nei = (int **)malloc(sizeof(int *)*nei_num)))
    return(NULL);
  if (!(ca->fnei = (int *)malloc(sizeof(int *)*nei_num)))
    return(NULL);


  /* aloccate and assign relative neighbors */
  for (i=0; i<nei_num; i++) {

    //printf("\t%i (%i)\n", i, sizeof(int)*dim);

    if (!(ca->nei[i] = (int *)malloc(sizeof(int)*dim)))
      return(NULL);
    for (j=0; j<dim; j++) {
      ca->nei[i][j] = nei[i][j];
    }
  }

  /* assign flat neighbors */
  for (i=0; i<nei_num; i++) {
    ca->fnei[i] = ca_point(ca, ca->nei[i]);
  }


  /* rule lookup */
  ca->r_num = 1;
  for (i=0; i<nei_num; i++) {
    ca->r_num *= k;
  }

  //printf("ca: allocating r (%i)\n", sizeof(int)*(ca->r_num));

  if (!(ca->r = (int *)malloc(sizeof(int)*(ca->r_num))))
    return(NULL);
  for (i=0; i<(ca->r_num); i++) {
    ca->r[i] = r[i];
  }

  /* basis */

  //printf("ca: allocating b (%i)\n", sizeof(float *)*dim);

  if (!(ca->b = (float **)malloc(sizeof(float *)*dim)))
    return(NULL);
  for (i=0; i<dim; i++) {

    //printf("ca:\t%i (%i)\n", i, sizeof(float)*4);

    if (!(ca->b[i] = (float *)malloc(sizeof(float)*4)))
      return(NULL);
    for (j=0; j<dim; j++) {
      ca->b[i][j] = b[i][j];
    }
    for (; j<4; j++) {
      ca->b[i][j] = 0.0;
    }
  }

  /* win vector */

  //printf("ca: allocating wb (%i)\n", sizeof(float)*(dim+1));

  if (!(ca->wb = (float *)malloc(sizeof(float)*(dim+1))))
    return(NULL);
  for (i=0; i<(dim+1); i++) {
    ca->wb[i] = wb[i];
  }

  ca->delay = 1000;

  ca->cur = 0;
  ca->next = 1;

  if (ca_cell_alloc(ca)<0) {
    return(NULL);
  }

  for (j=0; j<ca->wlen; j++) {
    for (i=0; i<ca->fd; i++) {
      ca->ca[j][i] = 0;
    }
  }

  ca->diff = 0;
  ca->slen = 1;

  return(ca);

}

/* outer totalistic init */
ca_t *ca_tot_init(ca_t *ca, int k, int dim, int wlen,
		  int *d, int nei_num, int **nei, int *w,
		  int *r, float **b, float *wb) {
  int i, j;
  /*
    ca_t *ca;
    ca = (ca_t *)malloc(sizeof(ca_t));
  */

  ca->k = k;
  ca->dim = dim;
  ca->wlen = wlen;
  ca->slen = wlen;

  ca->go = 0;

  /* dimensions of ca */
  if (!(ca->d = (int *)malloc(sizeof(int)*dim)))
    return(NULL);

  ca->fd=1;
  for (i=0; i<dim; i++) {
    ca->d[i] = d[i];
    ca->fd *= d[i];
  }


  ca->nei_num = nei_num;

  /* neigbors */
  if (!(ca->nei = (int **)malloc(sizeof(int *)*nei_num)))
    return(NULL);
  if (!(ca->fnei = (int *)malloc(sizeof(int *)*nei_num)))
    return(NULL);


  /* aloccate and assign relative neighbors */
  for (i=0; i<nei_num; i++) {
    if (!(ca->nei[i] = (int *)malloc(sizeof(int)*dim)))
      return(NULL);
    for (j=0; j<dim; j++) {
      ca->nei[i][j] = nei[i][j];
    }
  }

  /* weights */
  ca->w = (int *)malloc(sizeof(int)*nei_num);
  for (i=0; i<nei_num; i++) {
    ca->w[i] = w[i];
  }

  ca->wadd=1;
  for (i=0; i<ca->nei_num; i++) {
    ca->wadd += ca->w[i];
  }

  /* assign flat neighbors */
  for (i=0; i<nei_num; i++) {
    ca->fnei[i] = ca_point(ca, ca->nei[i]);
  }


  /* outer totalistic rule lookup */
  ca->tr_num = (nei_num+1)*k;
  if (!(ca->tr = (int *)malloc(sizeof(int)*(ca->tr_num))))
    return(NULL);
  for (i=0; i<(ca->tr_num); i++) {
    ca->tr[i] = r[i];
  }

  /* basis */
  if (!(ca->b = (float **)malloc(sizeof(float *)*dim)))
    return(NULL);
  for (i=0; i<dim; i++) {
    if (!(ca->b[i] = (float *)malloc(sizeof(float)*dim)))
      return(NULL);
    for (j=0; j<dim; j++) {
      ca->b[i][j] = b[i][j];
    }
  }

  /* win vector */
  if (!(ca->wb = (float *)malloc(sizeof(float)*(dim+1))))
    return(NULL);
  for (i=0; i<(dim+1); i++) {
    ca->wb[i] = wb[i];
  }

  ca->delay = 1000;

  ca->cur = 0;
  ca->next = 1;

  if (ca_cell_alloc(ca)<0) {
    return(NULL);
  }

  ca->slen = 1;
  ca->diff = 0;

  return(ca);

}

void ca_free(ca_t *ca) {
  int i, j, k;


  //printf("ca: freeing d\n");

  free(ca->d);

  for (i=0; i<ca->nei_num; i++) {

    //printf("ca:\t%i\n", i);

    free(ca->nei[i]);
  }

  //printf("ca: freeing nei & fnei\n");

  free(ca->nei);
  free(ca->fnei);

  for (i=0; i<ca->wlen; i++) {

    //printf("ca:\t%i\n", i);

    free(ca->ca[i]);
  }

  //printf("ca: freeing ca\n");

  free(ca->ca);

  if (ca->w) {

    //printf("ca: freeing w\n");

    free(ca->w);
  }
  for (i=0; i<ca->dim; i++) {

    //printf("ca:\t%i\n", i);

    free(ca->b[i]);
  }

  //printf("ca: freeing b\n");

  free(ca->b);

  if (ca->wb) {

    //printf("ca: freeing wb\n");

    free(ca->wb);
  }
  if (ca->r) {

    //printf("ca: freeing r\n");

    free(ca->r);
  }

  if (ca->tr) {

    //printf("ca: freeing tr\n");

    free(ca->tr);
  }

  //printf("ca: freeing mother ca\n");

  free(ca);
}

/**************************************