path: root/src/elser.c

/*
 * elser.c
 *
 * Elser Difference Map Algorithm
 *
 * (c) 2006-2007 Thomas White <taw27@cam.ac.uk>
 *
 *  synth2d - two-dimensional Fourier synthesis
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#define _GNU_SOURCE
#include <gtk/gtk.h>
#include <fftw3.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <assert.h>

#include "main.h"
#include "data.h"
#include "displaywindow.h"
#include "reflist.h"
#include "symmetry.h"
#include "gtk-valuegraph.h"

#define ELSER_START_PHASES  1

typedef enum {
	RELAXATION_LOCALLY_ORTHOGONAL,
	RELAXATION_FIENUP_HIO,
	RELAXATION_CHARGE_FLIPPING
} ElserRelaxationType;

typedef struct stuct_elserdialog {

	/* Dialog box bits */
	GtkWidget *window;
	GtkWidget *go;
	GtkWidget *stop;
	GtkWidget *reset;
	GtkWidget *solution;
	GtkWidget *beta_widget;
	GtkWidget *itnum;
	GtkWidget *graph;
	
	/* Thread control */
	GThread *work_thread;
	unsigned int running;
	unsigned int run_semaphore;
	GStaticMutex display_mutex;
	guint display_callback;
	
	/* Algorithm control */
	double beta;
	double gam1;
	double gam2;
	gint track_best;
	ElserRelaxationType relaxation;
	
	/* Algorithm guts */
	ReflectionList *reflections;		/* List of 'modulus constraints' */
	fftw_complex *in;			/* Holds structure factors when needed */
	fftw_complex *trans_out;		/* Array on the other end of the FFT of 'in' */
	fftw_plan plan_i2t;			/* Plan to get from 'in' to 'trans_out' */
	fftw_plan plan_t2i;			/* Plan to get from 'trans_out' to 'in', i.e. calculate structure factors */
	fftw_complex *out;			/* Current iterate */
	fftw_complex *display;			/* What's currently on the display */
	size_t array_size;			/* Size of the arrays in bytes */
	fftw_complex *f1;			/* f1 semi-iterate */
	fftw_complex *f2;			/* f2 semi-iterate */
	fftw_complex *scratch;			/* Scratchpad for calculations */
	unsigned int width;			/* Width of the array in pixels */
	unsigned int height;			/* Height of the array in pixels */
		
	/* Tracking */
	unsigned int n_iterations;		/* Number of iterations performed so far */
	double *e_evolution;			/* Array of values tracking ||e|| */
	size_t e_evolution_size;		/* Size of the e_evolution array */
	double best_norm;			/* Lowest value of ||e|| so far encountered */
	unsigned int best_norm_iteration;	/* Iteration number on which the best norm was achieved */
	fftw_complex *best_norm_grid;
	
} ElserDialog;

typedef struct struct_elserpeak {
	unsigned int x;
	unsigned int y;
} ElserPeak;

ElserDialog *elser_dialog = NULL;

static void elser_grid_multiply(fftw_complex *a, double b, size_t size) {
	size_t i;
	for ( i=0; i<size/sizeof(fftw_complex); i++ ) {
		a[i][0] *= b;
		a[i][1] *= b;
	}
}

static void elser_grid_add(fftw_complex *a, fftw_complex *b, size_t size) {
	size_t i;
	for ( i=0; i<size/sizeof(fftw_complex); i++ ) {
		a[i][0] += b[i][0];
		a[i][1] += b[i][1];
	}
}

static void elser_grid_subtract(fftw_complex *a, fftw_complex *b, size_t size) {
	size_t i;
	for ( i=0; i<size/sizeof(fftw_complex); i++ ) {
		a[i][0] -= b[i][0];
		a[i][1] -= b[i][1];
	}
}

static void elser_update_itnum(ElserDialog *elser_dialog) {

	char tmp[64];
	
	snprintf(tmp, 64, "<span weight=\"bold\">Iteration number: </span>%i", elser_dialog->n_iterations);
	gtk_label_set_markup(GTK_LABEL(elser_dialog->itnum), tmp);
	
}

/* Scale up from the (small) grid size to the (large) display size */
static void elser_do_display(fftw_complex *display, fftw_complex *grid, unsigned int grid_width, unsigned int grid_height,
								unsigned int disp_width, unsigned int disp_height) {

	double sx, sy;
	unsigned int dx, dy;
	
	sx = (double)disp_width/grid_width;
	sy = (double)disp_height/grid_height;
	
	for ( dx=0; dx<disp_width; dx++ ) {
		for ( dy=0; dy<disp_height; dy++ ) {
			unsigned int gx, gy;
			gx = dx / sx;	gy = dy / sy;
			display[dy + disp_height*dx][0] = grid[gy + grid_height*gx][0];
			display[dy + disp_height*dx][1] = grid[gy + grid_height*gx][1];
		}
	}

}

