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