/* Ising Model Simulation of Critical Point Phenomena
 *   Math 164 (Scientific Computing)
 *   Final Project - Spring 2003
 *
 *   Michael Vrable
 */

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>

// The actual data structures used in the simulation.
int lattice_site_count;
int **lattice_connections;
int *lattice_spins;

// Coupling constant; determines strength of spin interactions.
double coupling = 1.0;

// The normalized temperature of the system.  (kT)
double temperature = 2.0;

// Initial spin distribution: probability that the initial spin is up.
double initial_spin_up_probability = 0.0;

// Command-line option: whether to give a verbose display that includes visual
// representation of the system.  The default is summary data only.
int flag_verbose = 0;

// Size of square lattice, if used
const int row = 16, col = 64;

// Statistic: Number of iterations needed for convergence
int iterations_needed = 0;

/* Compute a random number in the range of [0.0, 1.0) and return it.
 *
 * This function should be used in most of the code elsewhere, so that it will
 * be easy to swap out one random number generator algorithm for another.  The
 * current implementation uses the BSD random() function, which gives better
 * random numbers than the rand() function.
 */
double
random_float()
{
    // We assume that generating a 32-bit random integer will give a good
    // enough random distribution.
    int rand_int = random();

    // Convert the integer to a floating point value.
    double rand_float = (float)rand_int / RAND_MAX;

    return rand_float;
}

/* Seed the random number generator with a value based on the current time or
 * other values, so that successive calls of the same program don't produce
 * identical results. */
void
seed_random()
{
    // It doesn't matter if time_t and unsigned int aren't compatible types, as
    // long as we get some unpredictable value or another out.
    unsigned int seed = (unsigned int)time(NULL);

    srandom(seed);
}

/* Computer average and standard deviation: Allow data to be added to a dataset
 * one point at a time, and then at the end, compute the average and standard
 * deviation of all the samples given. */

// Structure for storing summary statistics
struct dataset {
    int n;			// Number of data points
    double sum_x;		// Sum of values
    double sum_x2;		// Sum of squares of values
};

// Initialize a dataset structure; should be called before adding any data
// points.
void
dataset_init(struct dataset *ds)
{
    ds->n = 0;
    ds->sum_x = 0.0;
    ds->sum_x2 = 0.0;
}

// Add a data point to the set of values.
void
dataset_add(struct dataset *ds, double value)
{
    ds->n++;
    ds->sum_x += value;
    ds->sum_x2 += value * value;
}

// Compute the mean of all the values in the dataset currently.
double
dataset_mean(struct dataset *ds)
{
    return ds->sum_x / ds->n;
}

// Compute the (sample) standard deviation of all values currently stored.
double
dataset_stddev(struct dataset *ds)
{
    double mean = dataset_mean(ds);

    // Compute sum of squares of each value minus the mean
    double sum_sq_diff =
	ds->sum_x2 - 2.0 * mean * ds->sum_x + (ds->n) * (mean*mean);

    // Compute the sample (not population) standard deviation
    return sqrt(fabs(sum_sq_diff / (ds->n - 1)));
}

/* Test function: Initialize a lattice so that it stores an (m x n) grid where
 * the edges wrap around. */

// Small helper function for computing address of a cell in the list of lattice
// sites.  n and m are dimensions of grid, x and y are coordinates of desired
// cell.
static inline int
compute_address(int m, int n, int x, int y)
{
    // First, if coordinates have wrapped, correct for that.
    x = (x + m) % m;
    y = (y + n) % n;

    // Convert to linear position in list of lattice sites.  Lattices are
    // stored one row at a time.
    int location = (y * m) + x;
    assert(location >= 0 && location < lattice_site_count);

    return location;
}