static gboolean elser_update_display(gpointer data) {

	ElserDialog *elser_dialog = data;

	displaywindow_switchview();
	elser_update_itnum(elser_dialog);
	
	gtk_value_graph_set_data(GTK_VALUE_GRAPH(elser_dialog->graph), elser_dialog->e_evolution, elser_dialog->n_iterations);

	elser_dialog->display_callback = 0;
	g_static_mutex_unlock(&elser_dialog->display_mutex);
	
	return FALSE;
	
}

#if 0
/* Constrain the origin to prevent wandering in a conditional projection */
static void elser_projection_origin(fftw_complex *f, ElserDialog *elser_dialog) {

	int width, height;
	unsigned int x, y;
	double re, im, ph_offs;
	
	width = elser_dialog->width;
	height = elser_dialog->height;
	
	re = f[0][0];
	im = f[0][1];
	ph_offs = atan2(im, re);
	printf("ph_offs = %f deg\n", (ph_offs/M_PI)*180);
	
	for ( y=0; y<height; y++ ) {
		for ( x=0; x<width; x++ ) {
		
			double re, im, am, ph;
			
			re = f[y + height*x][0];
			im = f[y + height*x][1];
			am = sqrt(re*re + im*im);
			ph = atan2(im, re);
			ph -= ph_offs;
			
			f[y + height*x][0] = am * cos(ph);
			f[y + height*x][1] = am * sin(ph);
			
		}
	}

}
#endif

/* Find N atom supports of size 1+2ssx,1+2ssy in array f */
static ElserPeak *elser_projection_atomicity_findatoms(fftw_complex *f, unsigned int width, unsigned int height, unsigned int N,
								unsigned int ssx, unsigned int ssy) {

	ElserPeak *peaks;
	unsigned int i, x, y;
	unsigned int positivity = 1;
	
	peaks = malloc(sizeof(ElserPeak)*N);

	for ( i=0; i<N; i++ ) {
	
		unsigned int max_x, max_y;
		double max;
	
		max = 0; max_x = 0; max_y = 0;
		for ( y=0; y<height; y++ ) {
			for ( x=0; x<width; x++ ) {
				double re, im, am;
				re = f[y + height*x][0];
				im = f[y + height*x][1];
				am = sqrt(re*re + im*im);
				if ( fabs(am) > fabs(max) ) {
					unsigned int p, allowed;
					allowed = 1;
					for ( p=0; p<i; p++ ) {
						if ( (abs(x - peaks[p].x) < ssx) || (abs(y - peaks[p].y) < ssy) ) allowed = 0;
					}
					if ( positivity && (re < 0) ) allowed = 0;
					if ( allowed ) {
						max = am;
						max_x = x; max_y = y;
					}
				}
			}
		}
		
		peaks[i].x = max_x;
		peaks[i].y = max_y;
		
	}
	
	/* Fixed support for SnOx data */
//	peaks[0].x = 0;		peaks[0].y = 0;
//	peaks[1].x = 0;		peaks[1].y = 20;
//	peaks[2].x = 21;	peaks[2].y = 0;
//	peaks[3].x = 21;	peaks[3].y = 20;
//	peaks[4].x = 42;	peaks[4].y = 0;
//	peaks[5].x = 42;	peaks[5].y = 20;
	
	return peaks;

}

/* Find "atoms" in the grid and eliminate everything else */
static void elser_projection_atomicity(fftw_complex *f, ElserDialog *elser_dialog) {

	int width, height;
	unsigned int N, i;	/* Number of atoms sought */
	int ssx, ssy, css;
	ElserPeak *peaks;
	fftw_complex *grid_copy;
	
	width = elser_dialog->width;
	height = elser_dialog->height;
	N = 8; ssx=3; ssy=3; css=2;	/* Gives 2px radius circular atom supports */
	
	peaks = elser_projection_atomicity_findatoms(f, width, height, N, ssx, ssy);
	grid_copy = malloc(elser_dialog->array_size);
	bzero(grid_copy, elser_dialog->array_size);
	
	for ( i=0; i<N; i++ ) {
		signed int x, y;
		for ( x=(signed)peaks[i].x-ssx; x<=(signed)peaks[i].x+ssx; x++ ) {
			for ( y=(signed)peaks[i].y-ssy; y<=(signed)peaks[i].y+ssy; y++ ) {
				signed int xd, yd;
				if ( ((x-peaks[i].x)*(x-peaks[i].x))+((y-peaks[i].y)*(y-peaks[i].y)) <= css*css ) {
					xd = x % width;
					yd = y % height;
					if ( xd < 0 ) xd += width;
					if ( yd < 0 ) yd += height;
					grid_copy[yd + height*xd][0] = f[yd + height*xd][0];
					grid_copy[yd + height*xd][1] = f[yd + height*xd][1];
				}
			}
		}
	}
	free(peaks);
	memcpy(f, grid_copy, elser_dialog->array_size);
	free(grid_copy);
	
//	elser_projection_origin(f, elser_dialog);

}

