//  EQUIL.C
//  Functions concerned with equilibration
//  © 2021 Peter J. Meyer

#include "iss.h"

#define CHECK_MAG_ON_SPIN_RESTORATION  false
#define DISABLE_ESCAPE_KEY false
#define SHOW_WARNING false
#define DISPLAY_TRANSITION_PROBABILITIES false
#define DISPLAY_INITIAL_MAGNETIZATIONS   false
#define DISPLAY_MAP_FILE false

static int less_screen_output=false;
static int first_configuration_timed=false;
static int abs_magnetization_0, abs_internal_energy_0;

extern int num_spins_visited_this_run;      //  WOLFF.C

static void perform_equilibrations_for_configuration(int i);
static int perform_equilibrations_for_spin_assignment(int i, int j);
static void perform_one_equilibration(void);
static int escape_key_hit(int exit_if_escape_key_hit);
static void just_a_sec(void);

/*  
For each configuration
	Set up configuration (i.e. occupy sites and open bonds)
	For each spin assignment
		Assign initial spins
		(to same configuration if second or later spin assignment)
		If more than one repetition
			Save spin assignment and absolute magnetization
            (i.e. save sites[...] with initial spin assignment)
		End If
		For each repetition
			If second or later repetition
				Restore spin assignment and absolute magnetization
				(i.e. load saved sites[...] with initial spin assignment) 
			End If
			Do the equilibration and take measurements
		End For
	End For
End For
*/

/*-----------------------------*/
void perform_equilibrations(void)
{
int i;
clock_t start_clock, end_clock;
double time_for_first_configuration, time_elapsed;

if ( dynamics_is_metropolis || dynamics_is_glauber )
    calculate_transition_probabilities();       //  TRANSIT.C
else    //  dynamics is Swendsen-Wang or Wolff
    {
    if ( model_is_ising )
        //  sw_probability = wolff_probability = 1 - exp(-2*J/temperature); 
        v_bond_probability = 1 - exp(-2*J/temperature);  
    else    //  model is q-state Potts
        {
        #if SHOW_WARNING
        printf("\nSwendsen-Wang or Wolff algorithm "
               "with q-state Potts needs checking.  See %s.\n",__FILE__);
        #endif
        v_bond_probability = 1 - exp(-J/(temperature));
        //  See Wang & Swendsen, Physica A, 167, 565.
        }
    }

num_spin_visits = num_spin_visits_overflow = 0;
num_spin_flips = num_spin_flips_overflow = 0;
max_cluster_size = num_clusters_traced = 0;
num_sweeps = num_time_slices_done = 0;

zero_measurements();        //  MEAS.C

//  num_configurations is 1 in the case of a non-diluted (pure) lattice.
for ( i=0; i<num_configurations; i++ )
    {
    if ( num_configurations > 1 )
        {
        if ( !first_configuration_timed )
            start_clock = clock();

        printf("\n");
        if ( num_input_files > 1 )
            {
            printf("Input file %s",input_filepath);
            printf(" (%d/%d), ",input_file_num+1,num_input_files);
            }
        if ( multiple_temperatures_specified )
            printf("Temperature %d/%d, ",temperature_num+1,num_temperatures);
        printf("Configuration %d/%d",i+1,num_configurations);

        if ( i > 0 )
            {
            time_elapsed = (( (double)clock() - task_start ) / CLOCKS_PER_SEC );
            printf("\nTime elapsed this temperature = ");
            if ( time_elapsed < 60 )
                printf("%.2f seconds",time_elapsed);
            else   
                printf("%s",seconds_to_time_str((unsigned long)time_elapsed,temp));
            }
        }

    //  Occupy sites (all unless site-diluted)
    //  and open bonds (all between occupied sites unless bond-diluted).
    set_up_configuration();       //  S&B.C
    //  All spins are up-spins.

    perform_equilibrations_for_configuration(i);
    if ( escape_key_termination )
        break;

    if ( num_configurations > 1 )
        {
        if ( !first_configuration_timed )
            {
            end_clock = clock();
            time_for_first_configuration 
                = ( (double)end_clock - start_clock ) / CLOCKS_PER_SEC ;
            printf("Time for first configuration = %.2f seconds\n",
                time_for_first_configuration);
            if ( time_for_first_configuration < 30 )
                less_screen_output = true;
            first_configuration_timed = true;
            }
        }
    }

#if SAVE_MEASUREMENTS
//  Save measurements before analysis is done, since this changes 
//  values in the arrays used to hold the results of measurement.
printf("\n%d bytes written to %s.",
    save_measurements(),measurements_filepath);     //  FILE_IO.C
#endif

analyse_data();  
}

