//  RESULTS.C
//  © 2021 Peter J. Meyer

#include "iss.h"

#include <process.h>

//  #define WR_ST_DEV_AUTOCORRELATON     Now in ISS_DEF.H
#define WR_LOG_ST_DEV_MAGNETISM false
//  #define ALWAYS_WRITE_AUTOCORRELATION true
#define WRITE_SYSTEM_DATA false

static double mean_magn, std_dev_magn;
static double mean_bc, std_dev_bc;
static double mean_ie, std_dev_ie;
static double mean_sh, std_dev_sh;

#if TIMING_DIAGNOSTICS
extern long create_v_cluster_clock;
extern long flip_sw_clusters_clock;
extern long flip_sw_cluster_clock;
extern long another_sw_spin_clock;
extern long flip_sw_spin_clock;
#endif

extern unsigned int num_times_adjust_index_called;
extern unsigned int num_times_get_nn_sites_called;
extern unsigned int num_times_get_nn_site_in_dir_called;
extern unsigned int num_virtual_bonds_opened;
extern unsigned int total_num_opened_bonds;

static void write_lattice_info(FILE *of);
static void write_samples_info(FILE *of, int multiple_temperatures_summary);
static void calculate_magnetization_results(int pc, FILE *of);
static void calculate_binder_cumulant(int pc, FILE *of);
static void calculate_ie_sh_results(int pc, FILE *of);

extern void write_timeslice_values(FILE *of, int w);       //  TS_VALS.C