static void elser_projection_modulus(fftw_complex *f, ElserDialog *elser_dialog) {

	unsigned int width, height, i, x, y;
	
	width = elser_dialog->width;
	height = elser_dialog->height;
	
	/* Transform */
	memcpy(elser_dialog->trans_out, f, elser_dialog->array_size);
	fftw_execute(elser_dialog->plan_t2i);
	
	/* Scale */
	for ( x=0; x<width; x++ ) {
		for ( y=0; y<height; y++ ) {
			elser_dialog->in[y+height*x][0] /= (width*height);
			elser_dialog->in[y+height*x][1] /= (width*height);
		}
	}
	
	/* Apply constraints */
	for ( i=1; i<elser_dialog->reflections->n_reflections; i++ ) {
	
		double re, im;
		double am, ph;
		signed int h, k;
		
		h = elser_dialog->reflections->refs[i].h;
		k = elser_dialog->reflections->refs[i].k;
		
		//if ( reflist_inlist(elser_dialog->reflections, -h, -k, 0) == 0 ) {
		//	printf("Friedel error for %3i %3i\n", h, k);
		//}
		
		if ( h < 0 ) { h = width+h; }
		if ( k < 0 ) { k = height+k; }
		
		re = elser_dialog->in[k + height*h][0];
		im = elser_dialog->in[k + height*h][1];
		ph = atan2(im, re);		
		am = elser_dialog->reflections->refs[i].amplitude;
		
		elser_dialog->in[k + height*h][0] = am*cos(ph);
		elser_dialog->in[k + height*h][1] = am*sin(ph);
		
	}
	printf("F00 = %f\n", sqrt(elser_dialog->in[0][0]*elser_dialog->in[0][0] + elser_dialog->in[0][1]*elser_dialog->in[0][1]));
	symmetry_symmetrise_array(elser_dialog->in, width, height, SYMMETRY_FRIEDEL);
	
	/* Back transform */
	fftw_execute(elser_dialog->plan_i2t);
	memcpy(f, elser_dialog->trans_out, elser_dialog->array_size);

}

static double elser_calculate_norm(fftw_complex *delta, unsigned int array_size) {

	unsigned int i;
	double total, norm;
	
	total = 0;
	for ( i=0; i<array_size/sizeof(fftw_complex); i++ ) {
		double re, im;
		re = delta[i][0];
		im = delta[i][1];
		total += re*re + im*im;
	}
	norm = sqrt(total);
	
	//printf("||e||=%f\n", norm);
	
	return norm;

}

