// PCA.C
// DLL for Visual Basic 6 program, "Five Cellular Automata"
// Author: Peter Meyer
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "stdafx.h"
#define true 1
#define false 0
#define NUM_NEIGHBORS (8)
#define MAX_Q (256)
#define CONST_E (2.7182818)
#define MAX_POWER_E (709)
#define DISPLAY_SIZE_IN_PIXELS (528)
#define SWAP_STATES(i1,j1,i2,j2,state1,state2) { a1[i1][j1] = state2; a1[i2][j2] = state1; }
/*
Type passed from Visual Basic program:
Type PCAStruct
'Input:
lngType As Long
lngQ As Long
lngK1 As Long
lngK2 As Long
lngK3 As Long
lngK4 As Long
lngArraySize As Long
End Type
*/
typedef struct
{
int type;
int q;
int k1;
int k2;
int k3;
int k4;
int array_size;
} PCA;
unsigned char *d1;
unsigned char **a1;
unsigned char *d2;
unsigned char **a2;
int type, q, k1, k2, k3, k4, array_size, dynamics;
unsigned int num_bytes_in_row;
long sum, num_da_cells, num_fixed_da_cells;
int nnas[MAX_Q+1]; // holds numbers of immediate neighbours in all states
// Functions to export
__declspec(dllexport) int _stdcall allocate_memory(unsigned int max_array_size);
__declspec(dllexport) void _stdcall free_memory(void);
__declspec(dllexport) double _stdcall initialize_automaton(PCA *pca, int initial_state);
__declspec(dllexport) double _stdcall update_automaton(PCA *pca);
__declspec(dllexport) int _stdcall cell_state(int i, int j);
// Functions private to this module
int num_neighbors_in_state(int state, int i, int j, int *mm, int *nn, int array);
void num_neighbors_in_all_states(int i, int j, int array);
int sum_states(int i, int j, int include_cell);
double average_state(void);
double togetherness(void);
double proportion_healthy_cells(void);
double proportion_fixed_da_cells(void);
// Entry point for DLL
/*---------------------------------*/
BOOL APIENTRY DllMain(HANDLE hModule,
DWORD ul_reason_for_call,
LPVOID lpReserved)
{
return ( true );
}
/*-----------------------------------------------------*/
int _stdcall allocate_memory(unsigned int max_array_size)
{
unsigned int i;
unsigned char *ptr;
num_bytes_in_row = max_array_size;
if ( ( d1 = calloc(max_array_size*max_array_size,sizeof(char)) ) == NULL )
return ( false );
if ( ( a1 = malloc(max_array_size*sizeof(char *)) ) == NULL )
{
free(d1);
return ( false );
}
ptr = d1;
for ( i=0; i<max_array_size; i++ )
{
a1[i] = ptr;
ptr += max_array_size;
}
if ( ( d2 = calloc(max_array_size*max_array_size,sizeof(char)) ) == NULL )
{
free(a1);
free(d1);
return ( false );
}
if ( ( a2 = malloc(max_array_size*sizeof(char *)) ) == NULL )
{
free(a1);
free(d1);
free(d2);
return ( false );
}
ptr = d2;
for ( i=0; i<max_array_size; i++ )
{
a2[i] = ptr;
ptr += max_array_size;
}
srand(time(NULL));
return ( true );
}
/*---------------------------*/
void _stdcall free_memory(void)
{
free(a1);
free(d1);
free(a2);
free(d2);
}
// initial_state:
// 0 = random
// 1 = symmetric
// 2 = diagonal
/*-------------------------------------------------------------*/
double _stdcall initialize_automaton(PCA *pca, int initial_state)
{
int i, ii, j, jj, state, state0, m, k;
type = pca->type;
q = pca->q;
k1 = pca->k1;
k2 = pca->k2;
array_size = pca->array_size;
if ( type < 2 )
{
// q-state Life or BZ
switch ( initial_state )
{
case 0:
// random
for ( i=0; i<array_size; i++ )
{
for ( j=0; j<array_size; j++ )
a1[i][j] = (rand()%q) + 1;
}
break;
case 1:
// symmetric
// array_size is always even except for 33.
for ( i=0; i<(array_size+1)/2; i++ )
{
ii = array_size - i - 1;
for ( j=0; j<(array_size+1)/2; j++ )
{
jj = array_size - j - 1;
m = abs(j - i);
k = ( m*m*(rand()%q) + m*(rand()%q) + rand()%q)%q + 1;
a1[i][j] = a1[ii][j] = a1[i][jj] = a1[ii][jj] = k;
}
}
break;
case 2:
// diagonal
state0 = rand()%q;
for ( i=0; i<array_size; i++ )
{
state0 = (state0%q) + 1;
state = state0;
for ( j=0; j<array_size; j++ )
a1[i][j] = (state++%q) + 1;
}
break;
}
}
else if ( type == 2 )
{
// Togetherness
// c = (double)k1/100; // concentration of cells
for ( i=0; i<array_size; i++ )
{
for ( j=0; j<array_size; j++ )
{
if ( rand()*(unsigned long)100 > k1*(unsigned long)RAND_MAX )
a1[i][j] = 1; // no cell at this location
else
a1[i][j] = (rand()%(q-1)) + 2;
// Cells have states in the range 2 through q.
}
}
}
else if ( type == 3 )
{
// 50% of locations are occupied by cells.
// All cells are initially healthy.
for ( i=0; i<array_size; i++ )
{
for ( j=0; j<array_size; j++ )
{
if ( rand()%2 )
a1[i][j] = 1; // empty
else
a1[i][j] = q; // healthy cell
}
}
}
else if ( type == 4 )
{
num_da_cells = 0;
num_fixed_da_cells = 0;
memset(d1,0,(unsigned long)DISPLAY_SIZE_IN_PIXELS*DISPLAY_SIZE_IN_PIXELS);
// Seed cells are in state q.
// Of other locations, k1% are occupied by (invisible) mobile cells in state 1;
// the rest are empty (state 0).
m = array_size*array_size/4;
if ( k2 >= 6 )
{
// Place k2 seeds at random.
do {
i = rand()%array_size;
j = rand()%array_size;
ii = array_size/2 - i;
jj = array_size/2 - j;
if ( ii*ii + jj*jj <= m )
{
a1[i][j] = q;
num_da_cells++;
num_fixed_da_cells++;
}
} while ( num_fixed_da_cells < k2 );
}
for ( i=0; i<array_size; i++ )
{
for ( j=0; j<array_size; j++ )
{
if ( k2 >= 6 )
if ( a1[i][j] == q )
continue;
ii = array_size/2 - i;
jj = array_size/2 - j;
if ( ii*ii + jj*jj > m )
{
a1[i][j] = 0; // empty location
continue;
}
if ( k2 < 6 )
{
if ( ( k2 == 1 &&
( i == array_size/2 && j == array_size/2 ) )
|| ( k2 == 2 &&
( ( i == array_size/4 && j == array_size/2 )
|| ( i == 3*array_size/4 && j == array_size/2 ) ) )
|| ( k2 == 3 &&
( ( i == array_size/4 && j == array_size/2 )
|| ( i == array_size/2 + (int)(array_size*sqrt(3)/8) && j == 3*array_size/4 )
|| ( i == array_size/2 + (int)(array_size*sqrt(3)/8) && j == array_size/4 ) ) )
|| ( k2 == 4 &&
( ( i == array_size/2 && j == array_size/2 )
|| ( i == array_size/4 && j == array_size/2 )
|| ( i == array_size/2 + (int)(array_size*sqrt(3)/8) && j == 3*array_size/4 )
|| ( i == array_size/2 + (int)(array_size*sqrt(3)/8) && j == array_size/4 ) ) )
|| ( k2 == 5 &&
( ( i == array_size/2 && j == array_size/2 )
|| ( i == array_size/4 && j == array_size/4 )
|| ( i == array_size/4 && j == 3*array_size/4 )
|| ( i == 3*array_size/4 && j == array_size/4 )
|| ( i == 3*array_size/4 && j == 3*array_size/4 ) ) )
)
{
a1[i][j] = q;
num_da_cells++;
num_fixed_da_cells++;
}
else
{
if ( rand()*(unsigned long)100 < k1*(unsigned long)RAND_MAX )
{
a1[i][j] = 1; // mobile cell
num_da_cells++;
}
else
a1[i][j] = 0; // empty location
}
}
else // k2 >= 6, place mobile cells at random locations.
{
if ( rand()*(unsigned long)100 < k1*(unsigned long)RAND_MAX )
{
a1[i][j] = 1; // mobile cell
num_da_cells++;
}
else
a1[i][j] = 0; // empty location
}
}
}
}
if ( type <= 1 )
return ( average_state() );
else if ( type == 2 )
return ( togetherness() );
else if ( type == 3 )
return ( proportion_healthy_cells() );
else // if ( type == 4 )
return ( proportion_fixed_da_cells() );
}
/*--------------------------------------*/
double _stdcall update_automaton(PCA *pca)
{
int i, j, k, m, n, sxs;
int adjacent, gain1, gain2, gain;
int i1, i2, i3, j1, j2, j3, nns11, nns12, nns21, nns22;
int state, state1, state2;
int nns1, nnsq;
type = pca->type;
q = pca->q;
k1 = pca->k1;
k2 = pca->k2;
k3 = pca->k3;
k4 = pca->k4;
array_size = pca->array_size;
if ( type <= 1 )
#if false
// Copy 'array_size' rows of array a1 to array a2.
memcpy(d2,d1,num_bytes_in_row*array_size);
// Array a2 not used in Togetherness.
#endif
{
for ( i=0; i<array_size; i++ )
memcpy(d2+i*num_bytes_in_row,d1+i*num_bytes_in_row,array_size);
}
switch ( type )
{
case 0: // q-state Life
for ( i=0; i<array_size; i++ )
{
for ( j=0; j<array_size; j++ )
{
sum = sum_states(i,j,false);
// false = don't include state of cell itself
if ( a1[i][j] > q/2 )
{
if ( k1 <= sum && sum <= k2 )
a1[i][j] += 1;
else
a1[i][j] -= 1;
}
else
{
if ( k3 <= sum && sum <= k4 )
a1[i][j] += 1;
else
a1[i][j] -= 1;
}
if ( a1[i][j] < 1 )
a1[i][j] = 1;
else if ( a1[i][j] > q )
a1[i][j] = q;
}
}
break;
case 1: // Belouzov-Zhabotinsky
for ( i=0; i<array_size; i++ )
{
for ( j=0; j<array_size; j++ )
{
if ( a1[i][j] == 1 )
{
nns1 = num_neighbors_in_state(1,i,j,NULL,NULL,2);
if ( k1 == k2 )
state = ( NUM_NEIGHBORS - nns1 )/k1 + 1;
else
{
nnsq = num_neighbors_in_state(q,i,j,NULL,NULL,2);
state = ( NUM_NEIGHBORS - nns1 - nnsq )/k1 + nnsq/k2 + 1;
// 1 <= k1,k2 <= 8
}
if ( state > q )
state = q;
a1[i][j] = state;
}
else if ( a1[i][j] == q )
a1[i][j] = 1;
else
{
nns1 = num_neighbors_in_state(1,i,j,NULL,NULL,2);
// This also sets the global variable 'sum'
state = sum/(NUM_NEIGHBORS - nns1 + 1) + k3;
// k3 >= 1 so no need to check if state < 1.
if ( state > q )
state = q;
a1[i][j] = state;
}
}
}
break;
case 2: // Togetherness
sxs = array_size*array_size;
for ( k=0; k<sxs; k++ )
{
// Randomly choose a cell which is not surrounded by 8 neighbors of the same state.
do {
do {
i1 = rand()%array_size;
j1 = rand()%array_size;
} while ( a1[i1][j1] == 1 );
state1 = a1[i1][j1];
num_neighbors_in_all_states(i1,j1,1);
nns11 = nnas[state1];
} while ( nns11 == NUM_NEIGHBORS );
if ( k2 >= array_size )
{
// Choose another location at random.
do {
i2 = rand()%array_size;
j2 = rand()%array_size;
} while ( i2 == i1 && j2 == j1 );
}
else
{
// Choose at random a location within a distance k3
// which is other than the location of the cell itself.
do {
i3 = k2 - rand()%(2*k2+1); // -k2 <= i3,j3 <= k2
j3 = k2 - rand()%(2*k2+1);
} while ( ( i3 == 0 && j3 == 0 ) || ( abs(i3) + abs(j3) > k2 ) );
i2 = ( i1 + i3 + array_size )%array_size;
j2 = ( j1 + j3 + array_size )%array_size;
}
// If there is a cell at this location of the same state then give up.
if ( ( state2 = a1[i2][j2] ) == state1 )
continue;
nns12 = nnas[state2];
num_neighbors_in_all_states(i2,j2,1);
nns21 = nnas[state1];
adjacent = ( abs(i1-i2) <= 1 && abs(j1-j2) <= 1 );
if ( state2 == 1 )
{
// 2nd location is empty.
gain = ( nns21 - adjacent ) - nns11;
}
else
{
// 2nd location is occupied by a cell.
nns22 = nnas[state2];
if ( nns22 == NUM_NEIGHBORS )
continue;
// Check if exchange of states would lead to increase in adjacent cells with same state.
gain1 = ( nns21 - adjacent ) - nns11;
gain2 = ( nns12 - adjacent ) - nns22;
gain = gain1 + gain2;
}
if ( gain >= 0 )
SWAP_STATES(i1,j1,i2,j2,state1,state2)
else if ( ( state2 == 1 && gain == -1 ) || ( state2 != -1 && gain >= -2 ) )
{
if ( !(rand()%k3 ) )
SWAP_STATES(i1,j1,i2,j2,state1,state2)
}
}
break;
case 3: // Viral Replication
sxs = array_size*array_size;
for ( k=0; k<sxs; k++ )
{
// Choose a location at random
i = rand()%array_size;
j = rand()%array_size;
// If empty (state 1) do nothing
if ( a1[i][j] == 1 )
continue;
else if ( a1[i][j] == q )
{
// It's a healthy cell.
// It becomes infected with probability k2 chances in 100,000.
// If it does not become infected then
// if it is not completely surrounded by other cells
// then it divides with probability k3 percent,
// with one daughter cell occupying a random neighboring empty location.
if ( rand()*(unsigned long)100000 < k2*(unsigned long)RAND_MAX )
a1[i][j]--;
else
{
if ( num_neighbors_in_state(1,i,j,&i1,&j1,1) > 0 )
{
if ( rand()*(unsigned long)100 < k3*(unsigned long)RAND_MAX )
a1[i1][j1] = q;
}
}
}
else if ( a1[i][j] > 2 && a1[i][j] < q )
{
// If it's infected (but not fully infected) then its state decreases by 1.
a1[i][j]--;
}
else if ( a1[i][j] == 2 )
{
// If the cell is fully infected (state 2) then it disappears and all neighboring cells
// are infected with probability k1 percent.
a1[i][j] = 1;
for ( i1=-1; i1<=1; i1++ )
{
for ( j1=-1; j1<=1; j1++ )
{
if ( ! ( i1 == 0 && j1 == 0 ) )
{
i2 = ( i + i1 + array_size)%array_size;
j2 = ( j + j1 + array_size)%array_size;
if ( a1[i2][j2] > 2 )
{
if ( rand()*(unsigned long)100 < k1*(unsigned long)RAND_MAX )
a1[i2][j2]--;
}
}
}
}
}
}
break;
case 4: // Diffusion/aggregation
// Going through the array typewriter-fashion is a lot faster but it produces
// a definite bias in one direction, so must choose locations at random.
sxs = array_size*array_size;
i3 = sxs/4;
// Square of radius of largest circle within square display.
// Movable cells which pass beyond this disappear.
for ( k=0; k<sxs; k++ )
{
// Choose a location at random
i = rand()%array_size;
j = rand()%array_size;
// If not a movable cell at this location do nothing.
if ( a1[i][j] == 1 )
{
// If it has at least one fixed cell as a neighbor then it becomes a fixed cell.
// Unfortunately this code allows the structure to grow beyond the square boundary
// but not worth correcting.
adjacent = 0;
for ( m=-1; m<=1; m++ )
{
for ( n=-1; n<=1; n++ )
{
if ( ! ( m == 0 && n == 0 ) )
{
if ( a1[(i+m+array_size)%array_size][(j+n+array_size)%array_size] > 1 )
{
adjacent = 1;
break;
}
}
}
if ( adjacent == 1 )
break;
}
if ( adjacent == 1 )
{
num_fixed_da_cells++;
// set state of cell based on square root of proportion of fixed cells so that
// the central cells are white and there is gradual movement through the available colors
// the number of which is determined by the value of q.
a1[i][j] = (int)( ( q*q - 2 - q*(q-2)*sqrt(proportion_fixed_da_cells()) ) / (q-1) );
#if false
// Not needed.
if ( a1[i][j] > q )
a1[i][j] = q;
else if ( a1[i][j] < 2 )
a1[i][j] = 2;
#endif
}
else
{
// Move to a random empty neighboring location if possible
// (impossible only if all neighboring locations are occupied by movable cells).
if ( num_neighbors_in_state(0,i,j,&i1,&j1,1) > 0 )
{
a1[i][j] = 0;
// If cell crosses boundary then it disappears.
#if false
// This produces a squarish final image.
if ( ( i == 0 && i1 == array_size -1 ) || ( i1 == 0 && i == array_size -1 )
|| ( j == 0 && j1 == array_size -1 ) || ( j1 == 0 && j == array_size -1 ) )
#endif
// This gives a circular final image.
i2 = array_size/2 - i1;
j2 = array_size/2 - j1;
if ( i2*i2 + j2*j2 >= i3 )
num_da_cells--;
else
a1[i1][j1] = 1;
}
}
}
}
break;
}
if ( type <= 1 )
return ( average_state() );
else if ( type == 2 )
return ( togetherness() );
else if ( type == 3 )
return ( proportion_healthy_cells() );
else // if ( type == 4 )
return ( proportion_fixed_da_cells() );
}
// This used by Diffusion-Aggregation.
/*----------------------------------*/
double proportion_fixed_da_cells(void)
{
return ( (double)num_fixed_da_cells/num_da_cells );
}
// This used by Togetherness.
/*-----------------------------------------------------*/
void num_neighbors_in_all_states(int i, int j, int array)
{
int m, n;
memset(nnas,0,sizeof(nnas));
for ( m=-1; m<=1; m++ )
{
for ( n=-1; n<=1; n++ )
{
if ( ! ( m == 0 && n == 0 ) )
{
if ( array == 1 )
nnas[a1[(i+m+array_size)%array_size][(j+n+array_size)%array_size]]++;
else
nnas[a2[(i+m+array_size)%array_size][(j+n+array_size)%array_size]]++;
}
}
}
}
// Returns the number of immediate neighbours of cell at i,j in specified state using specified array.
// Also sets 'sum' to the sum of the states of the cells' neighbors plus the state of the cell itself.
// If there is a neighbor in state 'state' then on return *mm and *nn are the relative coordinates
// of a neighbor in that state, randomly selected.
// This used by BZ, VI and DA.
/*----------------------------------------------------------------------------*/
int num_neighbors_in_state(int state, int i, int j, int *i1, int *j1, int array)
{
int m, n, mm, nn, nns=0, neighbor_state, k1, k2;
sum = ( array == 1 ? a1[i][j] : a2[i][j] );
for ( m=-1; m<=1; m++ )
{
for ( n=-1; n<=1; n++ )
{
if ( ! ( m == 0 && n == 0 ) )
{
mm = (i+m+array_size)%array_size;
nn = (j+n+array_size)%array_size;
neighbor_state = ( array == 1 ? a1[mm][nn] : a2[mm][nn] );
sum += neighbor_state;
nns += ( neighbor_state == state );
}
}
}
if ( i1 != NULL && nns > 0 )
{
// Select at random a neighbor in state 'state'.
k1 = rand()%nns;
k2 = 0;
for ( m=-1; m<=1; m++ )
{
for ( n=-1; n<=1; n++ )
{
if ( ! ( m == 0 && n == 0 ) )
{
mm = (i+m+array_size)%array_size;
nn = (j+n+array_size)%array_size;
neighbor_state = ( array == 1 ? a1[mm][nn] : a2[mm][nn] );
if ( neighbor_state == state && k2++ == k1 )
{
*i1 = mm;
*j1 = nn;
return ( nns );
}
}
}
}
}
return ( nns );
}
// Returns the sum of the states of the cell's neighbours
// using periodic boundary conditions, and optionally includes
// the state of the cell itself. NB: Uses array2[][].
// This used only by q-state Life.
/*------------------------------------------*/
int sum_states(int i, int j, int include_cell)
{
int m, n;
sum = 0;
for ( m=-1; m<=1; m++ )
{
for ( n=-1; n<=1; n++ )
sum += a2[(i+m+array_size)%array_size][(j+n+array_size)%array_size];
}
return ( include_cell ? sum : sum - a2[i][j] );
}
// This used only by Togetherness.
/*---------------------*/
double togetherness(void)
{
int i, j, state, n=0, sum=0;
for ( i=0; i<array_size; i++ )
{
for ( j=0; j<array_size; j++ )
{
if ( ( state = a1[i][j] ) > 1 )
{
n++;
num_neighbors_in_all_states(i,j,1);
sum += nnas[state];
}
}
}
return ( (double)sum/n );
}
// This used only by q-state Life and BZ.
/*----------------------*/
double average_state(void)
{
int i, j;
sum = 0;
for ( i=0; i<array_size; i++ )
{
for ( j=0; j<array_size; j++ )
sum += a1[i][j];
}
return ( (double)sum/(array_size*array_size) );
}
// This used only by Viral Infection
/*---------------------------------*/
double proportion_healthy_cells(void)
{
long i, j, num_cells=0, num_healthy_cells=0;
for ( i=0; i<array_size; i++ )
{
for ( j=0; j<array_size; j++ )
{
if ( a1[i][j] > 1 )
{
num_cells++;
if ( a1[i][j] == q )
num_healthy_cells++;
}
}
}
if ( !num_cells )
return ( 0 );
else
return ( (double)num_healthy_cells/num_cells );
}
/*---------------------------------*/
int _stdcall cell_state(int i, int j)
{
return ( a1[i][j] );
}