/*----------------------------------*/
void write_equilibration_results(void)
{
int i, j;
int w=9;
int ch;
int temp_written=false;
double crit_temp;
FILE *of;
#if WRITE_SYSTEM_DATA
unsigned int nmcs;
#endif

while  ( !open_file(output_filepath,"a+t",&of) )
    //  "a+", append, in case output from multiple runs with a range of temperatures.
    //  file already truncated (if previously existing) in set_up_filenames() in FILE_IO.C.
    {
    printf("\n\aCannot open output file %s.",output_filepath);
    printf("\nPerhaps this file is already open?  Retry or quit? (R|Q) ");
    ch = getche();
    if ( ch == 'Q' || ch == 'q' ) 
        exit(1);
    }

fprintf(of,"Run began at %s %s\n",run_start_date,run_start_time);

if ( escape_key_termination )
    fprintf(of,"\nEscape key was pressed: run probably not completed.\n");

if ( stop_on_next_temperature && ( temperature_num+1 != num_temperatures ) )
    fprintf(of,"\nRequest was made to stop on next temperature."
               "\nSome temperatures not done.\n");

//  This run data is very similar to that given by display_startup_info().
fprintf(of,"\ninput file:         %s",input_filepath);
if ( num_input_files > 1 )
    fprintf(of," (#%d of %d)",input_file_num+1,num_input_files);
fprintf(of,"\noutput file:        %s",output_filepath);
#if SAVE_MEASUREMENTS
fprintf(of,"\nmeasurements file:  %s",measurements_filepath);
#endif

fprintf(of,"\n");

write_lattice_info(of);

sprintf(temp,"%.6f",temperature);
remove_trailing_zeros(temp);
fprintf(of,"\ntemperature:        %s ",temp);

//  If this is the critical temperature of the current model then show this.

if ( ( crit_temp = get_critical_temperature() ) != -1 )
    {
    if ( temperature == crit_temp )
        {
        if ( model_is_ising )
            fprintf(of,"(Ising %s cr. temp.)",lattice_types[l_type].name);
        else    //  model is q-state Potts
            fprintf(of,"(%dd %d-state Potts cr. temp.)",dimensionality,q);
        }
    else
        fprintf(of," = %.3f%% of Tc",100*temperature/crit_temp);
    }

if ( multiple_temperatures_specified )
    fprintf(of," (#%d of %d)",temperature_num+1,num_temperatures);

fprintf(of,"\ndynamics:           %s",dynamics_types[d_type][1]);

if ( dynamics_is_swendsen_wang || dynamics_is_wolff )
    fprintf(of," (p = %.5f)",v_bond_probability);
else
    fprintf(of,"\nspin selection:     %s",spin_selns[ss_type][1]);

fprintf(of,"\ninitial magnetization:      %.5f\n",initial_magnetization);

write_samples_info(of,false);

if ( absolute_magnetization_requested || !timeslice_values_requested )
    fprintf(of,"\n");

if ( absolute_magnetization_requested && timeslice_values_requested )
    fprintf(of,"\nAbsolute magnetization requested.");

if ( timeslice_values_requested )
    write_timeslice_values(of,w);       //  TS_VALS.C

calculate_magnetization_results(percentage_range_for_mean,of);
if (  multiple_temperatures_specified )
    {
    multiple_temperature_results[temperature_num][0][0] = mean_magn;
    multiple_temperature_results[temperature_num][0][1] = std_dev_magn;
    }

if ( binder_cumulant_measured )
    {
    calculate_binder_cumulant(percentage_range_for_mean,of);
    if (  multiple_temperatures_specified )
        {
        multiple_temperature_results[temperature_num][1][0] = mean_bc;
        multiple_temperature_results[temperature_num][1][1] = std_dev_bc;
        }
    }

if ( internal_energy_measured )
    {
    calculate_ie_sh_results(percentage_range_for_mean,of);    
    if (  multiple_temperatures_specified )
        {
        multiple_temperature_results[temperature_num][2][0] = mean_ie;
        multiple_temperature_results[temperature_num][2][1] = std_dev_ie;
        multiple_temperature_results[temperature_num][3][0] = mean_sh;
        multiple_temperature_results[temperature_num][3][1] = std_dev_sh;
        }
    }

fprintf(of,"\nRun ended at %s %s",run_end_date,run_end_time);

fprintf(of,"\nRun time = ");
if ( run_time < 60 )
    fprintf(of,"%.2f seconds",run_time);
else   
    {
    fprintf(of,"%s",seconds_to_time_str((unsigned long)run_time,temp));
    fprintf(of," (%lu seconds)",(unsigned long)run_time);
    }

if ( pause_time )
    fprintf(of,"\nexcluding pause time of %s seconds.",ultoa_commas(pause_time,temp));

#if TIMING_DIAGNOSTICS

fprintf(of,"\nTime spent in cluster-creation function = %.2f seconds (%.2f%%)",
    (double)create_v_cluster_clock/CLOCKS_PER_SEC,
    100*(double)create_v_cluster_clock/CLOCKS_PER_SEC/run_time);

fprintf(of,"\nTime spent in clusters-flip function = %.2f seconds (%.2f%%)",
    (double)flip_sw_clusters_clock/CLOCKS_PER_SEC,
    100*(double)flip_sw_clusters_clock/CLOCKS_PER_SEC/run_time);

fprintf(of,"\nTime spent in cluster-flip function = %.2f seconds (%.2f%%)",
    (double)flip_sw_cluster_clock/CLOCKS_PER_SEC,
    100*(double)flip_sw_cluster_clock/CLOCKS_PER_SEC/run_time);

fprintf(of,"\nTime spent in flip_sw_spin function = %.2f seconds (%.2f%%)",
    (double)flip_sw_spin_clock/CLOCKS_PER_SEC,
    100*(double)flip_sw_spin_clock/CLOCKS_PER_SEC/run_time);

fprintf(of,"\nTime spent in another_sw_spin = %.2f seconds (%.2f%%)",
    (double)another_sw_spin_clock/CLOCKS_PER_SEC,
    100*(double)another_sw_spin_clock/CLOCKS_PER_SEC/run_time);

fprintf(of,"\n");

#endif

#if WRITE_SYSTEM_DATA

if ( dynamics_is_swendsen_wang )
    {
    fprintf(of,"\nAverage number of SW virtual bonds opened per sweep per sample = %s",
        ultoa_commas((unsigned long)(num_virtual_bonds_opened/num_sweeps/num_samples),temp));
    fprintf(of,"\nAverage percentage of open bonds becoming SW virtual bonds = %.2f%%",
        100*(double)num_virtual_bonds_opened/num_sweeps/num_samples/total_num_opened_bonds);
    //  NEEDS CHECKING
    }

fprintf(of,"\nNumber of times get_nn_sites() called = %s",
    ultoa_commas(num_times_get_nn_sites_called,temp));
fprintf(of,"\nNumber of times get_nn_site_in_dir() called = %s",
    ultoa_commas(num_times_get_nn_site_in_dir_called,temp));
fprintf(of,"\nNumber of times adjust_index() called = %s",
    ultoa_commas(num_times_adjust_index_called,temp));

fprintf(of,"\n");

nmcs = (num_time_slices-1)*step_length + 1;

if ( !num_spin_visits && !num_spin_visits_overflow )
    fprintf(of,"\nNo spins were visited.");
else 
    {
    if ( !num_spin_visits_overflow )
        {
        fprintf(of,"\n%.2f%% of %s spin visits resulted in spin flips (%s spin flips).",
            (num_spin_flips*(double)100)/num_spin_visits,
            ultoa_commas((unsigned long)num_spin_visits,temp),
            ultoa_commas((unsigned long)num_spin_flips,temp2));
        fprintf(of,"\nSpin flips per spin per sample per timestep = %.4f",
            (double)num_spin_flips/num_samples/num_spins/nmcs);
        if ( run_time > 0 )
            fprintf(of,"\nNumber of spin visits per second = %s",
                ultoa_commas((unsigned long)(num_spin_visits/run_time),temp));
        }
    else
        {
        fprintf(of,"\n%.2f%% of %.0f million spins visited were flipped "
            "(%.0f million spins flipped).",
            ((num_spin_flips+pow(2,32)*num_spin_flips_overflow)*100)/
            (num_spin_visits+pow(2,32)*num_spin_visits_overflow),
            (num_spin_visits+pow(2,32)*num_spin_visits_overflow)/1e6,
            (num_spin_flips+pow(2,32)*num_spin_flips_overflow)/1e6);
        fprintf(of,"\nSpin flips per spin per sample per timestep = %.4f",
            (num_spin_flips+pow(2,32)*num_spin_flips_overflow)
            /num_samples/num_spins/nmcs);
        fprintf(of,"\nNumber of spin visits per second = %s",ultoa_commas(
                (unsigned long)((num_spin_visits+pow(2,32)*num_spin_visits_overflow)
                /run_time),temp));
        }

    if ( dynamics_is_swendsen_wang || dynamics_is_wolff )
        {
        fprintf(of,"\nTotal number of clusters traced = ");
        if ( num_clusters_traced_overflow )
            fprintf(of,"%.0f",num_clusters_traced+pow(2,32)*num_clusters_traced_overflow);
        else
            fprintf(of,"%s",ultoa_commas((unsigned long)num_clusters_traced,temp));

        fprintf(of,"\nMaximum cluster size = %s ",
            ultoa_commas((unsigned long)max_cluster_size,temp));
        fprintf(of,"(%.2f%% of occupied sites).",
            ((double)max_cluster_size*100)/num_spins);

        fprintf(of,"\nAverage size of clusters traced = ");
        fprintf(of,"%.4f",
            (num_spin_visits+pow(2,32)*num_spin_visits_overflow)
            /(num_clusters_traced+pow(2,32)*num_clusters_traced_overflow));

        fprintf(of,"\nAverage number of clusters per sweep = %.4f",
            (num_clusters_traced+pow(2,32)*num_clusters_traced_overflow)
            /(num_sweeps+pow(2,32)*num_sweeps_overflow));

        if ( dynamics_is_wolff )
            fprintf(of,"\nAverage number of spins visited per sweep per spin = %.4f",
                (num_spin_visits/(num_sweeps+pow(2,32)*num_sweeps_overflow))
                /num_spins);
        }

    if ( stack_size )
        write_stack_data(of);       //  STACK.C
    }

#endif

if ( num_spin_visits || num_spin_visits_overflow )
    {
    if ( (unsigned long)NUM_PRNG1_CALLS() < 2e9 )
        fprintf(of,"\n%s was called %s times.",PRNG1_NAME,
            ultoa_commas((unsigned long)NUM_PRNG1_CALLS(),temp));
    else
        fprintf(of,"\n%s was called %s million times.",PRNG1_NAME,
            ultoa_commas((unsigned long)(NUM_PRNG1_CALLS()/1e6),temp));

    if ( seed < 100000 )
        fprintf(of,"   Random seed = %u",seed);

    fprintf(of,"\n");

    if ( PRNG1_PERIOD_EXCEEDED() )
        fprintf(of,"\nThe period for %s was exceeded.\n",PRNG1_NAME);
    }

if ( multiple_temperatures_specified )
    fprintf(of,"\n--------------------------------------------------------\n\n");

if ( multiple_temperatures_specified && 
    ( final_temperature || escape_key_termination || stop_on_next_temperature ) )
    {
    write_lattice_info(of);
    fprintf(of,"\ndynamics: %s",dynamics_types[d_type][1]);
    write_samples_info(of,true);
    //  Print multiple results.
    fprintf(of,"\n     %-*s%-*s",w+1,"Temp.",2*(w+1),"Magnetization");
    if ( binder_cumulant_measured ) 
        fprintf(of,"%-*s",2*(w+1),"Binder cumulant");
    if ( internal_energy_measured ) 
        fprintf(of,"%-*s%-*s",2*(w+1),"Internal Energy",2*(w+1),"Specific heat");

    fprintf(of,"\n%*s",w+1,"");
    fprintf(of,"%*s%*s",w+1,"mean",w+1,"std dev");
    if ( binder_cumulant_measured )            
        fprintf(of,"%*s%*s",w+1,"mean",w+1,"std dev");
    if ( internal_energy_measured ) 
        for ( j=2; j<=3; j++ )           
            fprintf(of,"%*s%*s",w+1,"mean",w+1,"std dev");

    for ( i=0; i<num_temperatures; i++ )
        {
        fprintf(of,"\n%3d%*.6f",i+1,w+1,temperatur[i]);
        fprintf(of,"%*.6f",w+1,multiple_temperature_results[i][0][0]);
        if ( multiple_temperature_results[i][0][1] == NON_EXISTENT )
            fprintf(of,"%*s",w+1,"");
        else
            fprintf(of,"%*.6f",w+1,multiple_temperature_results[i][0][1]);

        if ( binder_cumulant_measured )            
            fprintf(of,"%*.6f%*.6f",w+1,multiple_temperature_results[i][1][0],
                                    w+1,multiple_temperature_results[i][1][1]);
        if ( internal_energy_measured ) 
            {
            for ( j=2; j<=3; j++ )           
                fprintf(of,"%*.6f%*.6f",w+1,multiple_temperature_results[i][j][0],
                                        w+1,multiple_temperature_results[i][j][1]);
            }
        }
    fprintf(of,"\n");
    }

fclose(of);

if ( final_temperature || escape_key_termination || stop_on_next_temperature )
    spawnlp(P_NOWAIT,"notepad.exe","notepad.exe",output_filepath,NULL);
}