static void *elser_work(void *data) {

	ElserDialog *elser_dialog;
	unsigned int width, height;
	fftw_complex *f1;
	fftw_complex *f2;
	
	elser_dialog = data;
	elser_dialog->running = 1;
	
	width = data_width();
	height = data_height();
	
	f1 = elser_dialog->f1;
	f2 = elser_dialog->f2;
	
	while ( elser_dialog->run_semaphore ) {
	
		double norm, gam1, gam2, beta;
		
		beta = elser_dialog->beta;
		gam1 = elser_dialog->gam1;
		gam2 = elser_dialog->gam2;
		
		memcpy(f1, elser_dialog->out, elser_dialog->array_size);
		memcpy(f2, elser_dialog->out, elser_dialog->array_size);
		
		elser_projection_atomicity(f1, elser_dialog);
		elser_grid_multiply(f1, 1+gam1, elser_dialog->array_size);
		memcpy(elser_dialog->scratch, elser_dialog->out, elser_dialog->array_size);
		elser_grid_multiply(elser_dialog->scratch, gam1, elser_dialog->array_size);
		elser_grid_subtract(f1, elser_dialog->scratch, elser_dialog->array_size);		/* f1 = (1+gam1)P1(Xn) - gam1*Xn */
		
		elser_projection_modulus(f2, elser_dialog);
		elser_grid_multiply(f2, 1+gam2, elser_dialog->array_size);
		memcpy(elser_dialog->scratch, elser_dialog->out, elser_dialog->array_size);
		elser_grid_multiply(elser_dialog->scratch, gam2, elser_dialog->array_size);
		elser_grid_subtract(f2, elser_dialog->scratch, elser_dialog->array_size);		/* f2 = (1+gam2)P2(Xn) - gam2*Xn */
		
		elser_projection_atomicity(f2, elser_dialog);
		elser_projection_modulus(f1, elser_dialog);
		elser_grid_subtract(f2, f1, elser_dialog->array_size);					/* 'f2' now contains delta = P1(f2) - P2(f1) */
		
		norm = elser_calculate_norm(f2, elser_dialog->array_size);
		elser_dialog->e_evolution[elser_dialog->n_iterations] = log(norm);
		if ( (elser_dialog->n_iterations % 100) == 99 ) {
			elser_dialog->e_evolution = realloc(elser_dialog->e_evolution, elser_dialog->e_evolution_size + 100*sizeof(double));
			elser_dialog->e_evolution_size += 100*sizeof(double);
		}
		
		elser_grid_multiply(f2, beta, elser_dialog->array_size);			/* 'f2' now contains beta(P1(f2)-P2(f1)) */
		elser_grid_add(elser_dialog->out, f2, elser_dialog->array_size);			/* Xn+1 = Xn + delta */
		
		if ( norm < elser_dialog->best_norm ) {
			elser_dialog->best_norm = norm;
			elser_dialog->best_norm_iteration = elser_dialog->n_iterations+1;
			memcpy(elser_dialog->best_norm_grid, elser_dialog->out, elser_dialog->array_size);
		}
		
		/* Tell the main thread to update the display */
		g_static_mutex_lock(&elser_dialog->display_mutex);
		elser_dialog->n_iterations++;
		if ( elser_dialog->track_best ) {
			elser_do_display(elser_dialog->display, elser_dialog->best_norm_grid, elser_dialog->width, elser_dialog->height, data_width(), data_height());
		} else {
			elser_do_display(elser_dialog->display, elser_dialog->out, elser_dialog->width, elser_dialog->height, data_width(), data_height());
		}
		elser_dialog->display_callback = g_idle_add(elser_update_display, elser_dialog);
		
	}

	elser_dialog->running = 0;
	return NULL;

}

static gint elser_control_reset(GtkWidget *reset, ElserDialog *elser_dialog) {

	unsigned int i, width, height;
	
	/* Create structures */
	width = 2*main_max_h() + 1;
	height = 2*main_max_k() + 1;
	printf("EL: Grid size is %i x %i\n", width, height);
	elser_dialog->array_size = width*height*sizeof(fftw_complex);
	if ( !elser_dialog->in ) elser_dialog->in = fftw_malloc(elser_dialog->array_size);
	if ( !elser_dialog->trans_out ) elser_dialog->trans_out = fftw_malloc(elser_dialog->array_size);
	if ( !elser_dialog->out ) elser_dialog->out = fftw_malloc(elser_dialog->array_size);
	if ( !elser_dialog->f1 ) elser_dialog->f1 = fftw_malloc(elser_dialog->array_size);
	if ( !elser_dialog->f2 ) elser_dialog->f2 = fftw_malloc(elser_dialog->array_size);
	if ( !elser_dialog->scratch ) elser_dialog->scratch = fftw_malloc(elser_dialog->array_size);
	if ( !elser_dialog->best_norm_grid ) elser_dialog->best_norm_grid = fftw_malloc(elser_dialog->array_size);
	elser_dialog->width = width;
	elser_dialog->height = height;
	elser_dialog->reflections = reflist_copy(main_reflist());
	if ( !elser_dialog->plan_i2t )
		elser_dialog->plan_i2t = fftw_plan_dft_2d(width, height, elser_dialog->in, elser_dialog->trans_out, FFTW_BACKWARD, FFTW_MEASURE | FFTW_PRESERVE_INPUT);
	if ( !elser_dialog->plan_t2i )
		elser_dialog->plan_t2i = fftw_plan_dft_2d(width, height, elser_dialog->trans_out, elser_dialog->in, FFTW_FORWARD, FFTW_MEASURE | FFTW_PRESERVE_INPUT);
	
	/* Set the initial iterate */
	for ( i=0; i<elser_dialog->array_size/sizeof(fftw_complex); i++ ) {
		elser_dialog->out[i][0] = random() - random();
		elser_dialog->out[i][1] = random() - random();
	}
#if ELSER_START_PHASES
//	elser_projection_modulus(elser_dialog->out, elser_dialog);
#endif
	elser_projection_atomicity(elser_dialog->out, elser_dialog);
	
	/* The old display array is going to get freed by displaywindow_set_realspace() */
	elser_dialog->display = fftw_malloc(data_width()*data_height()*sizeof(fftw_complex));
	elser_do_display(elser_dialog->display, elser_dialog->out, elser_dialog->width, elser_dialog->height, data_width(), data_height());
	displaywindow_set_realspace(elser_dialog->display, DWR_ELSER);
	displaywindow_forceview(DWV_REALSPACE);
	displaywindow_switchview();
	elser_dialog->n_iterations = 0;
	elser_update_itnum(elser_dialog);
	elser_dialog->best_norm_iteration = 0;
	elser_dialog->best_norm = INFINITY;
	
	/* Always exists - was set up when the dialog was created */
	free(elser_dialog->e_evolution);
	elser_dialog->e_evolution = malloc(100*sizeof(double));
	elser_dialog->e_evolution_size = 100*sizeof(double);
	gtk_value_graph_set_data(GTK_VALUE_GRAPH(elser_dialog->graph), elser_dialog->e_evolution, 0);

	gtk_widget_set_sensitive(elser_dialog->go, TRUE);
	gtk_widget_set_sensitive(elser_dialog->stop, FALSE);
	gtk_widget_set_sensitive(elser_dialog->reset, TRUE);
	gtk_widget_set_sensitive(elser_dialog->solution, FALSE);
	
	return 0;
	
}