void
init_lattice_2d(int m, int n)
{
    lattice_site_count = n * m;
    lattice_spins = (int *)malloc(lattice_site_count * sizeof(int));
    // Allocate all memory together: size for the array of points, plus space
    // for the sub-arrays themselves, each of which have five elements (4
    // neighbors and a terminator).
    const int max_neighbors = 4;
    lattice_connections = (int **)malloc(lattice_site_count * sizeof(int *));
    int *connection_arrays =
	(int *)malloc(lattice_site_count * (1 + max_neighbors) * sizeof(int));

    if ( lattice_spins == NULL || lattice_connections == NULL
	 || connection_arrays == NULL ) {
	fprintf(stderr, "Memory allocation failed.\n");
	abort();
    }

    // Fill in data for every lattice site.
    for ( int y = 0; y < n; y++ ) {
	for ( int x = 0; x < m; x++ ) {
	    int site = compute_address(m, n, x, y);
	    int *adjacency_list =
		&connection_arrays[(1 + max_neighbors) * site];

	    lattice_connections[site] = adjacency_list;

	    // Fill in connections in all four directions.
	    for ( int d = 0; d < 4; d++ ) {
		adjacency_list[d] = compute_address(m, n,
						    x + (d == 0) - (d == 2),
						    y + (d == 1) - (d == 3));
	    }
	    adjacency_list[4] = -1;

	    // Initialize the spin orientation randomly
	    lattice_spins[site] =
		(random_float() < initial_spin_up_probability) ? 1 : -1;
	}
    }
}

// Free memory associated with a 2D lattice created with the above function.
void
free_lattice_2d()
{
    free(lattice_spins);
    free(lattice_connections[0]);   // Base address for all adjacency lists
    free(lattice_connections);

    // Don't leave dangling pointers, which should make errors easier to
    // detect.
    lattice_site_count = 0;
    lattice_spins = NULL;
    lattice_connections = NULL;
}

/* A new topology!  Two-dimensional lattice with a single site ("super lattice
 * site") which is connected to all the others.  Gives a bit of long-range
 * interaction in the system, which should cause some new and interesting
 * behavior. */
void
init_lattice_2d_alt(int m, int n)
{
    lattice_site_count = n * m + 1;
    lattice_spins = (int *)malloc(lattice_site_count * sizeof(int));
    // Allocate all memory together: size for the array of points, plus space
    // for the sub-arrays themselves, each of which have six elements (4 normal
    // neighbors, super neighbor, and a terminator).  There will be a separate
    // adjacency array for the super lattice site, adjacent to all.
    const int max_neighbors = 5;
    lattice_connections = (int **)malloc(lattice_site_count * sizeof(int *));
    int *connection_arrays =
	(int *)malloc(lattice_site_count * (1 + max_neighbors) * sizeof(int));
    int *super_connection_array =
	(int *)malloc(lattice_site_count * sizeof(int));

    if ( lattice_spins == NULL || lattice_connections == NULL
	 || connection_arrays == NULL || super_connection_array == NULL ) {
	fprintf(stderr, "Memory allocation failed.\n");
	abort();
    }

    // Fill in data for every lattice site.
    int super_site = n * m;	// Index to "super" lattice site

    for ( int y = 0; y < n; y++ ) {
	for ( int x = 0; x < m; x++ ) {
	    int site = compute_address(m, n, x, y);
	    int *adjacency_list =
		&connection_arrays[(1 + max_neighbors) * site];

	    lattice_connections[site] = adjacency_list;

	    // Fill in connections in all four directions.
	    for ( int d = 0; d < 4; d++ ) {
		adjacency_list[d] = compute_address(m, n,
						    x + (d == 0) - (d == 2),
						    y + (d == 1) - (d == 3));
	    }
	    adjacency_list[4] = super_site;
	    adjacency_list[5] = -1;

	    // Initialize the spin orientation randomly
	    lattice_spins[site] =
		(random_float() < initial_spin_up_probability) ? 1 : -1;
	}
    }

    // Now, the super lattice site
    for ( int i = 0; i < n * m; i++ ) {
	super_connection_array[i] = i;
    }
    super_connection_array[n * m] = -1;

    lattice_connections[super_site] = super_connection_array;
    lattice_spins[super_site] =
	(random_float() < initial_spin_up_probability) ? 1 : -1;

    // DEBUG: Super lattice is only influenced by first few sites?
    super_connection_array[5] = -1;
}