/*------------------------------------*/
static void write_lattice_info(FILE *of)
{
fprintf(of,"\nmodel:              ");

if ( model_is_q_state_potts )
    fprintf(of,"%d-state Potts",q);
else
    fprintf(of,"%s",models[m_type][1]);    

fprintf(of,"\ndimensionality:     %d",dimensionality);
fprintf(of,"\nlattice type:       %s",lattice_types[l_type].name);
if ( direction_table_is_parity_dependent && MAJOR_AXIS_SPECIFICATION_PERMITTED )
    fprintf(of,"\nmajor axis:         %d",major_axis);
fprintf(of,"\nlattice size:       %d",size);

if ( site_diluted )
    fprintf(of,"\nsite concentration: %.3f",site_concentration);
else if ( bond_diluted )
    fprintf(of,"\nbond concentration: %.3f",bond_concentration);

fprintf(of,"\n");
}

/*-----------------------------------------------------------------------*/
static void write_samples_info(FILE *of, int multiple_temperatures_summary)
{
double num_requested;

fprintf(of,"\nnumber of configurations:   %d",num_configurations);
fprintf(of,"\nnumber of spin assignments: %d",num_spin_assignments);
fprintf(of,"\nnumber of repetitions:      %d",num_repetitions);   
fprintf(of,"\nnumber of samples:          ");
if ( num_samples < pow(2,31) )
    fprintf(of,"%s",ultoa_commas((unsigned long)num_samples,temp));
else
    fprintf(of,"%.0f",num_samples); 

if ( escape_key_termination )
    {
    if ( multiple_temperatures_summary )
        fprintf(of,"\nEscape key pressed during temperature #%d",temperature_num+1);
    else
        {
        num_requested = (double)num_configurations*num_spin_assignments*num_repetitions;
        if ( num_requested < pow(2,31) )
            fprintf(of," (of %s requested)",ultoa_commas((unsigned long)num_requested,temp));
        else
            fprintf(of," (of %.0f requested)",num_requested);
        }
    }
    
fprintf(of,"\n\nnumber of time slices:      %d",num_time_slices); 
fprintf(of,"\nstep length:                %d",step_length);
}