//  i is configuration number
/*-------------------------------------------------------*/
static void perform_equilibrations_for_configuration(int i) 
{
int j, k, map_file_displayed=false;

for ( j=0; j<num_spin_assignments; j++ )
    {
    if ( !less_screen_output )
        {
        printf("\n");
        if ( num_input_files > 1 )
            printf("%s (%d/%d)",input_filepath,input_file_num+1,num_input_files);
        if ( multiple_temperatures_specified )
            printf(" Temp. %d/%d.",temperature_num+1,num_temperatures);
        if ( num_configurations > 1 )
            printf(" Config'n %d/%d.",i+1,num_configurations);
        if ( num_spin_assignments > 1 )
            printf(" Spin ass't %d/%d",j+1,num_spin_assignments);
        if ( num_repetitions > 1 )
            printf(" (%d rep'ns)",num_repetitions);
        printf(": ");
        }

    //  For second spin assignment it is not the case that
    //  all spins are up (as on first assignment) so we need
    //  to re-assign the spins even if initial_magnetization is 1.0.
    assign_initial_spins();  //  INITSPIN.C   
    //  This function sets abs_magnetization and abs_internal_energy.

    if ( autocorrelation_measured || num_repetitions > 1 )
        copy_spin_assignment();     //  to initial_sites[...]
    
    #if DISPLAY_INITIAL_MAGNETIZATIONS
    printf("\nSpin model initialized with magnetization %f.\n",
        magnetization_per_spin());
    #endif

    if ( escape_key_hit(false) || escape_key_termination )
        break;        

    //  This is where the map file is written, if one is requested.
    //  A map file can be generated anytime by calling write_map().
    if ( map_file_requested && !map_file_displayed )
        {
        display_map(map_filepath,0);
        map_file_displayed = true;
        }

    if ( num_repetitions > 1 )
        {
        abs_magnetization_0 = abs_magnetization;
        if ( internal_energy_measured )
            abs_internal_energy_0 = abs_internal_energy;
        }

    k = perform_equilibrations_for_spin_assignment(i,j);

    if ( escape_key_termination )
        break;
    }

printf("\n");

if ( escape_key_termination )
    {
    num_samples = (double)i*num_spin_assignments*num_repetitions
        + j*num_repetitions + k + 1;
    printf("Terminated ");
    #if COMPLETE_SAMPLE_ON_ESCAPE
    printf("at end of ");
    #else
    printf("during ");
    #endif
    if ( multiple_temperatures_specified )
        printf("temperature #%d of %d (%.4f),\n",
            temperature_num+1,num_temperatures,temperatur[temperature_num]);
    printf("config'n %d, spin ass't %d, "
        "rep'n %d, MCS %d.\n",i+1,j+1,k+1,mcs_done);
    }
else
    num_samples = (double)num_configurations*num_spin_assignments*num_repetitions;
}

/*---------------------------------------------------------------*/
static int perform_equilibrations_for_spin_assignment(int i, int j)
                                                    //  i is configuration number
{                                                   //  j is spin assignment number
int k;
//  int time_slice_num;
clock_t this_clock, last_clock=0;

for ( k=0; k<num_repetitions; k++ )                 //  k is repetition number
    {
    if ( num_repetitions > 1 && !less_screen_output )
        {
        if ( k == 0 )
            {
            printf("Rep.1   ");
            last_clock = clock();
            }
        else
            {
            this_clock = clock();
            if ( this_clock > last_clock + CLOCKS_PER_SEC )
                {
                printf("Rep.%-4d",k+1);
                last_clock = this_clock;
                }
            }
        }

    if ( k > 0 )
        {
        restore_spin_assignment();                       //  SPINS.C

        //  And restore abs_magnetization.
        abs_magnetization = abs_magnetization_0;
        if ( internal_energy_measured )
            abs_internal_energy = abs_internal_energy_0;

        #if CHECK_MAG_ON_SPIN_RESTORATION
        if ( get_abs_magnetization(UP_SPIN) != abs_magnetization_0 )
            {
            printf("\nSpin assignment restoration failed.");
            printf("\nconfig'n %d/%d, spin ass't %d/%d, "
                "rep'n %d/%d",i+1,num_configurations,
                j+1,num_spin_assignments,k+1,num_repetitions);
            printf("\nget_abs_magnetization(UP_SPIN) = %d",
                   get_abs_magnetization(UP_SPIN));
            printf("\nabs_magnetization_0 = %d\n",abs_magnetization_0);
            exit(1);
            }
        #endif
        }

    mcs_num = 0;
    perform_one_equilibration();        //  below
    mcs_done = mcs_num;

    if ( escape_key_termination || escape_key_hit(false) )
        break;
    }

return ( k );   //  used only if escape key termination.
}