static gint elser_control_stop(GtkWidget *stop, ElserDialog *elser_dialog) {

	if ( !elser_dialog->running ) return 0;
	
	assert(elser_dialog->run_semaphore == 1);
	assert(elser_dialog->work_thread != NULL);
	elser_dialog->run_semaphore = 0;
	if ( elser_dialog->display_callback ) {
		g_idle_remove_by_data(elser_dialog);
		g_static_mutex_unlock(&elser_dialog->display_mutex);
	}
	g_thread_join(elser_dialog->work_thread);
	elser_dialog->work_thread = NULL;
	
	gtk_widget_set_sensitive(elser_dialog->go, TRUE);
	gtk_widget_set_sensitive(elser_dialog->stop, FALSE);
	gtk_widget_set_sensitive(elser_dialog->reset, TRUE);
	gtk_widget_set_sensitive(elser_dialog->solution, TRUE);
	
	return 0;
}

static gint elser_control_go(GtkWidget *go, ElserDialog *elser_dialog) {

	if ( elser_dialog->running ) return 0;
	
	assert(elser_dialog->run_semaphore == 0);
	elser_dialog->run_semaphore = 1;
	g_static_mutex_init(&elser_dialog->display_mutex);
	assert(elser_dialog->work_thread == NULL);
	elser_dialog->display_callback = 0;
	elser_dialog->work_thread = g_thread_create(elser_work, elser_dialog, TRUE, NULL);
	
	gtk_widget_set_sensitive(elser_dialog->go, FALSE);
	gtk_widget_set_sensitive(elser_dialog->stop, TRUE);
	gtk_widget_set_sensitive(elser_dialog->reset, FALSE);
	gtk_widget_set_sensitive(elser_dialog->solution, FALSE);
	
	displaywindow_disablephasegrabbing();
	
	return 0;
}

static gint elser_control_solution(GtkWidget *solution, ElserDialog *elser_dialog) {

	if ( elser_dialog->running ) return 0;
	
	if ( elser_dialog->track_best ) {
		memcpy(elser_dialog->out, elser_dialog->best_norm_grid, elser_dialog->array_size);
		elser_dialog->n_iterations = elser_dialog->best_norm_iteration;
		elser_update_itnum(elser_dialog);
	}
	
	memcpy(elser_dialog->f1, elser_dialog->out, elser_dialog->array_size);
	elser_projection_atomicity(elser_dialog->f1, elser_dialog);
	elser_grid_multiply(elser_dialog->f1, 1+elser_dialog->gam1, elser_dialog->array_size);
	memcpy(elser_dialog->scratch, elser_dialog->out, elser_dialog->array_size);
	elser_grid_multiply(elser_dialog->scratch, elser_dialog->gam1, elser_dialog->array_size);
	elser_grid_subtract(elser_dialog->f1, elser_dialog->scratch, elser_dialog->array_size);		/* f1 = (1+gam1)P1(Xn) - gam1*Xn */
	elser_projection_modulus(elser_dialog->f1, elser_dialog);
	memcpy(elser_dialog->out, elser_dialog->f1, elser_dialog->array_size);

	elser_do_display(elser_dialog->display, elser_dialog->out, elser_dialog->width, elser_dialog->height, data_width(), data_height());
	gtk_value_graph_set_data(GTK_VALUE_GRAPH(elser_dialog->graph), elser_dialog->e_evolution, elser_dialog->n_iterations);
	displaywindow_switchview();

	gtk_widget_set_sensitive(elser_dialog->go, FALSE);
	gtk_widget_set_sensitive(elser_dialog->stop, FALSE);
	gtk_widget_set_sensitive(elser_dialog->reset, TRUE);
	gtk_widget_set_sensitive(elser_dialog->solution, FALSE);
	
	displaywindow_enablephasegrabbing();

	return 0;

}