//  Average the mean and the standard deviation
//  of the magnetization over the last pc% of time slices.
/*---------------------------------------------------------*/
static void calculate_magnetization_results(int pc, FILE *of)
{
double sum_magn, sum_std_dev;
int j, pr;
int to_ts = num_time_slices - 1;
int from_ts = 1 + ((num_time_slices-1)*(100-pc))/100;
int num_values = to_ts-from_ts+1;

if ( from_ts > to_ts - 1 )
    return;

fprintf(of,
    "\nMean and standard deviation of the %smagnetization"
    "\nover the last %d%% of timeslices, i.e., timeslices %d through %d.",
    ( absolute_magnetization_requested ? "absolute " : "" ),
    pc,from_ts,to_ts);

fprintf(of,"\n%s per spin",
    ( absolute_magnetization_requested ? "Absolute magnetization" : "Magnetization" ));
fprintf(of,"\n%-9s%-9s%-9s%-9s","mean","std.dev.","m-sd","m+sd");

sum_magn = sum_std_dev = 0.0;
for ( j=from_ts; j<=to_ts; j++ )
    {
    sum_magn += magnetization[j][0];
    if ( magnetization[j][1] == NON_EXISTENT )
        std_dev_magn = NON_EXISTENT;
    else if ( std_dev_magn != NON_EXISTENT )
        sum_std_dev += magnetization[j][1];
    }

mean_magn = sum_magn/num_values;
if ( std_dev_magn != NON_EXISTENT )
    std_dev_magn = sum_std_dev/num_values;

for ( pr=5; pr>=2; pr--)
    {
    fprintf(of,"\n%-9.*f",pr,mean_magn);
    if ( std_dev_magn != NON_EXISTENT )
        fprintf(of,"%-9.*f%-9.*f%-9.*f",pr,std_dev_magn,
        pr,mean_magn-std_dev_magn,pr,mean_magn+std_dev_magn);
    }

fprintf(of,"\n");
}