// Free memory associated with a 2D lattice created with the above function.
void
free_lattice_2d_alt()
{
    free(lattice_spins);
    free(lattice_connections[lattice_site_count - 1]); // Super adjacency list
    free(lattice_connections[0]);   // Base address for other adjacency lists
    free(lattice_connections);

    // Don't leave dangling pointers, which should make errors easier to
    // detect.
    lattice_site_count = 0;
    lattice_spins = NULL;
    lattice_connections = NULL;
}

/* A third topology--particles arranged on a binary tree, so that each particle
 * has three neighbors, with the exception of the leaves and the root.  The
 * topology probably falls somewhere between one and two dimensional.  The
 * parameter l is the number of levels in the tree. */
void
init_lattice_tree(int l)
{
    assert(l > 1);

    // Number of lattice sites is (2^l - 1).
    lattice_site_count = (1 << l) - 1;
    lattice_spins = (int *)malloc(lattice_site_count * sizeof(int));

    // Allocate memory for neighbhor lists, given that we know each site has at
    // most three neighbors.
    const int max_neighbors = 3;
    lattice_connections = (int **)malloc(lattice_site_count * sizeof(int *));
    int *connection_arrays =
	(int *)malloc(lattice_site_count * (1 + max_neighbors) * sizeof(int));

    if ( lattice_spins == NULL || lattice_connections == NULL
	 || connection_arrays == NULL ) {
	fprintf(stderr, "Memory allocation failed.\n");
	abort();
    }

    // Fill in connection information.  First node (root is special).
    lattice_connections[0] = connection_arrays;
    connection_arrays[0] = 1;
    connection_arrays[1] = 2;
    connection_arrays[2] = -1;
    lattice_spins[0] = 1;

    for ( int i = 1; i < lattice_site_count; i++ ) {
	int *adjacency_list = &connection_arrays[(1 + max_neighbors) * i];
	lattice_connections[i] = adjacency_list;

	// Each site has pointers to the parent and two children, unless it is
	// on the last level, in which case only parent link.
	adjacency_list[0] = (i - 1) / 2;
	if ( i < (lattice_site_count)/2 ) {
	    adjacency_list[1] = 2 * i + 1;
	    adjacency_list[2] = 2 * i + 2;
	    adjacency_list[3] = -1;
	} else {
	    // Modification: on last level, connect nodes together in pairs
	    // with each other.
	    int pair = i + (l << 2);
	    if ( pair >= lattice_site_count )
		pair -= (l << 1);
	    adjacency_list[1] = pair;
	    adjacency_list[2] = -1;
	}

	lattice_spins[i] = 1;
    }
}

void
free_lattice_tree()
{
    free(lattice_spins);
    free(lattice_connections[0]);   // Base address for all adjacency lists
    free(lattice_connections);

    // Don't leave dangling pointers, which should make errors easier to
    // detect.
    lattice_site_count = 0;
    lattice_spins = NULL;
    lattice_connections = NULL;
}
/*** MAIN SIMULATION LOOP ***/

/* Compute the contribution to the total energy of the system from a single
 * lattice site. */
static double
energy_contribution(int n)
{
    // First, determine the sum of the spins of the neighboring sites.
    int *adjacency_list = lattice_connections[n];
    int spin_sum = 0;
    for ( int i = 0; adjacency_list[i] != -1; i++ ) {
	spin_sum += lattice_spins[adjacency_list[i]];
    }

    // Multiply by spin of current site, to get number of aligned spins minus
    // unaligned spins.
    spin_sum *= lattice_spins[n];

    // Aligned spins are good, and so give a negative energy contribution.
    // Also, factor in the coupling constant.
    return -coupling * spin_sum;
}

/* Run an update on a lattice site.  This computes the probability for the site
 * to switch, and updates it if it does. */
static void
metropolis_update(int n)
{
    // Determine the change in energy if the spin flipped.
    double E = -2.0 * energy_contribution(n);

    if ( E < 0.0 || random_float() < exp(-E / temperature) ) {
	lattice_spins[n] = -lattice_spins[n];
    }
}