static gint elser_control_close(GtkWidget *window, gint response, ElserDialog *this_elserdialog) {

	elser_control_stop(this_elserdialog->stop, this_elserdialog);
	
	gtk_widget_destroy(window);
	
	if ( this_elserdialog->in ) fftw_free(this_elserdialog->in);
	if ( this_elserdialog->out ) fftw_free(this_elserdialog->out);
	if ( this_elserdialog->f1 ) fftw_free(this_elserdialog->f1);
	if ( this_elserdialog->f2 ) fftw_free(this_elserdialog->f2);
	if ( this_elserdialog->scratch ) fftw_free(this_elserdialog->scratch);
	if ( this_elserdialog->best_norm_grid ) fftw_free(this_elserdialog->best_norm_grid);
	if ( this_elserdialog->trans_out ) fftw_free(this_elserdialog->trans_out);
	if ( this_elserdialog->plan_i2t ) fftw_destroy_plan(this_elserdialog->plan_i2t);
	if ( this_elserdialog->plan_t2i ) fftw_destroy_plan(this_elserdialog->plan_t2i);
	if ( this_elserdialog->reflections ) reflist_free(this_elserdialog->reflections);
	/* Don't free this_elserdialog->display, it's backing up the display */
	
	free(this_elserdialog);
	elser_dialog = NULL;

	return 0;
	
}

static gint elser_control_beta_changed(GtkWidget *beta, ElserDialog *elser_dialog) {

	elser_dialog->beta = gtk_range_get_value(GTK_RANGE(beta));
	
	switch ( elser_dialog->relaxation ) {
		case RELAXATION_LOCALLY_ORTHOGONAL : { elser_dialog->gam1 = -1/elser_dialog->beta; elser_dialog->gam2 = 1/elser_dialog->beta; break; }
		case RELAXATION_FIENUP_HIO : { elser_dialog->gam1 = -1; elser_dialog->gam2 = 1/elser_dialog->beta; break; }
		case RELAXATION_CHARGE_FLIPPING : { elser_dialog->gam1 = 1; elser_dialog->gam2 = -1; break; }
	}
	
	return 0;
	
}

static gint elser_control_trackbest_toggled(GtkWidget *track_best, ElserDialog *elser_dialog) {
	elser_dialog->track_best = gtk_toggle_button_get_active(GTK_TOGGLE_BUTTON(track_best));
	return 0;
}

static gint elser_control_relaxation_change(GtkWidget *relaxation_combo, ElserDialog *elser_dialog) {

	switch ( gtk_combo_box_get_active(GTK_COMBO_BOX(relaxation_combo)) ) {
		case 0 : elser_dialog->relaxation = RELAXATION_LOCALLY_ORTHOGONAL; break;
		case 1 : elser_dialog->relaxation = RELAXATION_FIENUP_HIO; break;
		case 2 : elser_dialog->relaxation = RELAXATION_CHARGE_FLIPPING; break;
		default : abort();	/* Only in the event of a screw-up */
	}
	
	elser_control_beta_changed(elser_dialog->beta_widget, elser_dialog);

	return 0;

}