//  Average the binder cumulant over the last pc% of time slices
//  and estimate the std. dev.
/*---------------------------------------------------*/
static void calculate_binder_cumulant(int pc, FILE *of)
{
double sum_bc, sum_sq_dev_bc;
int j, pr;
int to_ts = num_time_slices - 1;
int from_ts = 1 + ((num_time_slices-1)*(100-pc))/100;
int num_values = to_ts-from_ts+1;

if ( from_ts > to_ts - 1 )
    return;

fprintf(of,"\nMean and standard deviation of the Binder cumulant"
           "\nover the last %d%% of timeslices, i.e., timeslices %d through %d.",
           pc,from_ts,to_ts);

fprintf(of,"\nBinder cumulant");
fprintf(of,"\n%-9s%-9s%-9s%-9s","mean","std.dev.","m-sd","m+sd");

sum_bc = 0.0;
for ( j=from_ts; j<=to_ts; j++ )
    sum_bc += magnetization[j][3];

mean_bc = sum_bc/num_values;

sum_sq_dev_bc = 0.0;
for ( j=from_ts; j<=to_ts; j++ )
    sum_sq_dev_bc += (mean_bc - magnetization[j][3])*(mean_bc - magnetization[j][3]);   

std_dev_bc = sqrt(sum_sq_dev_bc/num_values);

for ( pr=5; pr>=2; pr--)
    {
    fprintf(of,"\n%-9.*f%-9.*f%-9.*f%-9.*f",
        pr,mean_bc,pr,std_dev_bc,
        pr,mean_bc-std_dev_bc,pr,mean_bc+std_dev_bc);
    }

fprintf(of,"\n");
}