/* Update one random lattice point. */
static void
metropolis_random_update()
{
    // Determine cell to update
    int n = (int)(random_float() * lattice_site_count);
    
    assert(n >= 0 && n < lattice_site_count);
    metropolis_update(n);
}

/* Compute summary statistics for the current state of the model. */
struct summary_results {
    double magnetization;
    double energy;
};

struct summary_results
compute_summary()
{
    struct summary_results results;
    
    int spin_sum = 0;
    results.energy = 0.0;
    for ( int i = 0; i < lattice_site_count; i++ ) {
	/* Magnetization is computed simply as the average spin. */
	spin_sum += lattice_spins[i];

	/* Energy is the sum of the computed energy, later divided by number of
	 * sites. */
	results.energy += energy_contribution(i);
    }
    results.magnetization = (double)spin_sum / (double)lattice_site_count;
    results.energy = results.energy / (double)lattice_site_count;

    return results;
}

/* Print out the lattice sites as if they were located on a rectangular array.
 * The number of columns and rows to the data must be specified. */
void
print_lattice_2d(int m, int n)
{
    assert(m * n <= lattice_site_count);

    puts("\033[1;1H");
    for ( int y = 0; y < n; y++ ) {
	for ( int x = 0; x < m; x++ ) {
	    putchar((lattice_spins[y * m + x] == 1) ? '.' : ' ');
	}
	putchar('\n');
    }

    /* Also print out summary statistics at the end. */
    struct summary_results results = compute_summary();
    printf("Temperature: %.2f    \n", temperature);
    printf("Magnetization: %-.3f    \n", results.magnetization);
    printf("Energy: %-.3f    \n", results.energy);
    printf("Iterations: %d    \n", iterations_needed);
}

void
print_lattice_2d_alt(int m, int n)
{
    assert(m * n <= lattice_site_count);

    puts("\033[1;1H");
    for ( int y = 0; y < n; y++ ) {
	for ( int x = 0; x < m; x++ ) {
	    putchar((lattice_spins[y * m + x] == 1) ? '.' : ' ');
	}
	putchar('\n');
    }

    /* Also print out summary statistics at the end. */
    struct summary_results results = compute_summary();
    printf("Super Site: %c\n",
	   lattice_spins[lattice_site_count - 1] > 0 ? '+' : '-');
    printf("Temperature: %.2f\n", temperature);
    printf("Magnetization: %-.3f\n", results.magnetization);
    printf("Energy: %-.3f\n", results.energy);
}

/* Save a snapshot of the lattice to a file on disk.  A raw dump is produced,
 * with 0x00 representing a spin of -1 and 0xff a spin of +1.  The data is
 * written to the specified file. */
void
dump_lattice(const char *filename)
{
    FILE *f = fopen(filename, "w");
    if ( f != NULL ) {
	for ( int i = 0; i < lattice_site_count; i++ ) {
	    fputc(lattice_spins[i] > 0 ? 0x00 : 0xff, f);
	}
	fclose(f);
    }
}

/* Run a series of updates (as long is appears necessary). */
void
do_updates()
{
    int has_converged = 0;
    const int samples = 10;	// Number of times to sample energy change

    iterations_needed = 0;
    while ( !has_converged ) {
	// Store the energy before each iteration, and compare afterwards to
	// see if the energy went up or down.
	double old_energy = compute_summary().energy;
	int increase_count = 0;
	int decrease_count = 0;

	for ( int i = 0; i < samples; i++ ) {
	    // Perform a round of updates
	    for ( int j = 0; j < lattice_site_count; j++ )
		metropolis_random_update();

	    // Now, get energy afterwards.  Count number of times energy
	    // increased.
	    double new_energy = compute_summary().energy;
	    if ( new_energy > old_energy )
		increase_count++;
	    else if ( new_energy < old_energy )
		decrease_count++;
	}

	// We would expect increase_count to be half of samples if at
	// equilbrium.  How close were we?
	int difference = increase_count - decrease_count;
	if ( difference < 0 ) difference = -difference;
	
	// If difference < 0.25 * samples
	if ( 4 * difference < samples )
	    has_converged = 1;

	// Update statistics
	iterations_needed++;
    }
}