void elser_dialog_open() {

	GtkWidget *vbox;
	GtkWidget *hbox;
	GtkWidget *buttons;
	GtkWidget *beta_label;
	GtkWidget *beta_hbox;
	GtkWidget *relaxation_hbox;
	GtkWidget *relaxation_label;
	GtkWidget *relaxation_combo;
	GtkWidget *track_best;
	
	if ( elser_dialog ) {
		return;
	}
	elser_dialog = malloc(sizeof(ElserDialog));
	elser_dialog->reflections = NULL;
	elser_dialog->in = NULL;
	elser_dialog->out = NULL;
	elser_dialog->trans_out = NULL;
	elser_dialog->f1 = NULL;
	elser_dialog->f2 = NULL;
	elser_dialog->scratch = NULL;
	elser_dialog->plan_i2t = NULL;
	elser_dialog->plan_t2i = NULL;
	elser_dialog->best_norm_grid = NULL;
	elser_dialog->display = NULL;
	elser_dialog->running = 0;
	elser_dialog->work_thread = NULL;
	elser_dialog->run_semaphore = 0;
	elser_dialog->running = 0;
	elser_dialog->beta = 0.7;
	elser_dialog->best_norm = INFINITY;
	elser_dialog->gam1 = -1/elser_dialog->beta;
	elser_dialog->gam2 = 1/elser_dialog->beta;
	elser_dialog->track_best = FALSE;
	
	elser_dialog->window = gtk_dialog_new_with_buttons("Elser Difference Map", GTK_WINDOW(displaywindow_gtkwindow()),
		GTK_DIALOG_DESTROY_WITH_PARENT, GTK_STOCK_CLOSE, GTK_RESPONSE_CLOSE, NULL);
	
	/* Structure */
	vbox = gtk_vbox_new(FALSE, 0);
	hbox = gtk_hbox_new(TRUE, 0);
	gtk_box_pack_start(GTK_BOX(GTK_DIALOG(elser_dialog->window)->vbox), GTK_WIDGET(hbox), TRUE, TRUE, 7);
	gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), TRUE, TRUE, 5);
	
	/* Parameters */
	beta_hbox = gtk_hbox_new(FALSE, 0);
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(beta_hbox), FALSE, FALSE, 5);
	beta_label = gtk_label_new("ß:");
	gtk_misc_set_padding(GTK_MISC(beta_label), 3, 0);
	gtk_box_pack_start(GTK_BOX(beta_hbox), GTK_WIDGET(beta_label), FALSE, FALSE, 0);
	elser_dialog->beta_widget = gtk_hscale_new_with_range(0, 1, 0.01);
	gtk_scale_set_value_pos(GTK_SCALE(elser_dialog->beta_widget), GTK_POS_RIGHT);
	gtk_range_set_value(GTK_RANGE(elser_dialog->beta_widget), elser_dialog->beta);
	gtk_box_pack_start(GTK_BOX(beta_hbox), GTK_WIDGET(elser_dialog->beta_widget), TRUE, TRUE, 0);
	g_signal_connect(G_OBJECT(elser_dialog->beta_widget), "value-changed", G_CALLBACK(elser_control_beta_changed), elser_dialog);
	
	relaxation_hbox = gtk_hbox_new(FALSE, 0);
	relaxation_label = gtk_label_new("Relaxation:");
	gtk_misc_set_alignment(GTK_MISC(relaxation_label), 1, 0.5);
	gtk_box_pack_start(GTK_BOX(relaxation_hbox), relaxation_label, FALSE, FALSE, 5);
	relaxation_combo = gtk_combo_box_new_text();
	gtk_combo_box_append_text(GTK_COMBO_BOX(relaxation_combo), "Locally Orthogonal Constraints");
	gtk_combo_box_append_text(GTK_COMBO_BOX(relaxation_combo), "Fienup Hybrid Input-Output");
	gtk_combo_box_append_text(GTK_COMBO_BOX(relaxation_combo), "Charge Flipping");
	gtk_box_pack_start(GTK_BOX(relaxation_hbox), relaxation_combo, TRUE, TRUE, 5);
	gtk_box_pack_start(GTK_BOX(vbox), relaxation_hbox, FALSE, FALSE, 5);
	gtk_combo_box_set_active(GTK_COMBO_BOX(relaxation_combo), 0);
	g_signal_connect(G_OBJECT(relaxation_combo), "changed", G_CALLBACK(elser_control_relaxation_change), elser_dialog);

	track_best = gtk_check_button_new_with_label("Track best solution");
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(track_best), FALSE, FALSE, 5);
	g_signal_connect(G_OBJECT(track_best), "toggled", G_CALLBACK(elser_control_trackbest_toggled), elser_dialog);
	
	/* Progress graph */
	elser_dialog->graph = gtk_value_graph_new();
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(elser_dialog->graph), TRUE, TRUE, 5);
	elser_dialog->e_evolution = malloc(100*sizeof(double));
	elser_dialog->e_evolution_size = 100*sizeof(double);
	gtk_value_graph_set_data(GTK_VALUE_GRAPH(elser_dialog->graph), elser_dialog->e_evolution, 0);
	
	/* Iteration number */
	elser_dialog->itnum = gtk_label_new("(uninitialised)");
	gtk_label_set_markup(GTK_LABEL(elser_dialog->itnum), "<span style=\"italic\">(uninitialised)</span>");
	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(elser_dialog->itnum), FALSE, FALSE, 5);
	
	/* Control buttons */
	buttons = gtk_hbox_new(FALSE, 0);
	gtk_box_pack_end(GTK_BOX(vbox), GTK_WIDGET(buttons), FALSE, FALSE, 5);
	/* Reset */
	elser_dialog->reset = gtk_button_new_with_label("Reset");
	gtk_button_set_image(GTK_BUTTON(elser_dialog->reset), gtk_image_new_from_stock(GTK_STOCK_MEDIA_PREVIOUS, GTK_ICON_SIZE_BUTTON));
	gtk_box_pack_start(GTK_BOX(buttons), GTK_WIDGET(elser_dialog->reset), TRUE, TRUE, 5);
	g_signal_connect(G_OBJECT(elser_dialog->reset), "clicked", G_CALLBACK(elser_control_reset), elser_dialog);
	/* Stop */
	elser_dialog->stop = gtk_button_new_from_stock(GTK_STOCK_MEDIA_STOP);
	gtk_box_pack_start(GTK_BOX(buttons), GTK_WIDGET(elser_dialog->stop), TRUE, TRUE, 5);
	g_signal_connect(G_OBJECT(elser_dialog->stop), "clicked", G_CALLBACK(elser_control_stop), elser_dialog);
	/* Go */
	elser_dialog->go = gtk_button_new_with_label("Run");
	gtk_button_set_image(GTK_BUTTON(elser_dialog->go), gtk_image_new_from_stock(GTK_STOCK_MEDIA_PLAY, GTK_ICON_SIZE_BUTTON));
	gtk_box_pack_start(GTK_BOX(buttons), GTK_WIDGET(elser_dialog->go), TRUE, TRUE, 5);
	g_signal_connect(G_OBJECT(elser_dialog->go), "clicked", G_CALLBACK(elser_control_go), elser_dialog);
	/* Iterate fixed point->solution */
	elser_dialog->solution = gtk_button_new_with_label("Solution");
	gtk_button_set_image(GTK_BUTTON(elser_dialog->solution), gtk_image_new_from_stock(GTK_STOCK_MEDIA_NEXT, GTK_ICON_SIZE_BUTTON));
	gtk_box_pack_start(GTK_BOX(buttons), GTK_WIDGET(elser_dialog->solution), TRUE, TRUE, 5);
	g_signal_connect(G_OBJECT(elser_dialog->solution), "clicked", G_CALLBACK(elser_control_solution), elser_dialog);
	
	g_signal_connect(G_OBJECT(elser_dialog->window), "response", G_CALLBACK(elser_control_close), elser_dialog);
	
	gtk_widget_set_sensitive(elser_dialog->go, FALSE);
	gtk_widget_set_sensitive(elser_dialog->stop, FALSE);
	gtk_widget_set_sensitive(elser_dialog->reset, TRUE);
	gtk_widget_set_sensitive(elser_dialog->solution, FALSE);
	
	gtk_widget_show_all(elser_dialog->window);
	displaywindow_forceview(DWV_REALSPACE);
	
	displaywindow_disablephasegrabbing();

}