//  Average the mean and the standard deviation of the internal energy 
//  and the mean of the specific heat over the last pc% of time slices.  
//  Estimate standard deviation of specific heat from std dev of mean value.
/*-------------------------------------------------*/
static void calculate_ie_sh_results(int pc, FILE *of)
{
double sum_ie, sum_std_dev_ie, sum_sh, sum_sq_dev_sh;
int j, pr;
int to_ts = num_time_slices - 1;
int from_ts = 1 + ((num_time_slices-1)*(100-pc))/100;
int num_values = to_ts-from_ts+1;

if ( from_ts > to_ts - 1 )
    return;

fprintf(of,"\nMean and standard deviation of the internal energy and specific heat"
           "\nover the last %d%% of timeslices, i.e., timeslices %d through %d.",
           pc,from_ts,to_ts);

fprintf(of,"\n%20s%35s","Internal energy per spin","Specific heat per spin");
fprintf(of,"\n%-9s%-9s%-9s%-9s%-3s%-9s%-9s%-9s%-9s",
           "mean","std.dev.","m-sd","m+sd","|",
           "mean","std.dev.","m-sd","m+sd");

sum_ie = sum_std_dev_ie = sum_sh = 0.0;
for ( j=from_ts; j<=to_ts; j++ )
    {
    sum_ie         += internal_energy[j][0];    //  internal energy
    sum_std_dev_ie += internal_energy[j][1];    //  std dev of internal energy
    sum_sh         += internal_energy[j][2];    //  specific heat
    }

mean_ie    = sum_ie/num_values;
std_dev_ie = sum_std_dev_ie/num_values;
mean_sh    = sum_sh/num_values;

sum_sq_dev_sh = 0.0;
for ( j=from_ts; j<=to_ts; j++ )
    sum_sq_dev_sh += (mean_sh - internal_energy[j][2])*(mean_sh - internal_energy[j][2]);   

std_dev_sh = sqrt(sum_sq_dev_sh/num_values);

for ( pr=5; pr>=2; pr--)
    {
    fprintf(of,"\n%-9.*f%-9.*f%-9.*f%-9.*f%-3s%-9.*f%-9.*f%-9.*f%-9.*f",
        pr,mean_ie,pr,std_dev_ie,
        pr,mean_ie-std_dev_ie,pr,mean_ie+std_dev_ie,
        "|",pr,mean_sh,pr,std_dev_sh,
        pr,mean_sh-std_dev_sh,pr,mean_sh+std_dev_sh);
    }

fprintf(of,"\n");
}