/* Run the simulation for a number of steps, compute results during the
 * simulation, then average and compare to determine accuracy.  Print out a
 * final result, including uncertainty. */
void
compute_and_print()
{
    int samples = 10;		// Number of samples to average

    // Number of updates to make between samples
    int approx_steps = lattice_site_count;

    // For computing summary statistics and error estimates
    struct dataset magnetization_set;
    struct dataset energy_set;
    dataset_init(&magnetization_set);
    dataset_init(&energy_set);

    // The first time this function is called, print a header line as well.
    static int first_call = 1;
    if ( first_call && !flag_verbose ) {
	printf("#Temperature\tMagnetization\t(stddev)\tEnergy\t(stddev)\tIterations\n");
	first_call = 0;
    }

    // Perform a series of updates before taking any samples, to allow the
    // system to come to equilibrium.
    do_updates();

    // Run through a loop--perform a series of updates, then sample.  Then
    // repeat, so that the specified number of samples is taken.
    for ( int i = 0; i < samples; i++ ) {
	// Do all the updates
	for ( int j = 0; j < approx_steps; j++ ) {
	    metropolis_random_update();
	}

	// Sample the current state of the system
	struct summary_results results = compute_summary();
	dataset_add(&magnetization_set, results.magnetization);
	dataset_add(&energy_set, results.energy);
    }

    // Now, print out results and precision
    if ( flag_verbose ) {
	print_lattice_2d(col, row);
    } else {
	printf("%.4f", temperature);
	printf("\t%-.6f\t%-.6f", dataset_mean(&magnetization_set),
	       dataset_stddev(&magnetization_set));
	printf("\t%-.6f\t%-.6f", dataset_mean(&energy_set),
	       dataset_stddev(&energy_set));
	printf("\t%d", iterations_needed);
	printf("\n");
	fflush(stdin);
    }
}

/* Main entry point.  Currently this just contains debugging code, and will be
 * replaced with something else later. */
int
main(int argc, char *argv[])
{
    // Process command-line options
    int c;
    while ( (c = getopt(argc, argv, "v")) != -1 ) {
	switch ( c ) {
	case 'v':
	    flag_verbose = 1;
	    break;
	case '?':
	    return 1;
	}
    }

    // Seed the random number generator so that results aren't always
    // identical.
    seed_random();

    // This should be moved elsewhere...
    // Clear the screen before printing anything else if we are giving a
    // verbose display.
    if ( flag_verbose )
	puts("\033[2J");

    // Initialize square lattice
    init_lattice_2d(col, row);

    // Initialize tree
    // init_lattice_tree(10);

    // Run through a range of temperatures and back, twice, then quit.
    const int steps = 800;
    const double max_temp = 4.0;
    for ( int time = 0; time <= steps; time++ ) {
	// Let rel_temp ramp up linearly to 1.0, then fall back to 0.0
	double rel_temp = 4.0 * time / steps;
	if ( rel_temp > 2.0 )
	    rel_temp -= 2.0;
	if ( rel_temp > 1.0 )
	    rel_temp = 2.0 - rel_temp;

	temperature = max_temp * rel_temp;

	// Either give a display of actual lattice sites, or if verbose was not
	// switched on, summary statistics.
	compute_and_print();
//	if ( flag_verbose ) {
//	    int iterations = 10;
//	    int update_steps = lattice_site_count;
//	    for ( int i = 0; i < iterations; i++ ) {
//		for ( int j = 0; j < update_steps; j++ ) {
//		    metropolis_random_update();
//		}
//		print_lattice_2d(col, row);
//	    }
//	} else {
//	    compute_and_print();
//	}

	// Produce occasional snapshots of the system for later analysis (and
	// pictures!).
	static char filename_buf[1024];
	if ( time % 4 == 0 ) {
	    snprintf(filename_buf, sizeof(filename_buf),
		     "images2/time%03d.raw", time);
	    dump_lattice(filename_buf);
	}
    }

    free_lattice_2d();
    // free_lattice_tree();

    return 0;
}