diff options
Diffstat (limited to 'src/elser.c')
| -rw-r--r-- | src/elser.c | 811 |
1 files changed, 811 insertions, 0 deletions
diff --git a/src/elser.c b/src/elser.c new file mode 100644 index 0000000..d562f9f --- /dev/null +++ b/src/elser.c @@ -0,0 +1,811 @@ +/* + * 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; + + } + +} + |