// PERC.C
// Determination of percolation thresholds
// © 2021 Peter J. Meyer
#include "iss.h"
#define NUM_ANOMALIES (5)
#define INITIAL_NUM_CONFIGS (500)
#define NUM_CONFIGURATONS_FACTOR (2)
#define MAX_NUM_CONFIGURATIONS (100000)
#define MAX_REPEATS 5
#define INITIAL_EPSILON (8)
#define EPSILON_FACTOR (2)
#define W 13
double min_p, max_p, precision2;
int total_num_configs;
char anomalies[NUM_ANOMALIES+1][64] =
{
/* 0 */ "",
/* 1 */ "must have first percentage < 50",
/* 2 */ "must have final percentage > 50",
/* 3 */ "a percentage decrease occurred",
/* 4 */ "percentage within epsilon of 50",
/* 5 */ "too many repeats for same range"
};
/* Outline of algorithm:
min_p = 0
max_p = 1
While ( max_p - min_p > precision )
For p=min_p, p<=max_p, P += (max_p-min_p)/5
Generate multiple configurations
and get percentage with spanning cluster.
End For
If percentages not increasing (or other anomaly)
repeat
Else
min_p = largest p with percentage < 50
max_p = smallest p with percentage > 50
End If
End While
*/
/*-------------------------------*/
int get_percolation_threshold(void)
{
int iter, p_num, i, level, num_times_exceeded_50;
int num_with_spanning_cluster, halfway_check_done;
int map_displayed;
double p, epsilon, range;
iter = -1;
level = 0;
min_p = 0.0;
max_p = 1.0;
epsilon = INITIAL_EPSILON;
num_configurations = INITIAL_NUM_CONFIGS;
total_num_configs = 0;
// free_boundary_conditions set to true in SETPARAM.C.
do {
level++;
map_displayed = false;
do {
if ( ++iter > 0 )
{
if ( iter == MAX_ITERATIONS )
goto done;
if ( iterations[iter-1].repeat )
printf("\nRepeating because %s.",anomalies[iterations[iter-1].anom_code]);
else
{
if ( iterations[iter-1].anom_code == 4 )
printf("\nAdjusting range because %s.",anomalies[4]);
epsilon /= EPSILON_FACTOR;
}
// To avoid percentage decreases at the later stages
// we need to increase the number of configurations.
num_configurations = (int)(num_configurations*NUM_CONFIGURATONS_FACTOR);
if ( num_configurations > MAX_NUM_CONFIGURATIONS )
num_configurations = MAX_NUM_CONFIGURATIONS;
}
printf("\nIteration #%d, level #%d, epsilon = %.*f",
iter+1,level,num_decimal_places+1,epsilon);
// num_decimal_places read from input file ("precision").
range = max_p - min_p;
printf("\nRange %.*f - %.*f, %.*f in steps of %.*f, %s configurations.\n",
num_decimal_places+1,min_p,
num_decimal_places+1,max_p,
num_decimal_places+1,range,
num_decimal_places+1,range/(NUM_CONCENTRATIONS-1),
ultoa_commas(num_configurations,temp));
// Copy data to array of structs for use
// when writing results to disk file.
iterations[iter].level = level;
iterations[iter].max_p = max_p;
iterations[iter].min_p = min_p;
iterations[iter].range = range;
iterations[iter].num_configurations = num_configurations;
iterations[iter].epsilon = epsilon;
for ( p_num=0; p_num<NUM_CONCENTRATIONS; p_num++ )
{
iterations[iter].perc_data[p_num][0] = p
= min_p + p_num*(range/(NUM_CONCENTRATIONS-1));
iterations[iter].perc_data[p_num][1] = 200; // 200 = no measurement made
printf("%*.*f",W,num_decimal_places+1,p);
}
printf("\n");
iterations[iter].repeat = false;
iterations[iter].anom_code = 0;
num_times_exceeded_50 = 0;
for ( p_num=0; p_num<NUM_CONCENTRATIONS; p_num++ )
{
p = iterations[iter].perc_data[p_num][0];
if ( goal_is_site_percolation_threshold )
{
site_concentration = p; // bond_concentration set in SETPARAM.C
site_diluted = ( site_concentration < 1.0 );
}
else // goal is bond_percolation;
{
bond_concentration = p; // site_concentration set in SETPARAM.C
bond_diluted = ( bond_concentration < 1.0 );
}
if ( p == 0.0 )
iterations[iter].perc_data[p_num][1] = 0.0; // never a spanning cluster
else
{
// Even if p = 1.0 there may not be a spanning cluster
// (e.g. if seeking bond percolation threshold with low site concentration).
num_with_spanning_cluster = 0;
halfway_check_done = false;
for ( i=0; i<num_configurations; i++ )
{
// Occupy (some or all) sites and open (some or all) bonds.
set_up_configuration(); // S&B.C
total_num_configs++;
ct_data.cluster_type = CT_OPEN_BOND;
ct_data.stop_on_spanning_cluster = true;
trace_all_clusters(&ct_data); // Look for spanning cluster.
if ( ct_data.spanning_cluster_found )
{
num_with_spanning_cluster++;
if ( map_file_requested )
{
if ( !map_displayed )
{
display_map(map_filepath,1);
map_displayed = true;
}
}
}
// If we are halfway through and percentage < 40 then give up.
if ( !halfway_check_done )
{
if ( i >= num_configurations/2 )
{
if ( ((double)num_with_spanning_cluster*200)
/num_configurations < 40 )
{
num_with_spanning_cluster *= 2;
break;
}
halfway_check_done = true;
}
}
// Check for Escape key after each 100 configurations.
if ( !(i%100) )
if ( kbhit() )
if ( getch() == ESCAPE )
{
printf("\nEscape pressed.");
escape_key_termination = true;
goto done;
}
}
iterations[iter].perc_data[p_num][1] =
((double)num_with_spanning_cluster*100)/num_configurations;
}
if ( iterations[iter].perc_data[p_num][1] == 0 )
printf("%*d",W,0);
else if ( iterations[iter].perc_data[p_num][1] == 100 )
printf("%*d",W,100);
else
printf("%*.2f",W,iterations[iter].perc_data[p_num][1]);
// Check that results show no anomalies.
if ( p_num == 0 && iterations[iter].perc_data[p_num][1] >= 50 )
{
iterations[iter].repeat = true;
iterations[iter].anom_code = 1;
break;
// because must have first percentage < 50
}
if ( p_num == NUM_CONCENTRATIONS - 1
&& iterations[iter].perc_data[p_num][1] <= 50 )
{
iterations[iter].repeat = true;
iterations[iter].anom_code = 2;
break;
// because must have final percentage > 50
}
if ( p_num > 0 )
{
if ( iterations[iter].perc_data[p_num][1]
< iterations[iter].perc_data[p_num-1][1] )
{
iterations[iter].repeat = true;
iterations[iter].anom_code = 3;
break;
// because a percentage decrease occurred
}
}
num_times_exceeded_50 += ( iterations[iter].perc_data[p_num][1] > 50 );
// If the percentage has exceeded 50% twice
// or has reached 80% then no need to continue;
if ( num_times_exceeded_50 == 2 || iterations[iter].perc_data[p_num][1] >= 80 )
break;
}
printf("\n");
if ( iter >= MAX_REPEATS-1 && iterations[iter].repeat )
{
if ( iterations[iter].level == iterations[iter-MAX_REPEATS+1].level )
{
// About to do a further repeat after MAX_REPEATS
// for the same range, so give up.
iterations[iter].anom_code = 5;
printf("\nQuit because %s.\n",anomalies[5]);
goto done;
}
}
} while ( iterations[iter].repeat );
// At this point percolation_data[iter][0][1] < 50 and
// percolation_data[iter][NUM_CONCENTRATIONS-1][1] > 50.
// Set min_p to the largest p with percentage less than 50%.
p_num = 1;
while ( iterations[iter].perc_data[p_num][1] < 50 )
p_num++;
min_p = iterations[iter].perc_data[--p_num][0];
// If percentage for min_p is too close to 50 then widen range.
if ( 50 - iterations[iter].perc_data[p_num][1] < epsilon )
{
// Subtract step length from min_p.
min_p -= range/(NUM_CONCENTRATIONS-1);
if ( min_p < 0 )
min_p = 0;
iterations[iter].anom_code = 4;
}
// Set max_p to the smallest p with percentage greater than 50%.
p_num = NUM_CONCENTRATIONS - 2;
while ( iterations[iter].perc_data[p_num][1] > 50 )
p_num--;
max_p = iterations[iter].perc_data[++p_num][0];
// If percentage for max_p is too close to 50 then widen range.
if ( iterations[iter].perc_data[p_num][1] - 50 < epsilon )
{
max_p += range/(NUM_CONCENTRATIONS-1);
if ( max_p > 1 )
max_p = 1;
iterations[iter].anom_code = 4;
}
} while ( max_p - min_p > precision );
done:
precision2 = max_p - min_p;
// The final value is ( max_p + min_p )/2 +/- precision2.
if ( precision < 0.5 )
{
printf("\nThe %s percolation threshold for this lattice"
"\n(with %s concentration %.3f)",
( goal_is_site_percolation_threshold ? "site" : "bond" ),
( goal_is_site_percolation_threshold ? "bond" : "site" ),
( goal_is_site_percolation_threshold ?
bond_concentration : site_concentration ));
printf("\nis %.*f +/- %.*f (%.*f - %.*f)\n",
num_decimal_places+1, (max_p+min_p)/2,
num_decimal_places+1, precision2,
num_decimal_places+1, (max_p+min_p)/2-precision2,
num_decimal_places+1, (max_p+min_p)/2+precision2);
}
printf("\n%s configurations generated.\n",ultoa_commas(total_num_configs,temp));
return ( ++iter );
}