/*---------------------------------------*/
static void perform_one_equilibration(void)
{
int time_slice_num=0;
unsigned long num_mcs_to_do = (num_time_slices-1)*step_length;
double magn;

num_spins_visited_this_run = 0;     //  for use with Wolff algorithm

//  Once through the loop is a Monte Carlo step per spin, a.k.a. a timestep.
loop
    {
    if ( !(mcs_num%step_length) )  
        {
        update_measurements(time_slice_num++);  //  MEAS.C

        #if DISPLAY_MAP_FILE
        write_map(map_filepath,true,true);
        getch();
        //  If DISPLAY_MAP_FILE is true and ISM is run in release mode 
        //  without this getch() something bad will happen.
        #endif
        }

    if ( mcs_num == num_mcs_to_do )
        break;
    
    #if !COMPLETE_SAMPLE_ON_ESCAPE
    //  This is #defined in ISS_DEF.H
    if ( escape_key_hit(false) )
        break;
    #endif

    num_spins_visited_this_sweep = num_spins_flipped_this_sweep
        = num_clusters_traced_this_sweep = 0;
                
    if ( dynamics_is_metropolis || dynamics_is_glauber )  
        {
        if ( spin_seln_is_random )
            do_random_ssf_sweep();          //  SINGLE.C
        else
            do_checkerboard_ssf_sweep();    //  SINGLE.C
        //  Random single spin flip is significantly more time-consuming
        //  than is the checkerboard, due to need to generate random numbers.
        }
    else    //  dynamics is Swendsen-Wang or Wolff
        {
        if ( dynamics_is_swendsen_wang )
            perform_swendsen_wang_sweep();      //  SW-WANG.C
        else    //  dynamics is wolff
            perform_wolff_sweep();              //  WOLFF.C
        }

    if ( num_spin_visits + num_spins_visited_this_sweep < num_spin_visits )
        num_spin_visits_overflow++;

    num_spin_visits += num_spins_visited_this_sweep;

    if ( num_spin_flips + num_spins_flipped_this_sweep < num_spin_flips )
        num_spin_flips_overflow++;

    num_spin_flips += num_spins_flipped_this_sweep;

    if ( ++num_sweeps == 0 )
        num_sweeps_overflow++;

    if ( dynamics_is_swendsen_wang || dynamics_is_wolff )
        {
        if ( num_clusters_traced + num_clusters_traced_this_sweep 
            < num_clusters_traced )
            num_clusters_traced_overflow++;

        num_clusters_traced += num_clusters_traced_this_sweep;
        }

    mcs_num++;

    magn = magnetization_per_spin();
    if ( magn < -1.0 || magn > 1.0 )
        {
        printf("\nError! magnetization = %f abs_magnetization = %d,"
            "\nmcs_num = %d  num_spins = %d\n",
            magn,abs_magnetization,mcs_num,num_spins);
        exit(1);
        }
    }       //  end loop 
}

//  This also handles pause/resume and accummulates pause time.
//  Currently always called with exit_if_escape_key_hit = false.
/*-------------------------------------------------*/
static int escape_key_hit(int exit_if_escape_key_hit)
{
int ch, paused=false;
clock_t pause_start;

if ( kbhit() )
    {
    ch = getch();
    if ( ch == SPACE )
        {
        printf("\nPaused  ...  ");
        paused = true;
        pause_start = clock();
        //  Check the keyboard once per second.
        while ( !kbhit() )
            just_a_sec();
        ch = getch();
        }
    if ( paused )
        pause_time += (unsigned long)(( clock() - pause_start ) / CLOCKS_PER_SEC );
    if ( ch == CONTROL_T )
        stop_on_next_temperature = true;
        //  Not control-S, since that causes pause in screen output.
    #if !DISABLE_ESCAPE_KEY
    if ( ch == ESCAPE )
        {
        if ( exit_if_escape_key_hit )
            {
            printf("\nEscape key hit. Program terminated ");
            system_date(run_end_date,ISO);
            system_time(run_end_time,true);
            printf(" at %s %s",run_end_date,run_end_time);
            if ( num_input_files > 1 )
                printf(" while processing %s",input_filepath);
            printf(".\n");   
            exit(1);
            }
        printf("\nEscape key pressed.");
        return ( escape_key_termination = true );
        }
    #endif
    if ( paused )
        printf("Resumed    ");
    }

return ( false );
}

//  Pause for one second.
/*------------------------*/
static void just_a_sec(void)
{
clock_t start = clock();

while ( clock() < start + CLOCKS_PER_SEC )
    ;
}