void elser_grab_phases(ReflectionList *reflections) {

	size_t i;
	
	if ( elser_dialog == NULL ) {
		error_report("You need to keep the Elser dialog open");
		return;
	}
	
	memcpy(elser_dialog->trans_out, elser_dialog->out, elser_dialog->array_size);
	fftw_execute(elser_dialog->plan_t2i);
	
	for ( i=0; i<reflections->n_reflections; i++ ) {
		
		signed int h, k, l;
		double re, im;
		
		h = reflections->refs[i].h;
		k = reflections->refs[i].k;
		l = 0;
		
		if ( abs(h) > elser_dialog->width/2 ) {
			printf("Index %i %i (%i) is above the Nyquist frequency!\n", h, k, reflections->refs[i].l);
			reflections->refs[i].phase_known_set = 0;
			continue;
		}
		if ( abs(k) > elser_dialog->height/2 ) {
			printf("Index %i %i (%i) is above the Nyquist frequency!\n", h, k, reflections->refs[i].l);
			reflections->refs[i].phase_known_set = 0;
			continue;
		}
		if ( h < 0 ) { h = elser_dialog->width+h; }
		if ( k < 0 ) { k = elser_dialog->height+k; }
		
		re = elser_dialog->in[k+elser_dialog->height*h][0];
		im = elser_dialog->in[k+elser_dialog->height*h][1];
		
		reflections->refs[i].phase_known = atan2(im, re);
		reflections->refs[i].phase_known_set = 1;
		
	}

}