diff options
Diffstat (limited to 'src/cdm.c')
| -rw-r--r-- | src/cdm.c | 1202 |
1 files changed, 1202 insertions, 0 deletions
diff --git a/src/cdm.c b/src/cdm.c new file mode 100644 index 0000000..d0a1a91 --- /dev/null +++ b/src/cdm.c @@ -0,0 +1,1202 @@ +/* + * cdm.c + * + * "Conventional Direct Methods" - triplet and tangent stuff + * + * (c) 2006-2009 Thomas White <taw27@cam.ac.uk> + * + * synth2d - Two-Dimensional Crystallographic Fourier Synthesis + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdlib.h> +#include <stdio.h> +#include <gtk/gtk.h> +#include <math.h> +#include <string.h> +#include <gsl/gsl_sf_bessel.h> +#include <assert.h> +#include <stdint.h> +#include <fftw3.h> + +#include "displaywindow.h" +#include "symmetry.h" +#include "main.h" +#include "gtk-symmetry.h" +#include "data.h" + +typedef struct { + GtkWidget *emin; + GtkWidget *gmin; + GtkWidget *amin; + GtkWidget *symmetry; + GtkWidget *refine; + GtkWidget *auto_iterate; +} CDMWindow; + +static CDMContext *the_cdmctx = NULL; + +static void cdm_solutionwindow_open(CDMContext *cdm); + +/* ........................ The Maths ........................ */ + +static ReflectionList *cdm_findstrongreflections(ReflectionList *reflections_orig, Symmetry sym, double emin) { + + ReflectionList *strongest_reflections; + ReflectionList *reflections; +// unsigned int target; + unsigned int i; + + /* Initialise reflection list */ + reflections = reflist_copy(reflections_orig); + strongest_reflections = reflist_new_parent(reflections_orig); + + for ( i=0; i<reflections->n_reflections; i++ ) { + reflections->refs[i].parent_index = i; + } + + for ( i=1; i<reflections->n_reflections; i++ ) { + if ( (reflections->refs[i].amplitude >= emin) && (reflections->refs[i].h != 0) && (reflections->refs[i].k != 0) ) { + reflist_addref_parent(strongest_reflections, reflections->refs[i].h, reflections->refs[i].k, reflections->refs[i].l, + reflections->refs[i].parent_index); + } + } + +#if 0 + while ( reflections->n_reflections > 1 ) { + + unsigned int j; + double am_max = 0; + unsigned int j_max = 0; + + /* Find the next strongest reflection */ + for ( j=1; j<reflections->n_reflections; j++ ) { + + double am = reflections->refs[j].amplitude; + signed int h, k, l; + + if ( am == 0 ) continue; /* Avoid Trouble if emin is negative... */ + + h = reflections->refs[j].h; + k = reflections->refs[j].k; + l = reflections->refs[j].l; + + if ( am >= am_max ) { + am_max = am; + j_max = j; + } + + } + + /* Add reflection to the new list, with a link back to the original copy */ + reflist_addref_parent(strongest_reflections, reflections->refs[j_max].h, reflections->refs[j_max].k, reflections->refs[j_max].l, + reflections->refs[j_max].parent_index); + reflist_delref(reflections, reflections->refs[j_max].h, reflections->refs[j_max].k, reflections->refs[j_max].l); + + } + + /* Chop off everything but the first 40% of the reflections */ + target = 4*strongest_reflections->n_reflections/10; + while ( strongest_reflections->n_reflections > target ) { + reflist_delref(strongest_reflections, strongest_reflections->refs[strongest_reflections->n_reflections-1].h, + strongest_reflections->refs[strongest_reflections->n_reflections-1].k, + strongest_reflections->refs[strongest_reflections->n_reflections-1].l); + } +#endif + + return strongest_reflections; + +} + +#if 0 +static TripletList *cdm_sort_triplets(TripletList *triplet_list) { + + TripletList *sorted_triplets; + TripletList *triplet_copy; + unsigned int target; + + sorted_triplets = malloc(sizeof(TripletList)); + + triplet_copy = malloc(sizeof(TripletList)); + memcpy(triplet_copy, triplet_list, sizeof(TripletList)); + + target = 9*triplet_copy->n_triplets/10; + while ( triplet_copy->n_triplets > target ) { + + unsigned int i; + double G_max = 0; + unsigned int i_max = 0; + + /* Locate most concentrated triplet */ + for ( i=0; i<triplet_copy->n_triplets; i++ ) { + if ( triplet_copy->triplets[i].G > G_max ) { + i_max = i; + G_max = triplet_copy->triplets[i].G; + } + } + + /* Copy the best triplet into the new list */ + sorted_triplets->triplets[sorted_triplets->n_triplets].G = triplet_copy->triplets[i_max].G; + sorted_triplets->triplets[sorted_triplets->n_triplets].p.h = triplet_copy->triplets[i_max].p.h; + sorted_triplets->triplets[sorted_triplets->n_triplets].p.k = triplet_copy->triplets[i_max].p.k; + sorted_triplets->triplets[sorted_triplets->n_triplets].p.l = triplet_copy->triplets[i_max].p.l; + sorted_triplets->triplets[sorted_triplets->n_triplets].q.h = triplet_copy->triplets[i_max].q.h; + sorted_triplets->triplets[sorted_triplets->n_triplets].q.k = triplet_copy->triplets[i_max].q.k; + sorted_triplets->triplets[sorted_triplets->n_triplets].q.l = triplet_copy->triplets[i_max].q.l; + sorted_triplets->n_triplets++; + + /* Remove the best triplet from the original list */ + for ( i=i_max+1; i<triplet_copy->n_triplets; i++ ) { + triplet_copy->triplets[i-1].G = triplet_copy->triplets[i].G; + triplet_copy->triplets[i-1].p.h = triplet_copy->triplets[i].p.h; + triplet_copy->triplets[i-1].p.k = triplet_copy->triplets[i].p.k; + triplet_copy->triplets[i-1].p.l = triplet_copy->triplets[i].p.l; + triplet_copy->triplets[i-1].q.h = triplet_copy->triplets[i].q.h; + triplet_copy->triplets[i-1].q.k = triplet_copy->triplets[i].q.k; + triplet_copy->triplets[i-1].q.l = triplet_copy->triplets[i].q.l; + } + triplet_copy->n_triplets--; + + } + free(triplet_copy); + + return sorted_triplets; + +} +#endif + +static TripletList *cdm_findtriplets(ReflectionList *strongest_reflections, Symmetry sym, double emin, double gmin) { + + TripletList *triplet_list; + unsigned int p; + unsigned int q; +// TripletList *sorted_triplets; + unsigned int n_total_triplets; + + /* Initialise triplet list */ + triplet_list = malloc(sizeof(TripletList)); + triplet_list->n_triplets = 0; + + /* Iterate doubly over the reflections */ + n_total_triplets = 0; + for ( p=1; p<strongest_reflections->n_reflections; p++ ) { + for ( q=1; q<strongest_reflections->n_reflections; q++ ) { + + double E1; + double E2; + double E3; + double G; + unsigned int r; + signed int ph, pk, pl, qh, qk, ql, rh, rk, rl; + + ph = strongest_reflections->refs[p].h; + pk = strongest_reflections->refs[p].k; + pl = strongest_reflections->refs[p].l; + qh = strongest_reflections->refs[q].h; + qk = strongest_reflections->refs[q].k; + ql = strongest_reflections->refs[q].l; + + assert(ph || pk || pl); + assert(qh || qk || ql); + + /* Determine the third reflection in the triplet */ + rh = -ph-qh; + rk = -pk-qk; + if ( pl != ql ) { + printf("CDM: pl != ql: %3i %3i\n", pl, ql); + continue; + } + //assert(pl == ql); /* This failure mode should have been weeded out ages ago */ + rl = pl; + + /* See if the third reflection is available */ + if ( (r = reflist_inlist(strongest_reflections, rh, rk, rl)) ) { + + /* Trust me on this ... */ + E1 = strongest_reflections->parent->refs[strongest_reflections->refs[p].parent_index].amplitude; + E2 = strongest_reflections->parent->refs[strongest_reflections->refs[q].parent_index].amplitude; + E3 = strongest_reflections->parent->refs[strongest_reflections->refs[r].parent_index].amplitude; + + G = E1 * E2 * E3; + + /* Convincing enough? */ + if ( G > gmin ) { + triplet_list->triplets[triplet_list->n_triplets].p.h = ph; + triplet_list->triplets[triplet_list->n_triplets].p.k = pk; + triplet_list->triplets[triplet_list->n_triplets].p.l = pl; + triplet_list->triplets[triplet_list->n_triplets].q.h = qh; + triplet_list->triplets[triplet_list->n_triplets].q.k = qk; + triplet_list->triplets[triplet_list->n_triplets].q.l = ql; + triplet_list->triplets[triplet_list->n_triplets].G = G; + triplet_list->n_triplets++; + if ( triplet_list->n_triplets == MAX_TRIPLETS ) { + printf("DM: Too many triplets!\n"); + return NULL; + } + //printf("DM: Found triplet %3i,%3i %3i : %3i,%3i %3i : %3i %3i,%3i G = %f\n", + // ph, pk, pl, qh, qk, ql, rh, rk, rl, G); + } + n_total_triplets++; + + } + + } + } + + printf("DM: Have %i concentrated triplets (%0.1f%%)\n", triplet_list->n_triplets, ((double)triplet_list->n_triplets/n_total_triplets)*100); + + //sorted_triplets = cdm_sort_triplets(triplet_list); + //free(triplet_list); + //return sorted_triplets; + + return triplet_list; + +} + +static ReflectionList *cdm_definebasis(ReflectionList *reflections_orig, TripletList *triplet_list, Symmetry sym, double min_alpha) { + + signed int j, jmin; + ReflectionList *convergence; + ReflectionList *basis; + ReflectionList *reflections; + + if ( reflections_orig->n_reflections <= 1 ) { + printf("DM: No reflections!\n"); + return NULL; + } + + /* Circumvent normal creation mechanism in order to enact nastiness */ + reflections = malloc(sizeof(ReflectionList)); + if ( !reflections ) return NULL; + /* Start with the entire list of reflections and work backwards */ + memcpy(reflections, reflections_orig, sizeof(ReflectionList)); + + //printf("DM: Finding maximally connected starting set among %i reflections\n", reflections->n_reflections); + + convergence = reflist_new_parent(reflections_orig->parent); + + while ( reflections->n_reflections > 1 ) { /* Remember 0 0 0 */ + + double alpha_min = 999999; + signed int h_min = 999999; + signed int k_min = 999999; + signed int l_min = 999999; + unsigned int i; + unsigned int i_min = 0; + unsigned int p; + ReflectionList *equivalents; + + for ( i=1; i<reflections->n_reflections; i++ ) { + + signed int h, k, l; + unsigned int t; + double alpha = 0; + + h = reflections->refs[i].h; + k = reflections->refs[i].k; + l = reflections->refs[i].l; + assert( (h != 0) || (k != 0) || (l != 0) ); + + for ( t=0; t<triplet_list->n_triplets; t++ ) { + + signed int ph = triplet_list->triplets[t].p.h; + signed int pk = triplet_list->triplets[t].p.k; + + /* Does this triplet concern this reflection? */ + if ( (ph == h) && (pk == k) ) { + + signed int qh = triplet_list->triplets[t].q.h; + signed int qk = triplet_list->triplets[t].q.k; + signed int rh = -ph-qh; + signed int rk = -pk-qk; + double G = triplet_list->triplets[t].G; + + /* Are the other reflections in the triplet still in the list? */ + if ( reflist_inlist_2d(reflections, qh, qk) && reflist_inlist_2d(reflections, rh, rk) ) { + alpha += G * (gsl_sf_bessel_I1(G) / gsl_sf_bessel_I0(G)); + //printf("%3i %3i %3i alpha=%f * (%f / %f) = %f\n", h, k, l, G, gsl_sf_bessel_I1(G), gsl_sf_bessel_I0(G), alpha); + } + + } + + } + + if ( alpha < alpha_min ) { + alpha_min = alpha; + h_min = h; + k_min = k; + l_min = l; + i_min = reflections->refs[i].parent_index; + } + + } + + if ( alpha_min == 99999 ) { + printf("DM: Convergence procedure isn't working sensibly: alpha is too high for all reflections.\n"); + printf("DM: Have you normalised properly?\n"); + return NULL; + } + + /* Remove the reflection h_min,k_min from the basis list */ + reflist_delref(reflections, h_min, k_min, l_min); + //printf("DM: %3i %3i %3i removed from convergence list, alpha = %f, %i reflections remaining\n", + // h_min, k_min, l_min, alpha_min, reflections->n_reflections); + + /* Remove all equivalent reflections from the basis list */ + equivalents = symmetry_generate_equivalent_reflections(sym, h_min, k_min, l_min); + + for ( p=1; p<equivalents->n_reflections; p++ ) { + + signed int h, k, l; + + h = equivalents->refs[p].h; + k = equivalents->refs[p].k; + l = equivalents->refs[p].l; + + reflist_delref(reflections, h, k, l); + //printf("DM: Removing symmetry equivalent %3i %3i %3i\n", h, k, l); + + } + + reflist_free(equivalents); + + /* Add the first reflections (but not the equivalents) */ + reflist_addref_alpha_parent(convergence, h_min, k_min, l_min, alpha_min, i_min); + + } + + /* "Convergence" is in increasing order of connectivity. Reverse the order, examining + the last 50 reflections on the convergence list */ + basis = reflist_new_parent(reflections_orig->parent); + jmin = (convergence->n_reflections)-13; + if ( jmin < 1 ) jmin = 1; + for ( j=(convergence->n_reflections)-1; j>=jmin; j-- ) { + if ( (convergence->refs[j].alpha < min_alpha) || ( basis->n_reflections<10) ) { + unsigned int parent; + //printf("DM: Adding %3i %3i %3i to starting set\n", convergence->refs[j].h, convergence->refs[j].k, convergence->refs[j].l); + parent = convergence->refs[j].parent_index; + reflist_addref_parent(basis, convergence->refs[j].h, convergence->refs[j].k, convergence->refs[j].l, parent); + } + } + + return basis; + +} + +static int cdm_phaserelationships(CDMContext *cdm) { + + if ( cdm->strongest_reflections == NULL ) { + /* Get the strongest reflections */ + cdm->strongest_reflections = cdm_findstrongreflections(cdm->reflections, cdm->sym, cdm->emin); + if ( !cdm->strongest_reflections ) { + printf("DM: Error while searching for strongest reflections.\n"); + return -1; + } else { + printf("DM: Have %i strong reflections (%0.1f%%).\n", cdm->strongest_reflections->n_reflections-1, + (((double)cdm->strongest_reflections->n_reflections-1)/(cdm->reflections->n_reflections-1))*100); + if ( cdm->strongest_reflections->n_reflections == 1 ) { + printf("DM: No strong reflections found!\n"); + return -1; + } + } + } else { + printf("DM: Already have strongest reflections.\n"); + } + + if ( cdm->triplet_list == NULL ) { + /* Find the triplets */ + cdm->triplet_list = cdm_findtriplets(cdm->strongest_reflections, cdm->sym, cdm->emin, cdm->gmin); + if ( !cdm->triplet_list ) { + printf("DM: Error while searching for triplets.\n"); + return -1; + } + } else { + printf("DM: Already have triplet list.\n"); + } + + if ( cdm->basis_list == NULL ) { + /* Determine an appropriate starting point */ + cdm->basis_list = cdm_definebasis(cdm->strongest_reflections, cdm->triplet_list, cdm->sym, cdm->amin); + if ( !cdm->basis_list ) { + printf("DM: Error while defining basis set.\n"); + return -1; + } + printf("DM: %i reflections in the starting set\n", cdm->basis_list->n_reflections-1); + } else { + printf("DM: Already have basis set.\n"); + } + + return 0; + +} + +static unsigned int cdm_phasebasis(ReflectionList *basis_list, ReflectionList *reflections, Symmetry sym) { + + unsigned int i; + unsigned int n_assigned = 0; + + /* Unset all phases */ + for ( i=0; i<reflections->n_reflections; i++ ) { + reflections->refs[i].phase_calc_set = 0; + } + + /* Assign random phases to the basis set */ + for ( i=1; i<basis_list->n_reflections; i++ ) { + + ReflectionList *equivalents; + unsigned int p, j; + signed int h, k, l; + + j = basis_list->refs[i].parent_index; + h = reflections->refs[j].h; + k = reflections->refs[j].k; + l = reflections->refs[j].l; + + reflections->refs[j].phase_calc = (((double)random()/RAND_MAX) * (M_PI)) - (((double)random()/RAND_MAX) * (M_PI)); + reflections->refs[j].phase_calc_set = 1; + n_assigned++; + + /* Centricity */ + symmetry_centricity(reflections, j, sym, SYMFLAG_PHASES_CALC); + //printf("Set phase of %3i %3i %3i to %f\n", h, k, l, reflections->refs[j].phase_calc); + + /* Assign all symmetry equivalents as well */ + equivalents = symmetry_generate_equivalent_reflections(sym, h, k, l); + + for ( p=1; p<equivalents->n_reflections; p++ ) { + + unsigned int q; + + h = equivalents->refs[p].h; + k = equivalents->refs[p].k; + l = equivalents->refs[p].l; + + if ( (q = reflist_inlist(reflections, h, k, l)) ) { + + double equivalent_phase; + + equivalent_phase = (reflections->refs[j].phase_calc + equivalents->refs[p].delta_theta) * equivalents->refs[p].multiplier; + + reflections->refs[q].phase_calc = equivalent_phase; + reflections->refs[q].phase_calc_set = 1; + symmetry_centricity(reflections, q, sym, SYMFLAG_PHASES_CALC); + + //printf("DM: Setting equivalent reflection %3i %3i %3i to %f\n", h, k, l, equivalent_phase); + n_assigned++; + + } else { + //printf("DM: Need to set phase of %3i %3i %3i\n", h, k, l); + } + + } + + reflist_free(equivalents); + + } + + return n_assigned; + +} + +static unsigned int cdm_assign_code(ReflectionList *reflections, ReflectionList *basis_list, Symmetry sym, unsigned int phasing_code) { + + unsigned int i; + unsigned int n_assigned = 0; + + /* Unset all phases */ + for ( i=0; i<reflections->n_reflections; i++ ) { + reflections->refs[i].phase_calc_set = 0; + } + + /* Assign phases to the basis set */ + for ( i=1; i<basis_list->n_reflections; i++ ) { + + ReflectionList *equivalents; + unsigned int p, j; + signed int h, k, l; + + j = basis_list->refs[i].parent_index; + h = reflections->refs[j].h; + k = reflections->refs[j].k; + l = reflections->refs[j].l; + + if ( phasing_code & (1<<(i-1)) ) { + reflections->refs[j].phase_calc = M_PI; + } else { + reflections->refs[j].phase_calc = 0; + } + reflections->refs[j].phase_calc_set = 1; + n_assigned++; + + /* Centricity */ + symmetry_centricity(reflections, j, sym, SYMFLAG_PHASES_CALC); + //printf("Set phase of %3i %3i %3i to %f\n", h, k, l, reflections->refs[j].phase_calc); + + /* Assign all symmetry equivalents as well */ + equivalents = symmetry_generate_equivalent_reflections(sym, h, k, l); + + for ( p=1; p<equivalents->n_reflections; p++ ) { + + unsigned int q; + + h = equivalents->refs[p].h; + k = equivalents->refs[p].k; + l = equivalents->refs[p].l; + + if ( (q = reflist_inlist(reflections, h, k, l)) ) { + + double equivalent_phase; + + equivalent_phase = (reflections->refs[j].phase_calc + equivalents->refs[p].delta_theta) * equivalents->refs[p].multiplier; + + reflections->refs[q].phase_calc = equivalent_phase; + reflections->refs[q].phase_calc_set = 1; + symmetry_centricity(reflections, q, sym, SYMFLAG_PHASES_CALC); + + //printf("DM: Setting equivalent reflection %3i %3i %3i to %f\n", h, k, l, equivalent_phase); + n_assigned++; + + } else { + //printf("DM: Need to set phase of %3i %3i %3i\n", h, k, l); + } + + } + + reflist_free(equivalents); + + } + + return n_assigned; + +} + +static unsigned int cdm_expand(CDMContext *cdm) { + + unsigned int total_assigned = 0; + unsigned int total_refined = 0; + unsigned int n_assigned, n_refined; + ReflectionList *reflections = cdm->reflections; + ReflectionList *strongest_reflections = cdm->strongest_reflections; + TripletList *triplet_list = cdm->triplet_list; + Symmetry sym = cdm->sym; + double min_alpha = cdm->amin; + unsigned int refinement = cdm->refine; + + /* Expand as far as possible */ + do { + + unsigned int i; + n_assigned = 0; + n_refined = 0; + + for ( i=0; i<strongest_reflections->n_reflections; i++ ) { + + double beta; + unsigned int t; + unsigned int n_tri = 0; + double alpha; + double sigma_sin = 0; + double sigma_cos = 0; + signed int h = strongest_reflections->refs[i].h; + signed int k = strongest_reflections->refs[i].k; + signed int l = strongest_reflections->refs[i].l; + unsigned int j = strongest_reflections->refs[i].parent_index; + unsigned int refine, assign; + + if ( !refinement && reflections->refs[j].phase_calc_set ) continue; /* Already phased? */ + + /* Calculate the most likely phase for this reflection based on the active triplets */ + for ( t=0; t<triplet_list->n_triplets; t++ ) { + + signed int ph = triplet_list->triplets[t].p.h; + signed int pk = triplet_list->triplets[t].p.k; + + /* Does this triplet concern this reflection? */ + if ( (ph == h) && (pk == k) ) { + + signed int qh = triplet_list->triplets[t].q.h; + signed int qk = triplet_list->triplets[t].q.k; + signed int rh = -ph-qh; + signed int rk = -pk-qk; + unsigned int q = reflist_inlist_2d(reflections, qh, qk); + unsigned int r = reflist_inlist_2d(reflections, rh, rk); + + if ( !q ) { + printf("DM: q: %3i %3i not in list!\n", qh, qk); + continue; + } + if ( !r ) { + printf("DM: r: %3i %3i not in list!\n", rh, rk); + continue; + } + + /* Does this triplet concern already phased reflections? */ + if ( reflections->refs[q].phase_calc_set && reflections->refs[r].phase_calc_set ) { + + double omega = reflections->refs[q].phase_calc + reflections->refs[r].phase_calc; + + sigma_sin += triplet_list->triplets[t].G * sin(omega); + sigma_cos += triplet_list->triplets[t].G * cos(omega); + + n_tri++; + + + } + + } + + } + + alpha = sqrt(sigma_sin*sigma_sin + sigma_cos*sigma_cos); + assign = 0; refine = 0; + if ( !reflections->refs[j].phase_calc_set && (alpha > min_alpha) ) assign = 1; + if ( refinement && reflections->refs[j].phase_calc_set && (alpha > reflections->refs[j].alpha) ) refine = 1; + /* Assign phase if it's reliable enough */ + if ( assign || refine ) { + + ReflectionList *equivalents; + unsigned int p; + + beta = atan2(sigma_sin, sigma_cos); + if ( assign ) { + reflections->refs[j].phase_calc = beta; + reflections->refs[j].phase_calc_set = 1; + reflections->refs[j].alpha = alpha; + n_assigned++; + //printf("DM: Assigned %3i %3i phase %f (alpha=%f from %i triplets)\n", h, k, beta, alpha, n_tri); + } else { + reflections->refs[j].phase_calc = beta; + reflections->refs[j].alpha = alpha; + n_refined++; + //printf("DM: Refined %3i %3i phase %f (alpha=%f from %i triplets)\n", h, k, beta, alpha, n_tri); + } + + /* Centricity */ + symmetry_centricity(reflections, j, sym, SYMFLAG_PHASES_CALC); + + /* Assign all symmetry equivalents as well */ + equivalents = symmetry_generate_equivalent_reflections(sym, h, k, l); + + for ( p=1; p<equivalents->n_reflections; p++ ) { + + signed int h, k, l; + unsigned int q; + + h = equivalents->refs[p].h; + k = equivalents->refs[p].k; + l = equivalents->refs[p].l; + + if ( (q = reflist_inlist(reflections, h, k, l)) ) { + + double equivalent_phase; + + equivalent_phase = (beta + equivalents->refs[p].delta_theta) * equivalents->refs[p].multiplier; + + reflections->refs[q].phase_calc = equivalent_phase; + reflections->refs[q].phase_calc_set = 1; + reflections->refs[q].alpha = alpha; + symmetry_centricity(reflections, q, sym, SYMFLAG_PHASES_CALC); + + //printf("DM: Setting equivalent reflection %3i %3i %3i to %f\n", h, k, l, equivalent_phase); + n_assigned++; + + } else { + //printf("DM: Need to set phase of %3i %3i %3i, but it isn't in the list\n", h, k, l); + } + + } + + reflist_free(equivalents); + + } + + } + //printf("DM: This cycle: %i phases assigned, %i refined\n", n_assigned, n_refined); + total_assigned += n_assigned; + total_refined += n_refined; + + } while (n_refined+n_assigned > 0); + + cdm->n_assigned_expansion = total_assigned; + cdm->n_refined = total_refined; + return total_assigned; + +} + +static void cdm_sortsolutions(CDMContext *cdm) { + + + +} + +unsigned int cdm_systematic_expansion(CDMContext *cdm) { + + unsigned int i; + unsigned int n_solutions; + PhasingSolution *prev; + fftw_complex *cdm_in; + fftw_complex *cdm_out; + fftw_plan cdm_plan; + signed int h, k; + ReflectionList *reflections; + PhasingSolution *sol; + + if ( !(cdm->sym & SYMMETRY_CENTRE) ) { + error_report("Can only automatically iterate in a centrosymmetric planegroup"); + return 0; + } + + if ( cdm->basis_list->n_reflections > (sizeof(int)*8) ) { + error_report("Too many solutions to enumerate"); + return 0; + } + + /* Enumerate the solutions */ + n_solutions = 1<<(cdm->basis_list->n_reflections-1); + prev = NULL; + for ( i=1; i<=n_solutions; i++ ) { + + PhasingSolution *new; + + new = malloc(sizeof(PhasingSolution)); + new->n = i; + new->phasing_code = i-1; /* phasing_code might not always equal n */ + + new->next = NULL; + if ( i == 1 ) { + cdm->solutions = new; + } else { + prev->next = new; + } + prev = new; + + } + printf("DM: Enumerated %i solutions\n", n_solutions); + + /* Prepare */ + cdm_in = fftw_malloc(data_width()*data_height()*sizeof(fftw_complex)); + cdm_out = fftw_malloc(data_width()*data_height()*sizeof(fftw_complex)); + cdm_plan = fftw_plan_dft_2d(data_width(), data_height(), cdm_in, cdm_out, + FFTW_BACKWARD, FFTW_MEASURE | FFTW_PRESERVE_INPUT); + + for ( h=0; h<data_width(); h++ ) { + for ( k=0; k<data_height(); k++ ) { + cdm_in[k+data_height()*h][0] = 0; + cdm_in[k+data_height()*h][1] = 0; + } + } + + /* Step through each basis set in turn, and expand phases */ + reflections = cdm->reflections; + sol = cdm->solutions; + do { + + unsigned int basis_assigned, expnd_assigned; + int x, y; + double max, unflatness, entropy; + + printf("%7i, phasing code %06X", sol->n, sol->phasing_code); + basis_assigned = cdm_assign_code(reflections, cdm->basis_list, cdm->sym, sol->phasing_code); + expnd_assigned = cdm_expand(cdm); + + /* Transfer to Fourier array */ + for ( i=0; i<reflections->n_reflections; i++ ) { + if ( reflections->refs[i].phase_calc_set ) { + h = reflections->refs[i].h; + k = reflections->refs[i].k; + if ( h < 0 ) { h = data_width()+h; } + if ( k < 0 ) { k = data_height()+k; } + cdm_in[k + data_height()*h][0] = reflections->refs[i].amplitude * cos(reflections->refs[i].phase_calc); + cdm_in[k + data_height()*h][1] = reflections->refs[i].amplitude * sin(reflections->refs[i].phase_calc); + } + } + + fftw_execute(cdm_plan); + displaywindow_switchview(); + while ( gtk_events_pending() ) gtk_main_iteration(); + + /* Calculate "unflatness" and entropy of map, and record */ + max = 0; + for ( y=0; y<data_height(); y++ ) { + for ( x=0; x<data_width(); x++ ) { + double re = cdm_out[y+data_height()*x][0]; + double im = cdm_out[y+data_height()*x][1]; + double am = sqrt(re*re + im*im); + if ( am > max ) max = am; + } + } + unflatness = 0; entropy = 0; + for ( y=0; y<data_height(); y++ ) { + for ( x=0; x<data_width(); x++ ) { + double re = cdm_out[y+data_height()*x][0]; + double im = cdm_out[y+data_height()*x][1]; + double am = sqrt(re*re + im*im); + unflatness += (am/max) * (am/max) * (am/max) * (am/max); + entropy += (am/max) * log(am/max); + } + } + unflatness /= (data_height() * data_width()); + entropy = -entropy; + sol->unflatness = unflatness; + sol->entropy = entropy; + printf(" : %i + %i (+%i) phased, unflatness=%.5f, entropy=%.5f\n", basis_assigned, expnd_assigned, cdm->n_refined, unflatness, entropy); + + sol = sol->next; + + } while ( sol ); + + /* Finalise */ + fftw_destroy_plan(cdm_plan); + fftw_free(cdm_in); + fftw_free(cdm_out); + + cdm_sortsolutions(cdm); + cdm_solutionwindow_open(cdm); + + return 0; + +} + +unsigned int cdm_tangentexpansion(CDMContext *cdm) { + + int res; + + /* Set up the phase relationships */ + res = cdm_phaserelationships(cdm); + if ( res ) return res; + + if ( !cdm->auto_iterate ) { + + unsigned int basis_assigned, expnd_assigned; + + /* Phase the initial reflections */ + basis_assigned = cdm_phasebasis(cdm->basis_list, cdm->reflections, cdm->sym); + + /* Expand */ + expnd_assigned = cdm_expand(cdm); + + /* Report */ + printf("DM: %i + %i (+%i) phases assigned\n", basis_assigned, expnd_assigned, cdm->n_refined); + + return basis_assigned + expnd_assigned; + + } else { + + return cdm_systematic_expansion(cdm); + + } + +} + +/* ..................... The Other Stuff ..................... */ + +static int cdm_window_open = 0; + +static void cdm_execute(CDMWindow *cdm) { + + char *string; + int total_assigned; + Symmetry sym; + float emin, gmin, amin; + const char *str; + CDMContext *cdmctx; + + sym = gtk_symmetry_get_symmetry(GTK_SYMMETRY(cdm->symmetry)); + + str = gtk_entry_get_text(GTK_ENTRY(cdm->emin)); + if ( !sscanf(str, "%f", &emin) ) { + error_report("Plase specify minimum |E|"); + return; + } + str = gtk_entry_get_text(GTK_ENTRY(cdm->gmin)); + if ( !sscanf(str, "%f", &gmin) ) { + error_report("Plase specify minimum G"); + return; + } + str = gtk_entry_get_text(GTK_ENTRY(cdm->amin)); + if ( !sscanf(str, "%f", &amin) ) { + error_report("Plase specify minimum α"); + return; + } + + if ( the_cdmctx == NULL ) { + printf("DM: Creating new context.\n"); + the_cdmctx = malloc(sizeof(CDMContext)); + if ( the_cdmctx == NULL ) { + error_report("Couldn't create CDM control structure"); + return; + } + the_cdmctx->auto_iterate = gtk_toggle_button_get_active(GTK_TOGGLE_BUTTON(cdm->auto_iterate)); + the_cdmctx->refine = gtk_toggle_button_get_active(GTK_TOGGLE_BUTTON(cdm->refine)); + the_cdmctx->sym = sym; + the_cdmctx->emin = emin; + the_cdmctx->gmin = gmin; + the_cdmctx->amin = amin; + the_cdmctx->strongest_reflections = NULL; + the_cdmctx->triplet_list = NULL; + the_cdmctx->basis_list = NULL; + } else { + int restart = 0; + if ( the_cdmctx->sym != sym ) { + restart = 1; + } + if ( the_cdmctx->emin != emin ) { + restart = 1; + } + if ( the_cdmctx->gmin != gmin ) { + restart = 1; + } + /* The others do not require a re-calculation */ + if ( restart ) { + printf("DM: Conditions have changed. Recreating context...\n"); + /* Tidy up */ + free(the_cdmctx->triplet_list); + free(the_cdmctx->strongest_reflections); + free(the_cdmctx->basis_list); + free(the_cdmctx); + the_cdmctx = malloc(sizeof(CDMContext)); + if ( the_cdmctx == NULL ) { + error_report("Couldn't create CDM control structure"); + return; + } + the_cdmctx->auto_iterate = gtk_toggle_button_get_active(GTK_TOGGLE_BUTTON(cdm->auto_iterate)); + the_cdmctx->refine = gtk_toggle_button_get_active(GTK_TOGGLE_BUTTON(cdm->refine)); + the_cdmctx->sym = sym; + the_cdmctx->emin = emin; + the_cdmctx->gmin = gmin; + the_cdmctx->amin = amin; + the_cdmctx->strongest_reflections = NULL; + the_cdmctx->triplet_list = NULL; + the_cdmctx->basis_list = NULL; + } + + } + cdmctx = the_cdmctx; + + if ( !cdmctx->auto_iterate ) { + + /* Do one expansion based on a random starting set, display the result */ + total_assigned = main_cdm_tangentexpansion(cdmctx); + displaywindow_switchview(); + + string = malloc(64); + snprintf(string, 63, "%i phases assigned", total_assigned); + displaywindow_statusbar(string); + free(string); + + } else { + + /* Systematically iterate over all possible starting sets, rank by FoM */ + main_cdm_tangentexpansion(cdmctx); + + } + +} + +static gint cdm_window_response(GtkWidget *widget, gint response, CDMWindow *cdm) { + + if ( response == GTK_RESPONSE_OK ) { + cdm_execute(cdm); + } + + if ( response != GTK_RESPONSE_OK ) { + cdm_window_open = 0; + gtk_widget_destroy(widget); + } + + return 0; + +} + +void cdm_dialog_open() { + + CDMWindow *cdm; + GtkWidget *cdm_window; + GtkWidget *vbox; + GtkWidget *hbox; + GtkWidget *table; + GtkWidget *emin_label; + GtkWidget *gmin_label; + GtkWidget *amin_label; + + if ( cdm_window_open ) { + return; + } + cdm_window_open = 1; + + cdm = malloc(sizeof(CDMWindow)); + + cdm_window = gtk_dialog_new_with_buttons("Tangent Formula", GTK_WINDOW(displaywindow_gtkwindow()), + GTK_DIALOG_DESTROY_WITH_PARENT, GTK_STOCK_CLOSE, GTK_RESPONSE_CANCEL, GTK_STOCK_EXECUTE, GTK_RESPONSE_OK, NULL); + + vbox = gtk_vbox_new(FALSE, 0); + hbox = gtk_hbox_new(TRUE, 0); + gtk_box_pack_start(GTK_BOX(GTK_DIALOG(cdm_window)->vbox), GTK_WIDGET(hbox), FALSE, FALSE, 7); + gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, FALSE, 5); + + cdm->symmetry = gtk_symmetry_new(2, 2, TRUE); + gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(cdm->symmetry), TRUE, TRUE, 5); + + table = gtk_table_new(5, 2, FALSE); + gtk_table_set_row_spacings(GTK_TABLE(table), 5); + gtk_table_set_col_spacings(GTK_TABLE(table), 5); + gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(table), FALSE, FALSE, 0); + + emin_label = gtk_label_new("Minimum |E|:"); + gtk_misc_set_alignment(GTK_MISC(emin_label), 1, 0.5); + gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(emin_label), 1, 2, 1, 2); + cdm->emin = gtk_entry_new(); + gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(cdm->emin), 2, 3, 1, 2); + gtk_entry_set_text(GTK_ENTRY(cdm->emin), "0.0001"); + + gmin_label = gtk_label_new("Minimum G:"); + gtk_misc_set_alignment(GTK_MISC(gmin_label), 1, 0.5); + gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(gmin_label), 1, 2, 2, 3); + cdm->gmin = gtk_entry_new(); + gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(cdm->gmin), 2, 3, 2, 3); + gtk_entry_set_text(GTK_ENTRY(cdm->gmin), "0.00024"); + + amin_label = gtk_label_new("Minimum α:"); + gtk_misc_set_alignment(GTK_MISC(amin_label), 1, 0.5); + gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(amin_label), 1, 2, 3, 4); + cdm->amin = gtk_entry_new(); + gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(cdm->amin), 2, 3, 3, 4); + gtk_entry_set_text(GTK_ENTRY(cdm->amin), "0.0001"); + + cdm->refine = gtk_check_button_new_with_label("Refine Phases"); + gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(cdm->refine), 1, 3, 4, 5); + gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(cdm->refine), TRUE); + + cdm->auto_iterate = gtk_check_button_new_with_label("Automatically Iterate"); + gtk_table_attach_defaults(GTK_TABLE(table), GTK_WIDGET(cdm->auto_iterate), 1, 3, 5, 6); + gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(cdm->auto_iterate), FALSE); + + g_signal_connect(G_OBJECT(cdm_window), "response", G_CALLBACK(cdm_window_response), cdm); + gtk_widget_show_all(cdm_window); + +} + +enum { + SOLUTION_COLUMN_NUMBER, + SOLUTION_COLUMN_PHASECODE, + SOLUTION_COLUMN_UNFLATNESS, + SOLUTION_COLUMN_ENTROPY, + SOLUTION_COLUMNS +}; + +void cdm_display_phasing_solution(CDMContext *cdm, PhasingSolution *sol, ReflectionList *reflections) { + + cdm_assign_code(reflections, cdm->basis_list, cdm->sym, sol->phasing_code); + cdm_expand(cdm); + + displaywindow_switchview(); + +} + +static gint cdm_solutionwindow_response(GtkWidget *solution_window, gint response, gpointer data) { + gtk_widget_destroy(solution_window); + return 0; +} + +static void cdm_solutionwindow_selectionchanged(GtkTreeSelection *selection, CDMContext *cdm) { + + GtkTreePath *path; + GtkTreeViewColumn *column; + GtkTreeIter iter; + gint snum; + PhasingSolution *sol_it; + PhasingSolution *sol; + + gtk_tree_view_get_cursor(GTK_TREE_VIEW(cdm->solution_tree_view), &path, &column); + gtk_tree_model_get_iter(GTK_TREE_MODEL(cdm->solution_list_store), &iter, path); + gtk_tree_model_get(GTK_TREE_MODEL(cdm->solution_list_store), &iter, SOLUTION_COLUMN_NUMBER, &snum, -1); + + sol_it = cdm->solutions; sol = NULL; + do { + if ( sol_it->n == snum ) { + sol = sol_it; + break; + } + sol_it = sol_it->next; + } while ( sol_it ); + if ( sol ) { + printf("DM: Displaying solution %7i, phasing code %06X\n", sol->n, sol->phasing_code); + main_display_phasing_solution(cdm, sol); + } else { + printf("DM: Can't find solution %i\n", snum); + error_report("Can't find solution!\n"); + } + +} + +static void cdm_solutionwindow_open(CDMContext *cdm) { + + GtkWidget *solution_window; + GtkTreeViewColumn *column; + GtkCellRenderer *renderer; + GtkWidget *scrolledwindow; + GtkWidget *solution_tree; + GtkListStore *solution_list_store; + GtkWidget *hbox; + GtkWidget *vbox; + GtkTreeIter iter; + PhasingSolution *sol; + + solution_window = gtk_dialog_new_with_buttons("Phasing Solutions", GTK_WINDOW(displaywindow_gtkwindow()), + GTK_DIALOG_DESTROY_WITH_PARENT, GTK_STOCK_CLOSE, GTK_RESPONSE_CLOSE, NULL); + gtk_window_set_default_size(GTK_WINDOW(solution_window), 400, 500); + + /* Window layout */ + vbox = gtk_vbox_new(FALSE, 0); + hbox = gtk_hbox_new(FALSE, 0); + gtk_box_pack_start(GTK_BOX(GTK_DIALOG(solution_window)->vbox), GTK_WIDGET(hbox), TRUE, TRUE, 7); + gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), TRUE, TRUE, 5); + + /* List and Tree */ + solution_list_store = gtk_list_store_new(SOLUTION_COLUMNS, G_TYPE_INT, G_TYPE_STRING, G_TYPE_DOUBLE, G_TYPE_DOUBLE); + solution_tree = gtk_tree_view_new_with_model(GTK_TREE_MODEL(solution_list_store)); + g_object_unref(G_OBJECT(solution_list_store)); + cdm->solution_tree_view = solution_tree; + cdm->solution_list_store = solution_list_store; + + /* Set up columns */ + renderer = gtk_cell_renderer_text_new(); + column = gtk_tree_view_column_new_with_attributes("#", renderer, "text", SOLUTION_COLUMN_NUMBER, NULL); + gtk_tree_view_append_column(GTK_TREE_VIEW(solution_tree), column); + column = gtk_tree_view_column_new_with_attributes("Phase", renderer, "text", SOLUTION_COLUMN_PHASECODE, NULL); + gtk_tree_view_append_column(GTK_TREE_VIEW(solution_tree), column); + column = gtk_tree_view_column_new_with_attributes("Unflatness", renderer, "text", SOLUTION_COLUMN_UNFLATNESS, NULL); + gtk_tree_view_append_column(GTK_TREE_VIEW(solution_tree), column); + column = gtk_tree_view_column_new_with_attributes("Entropy", renderer, "text", SOLUTION_COLUMN_ENTROPY, NULL); + gtk_tree_view_append_column(GTK_TREE_VIEW(solution_tree), column); + + /* Scrollbar */ + scrolledwindow = gtk_scrolled_window_new(gtk_tree_view_get_hadjustment(GTK_TREE_VIEW(solution_tree)), + gtk_tree_view_get_vadjustment(GTK_TREE_VIEW(solution_tree))); + gtk_scrolled_window_set_policy(GTK_SCROLLED_WINDOW(scrolledwindow), GTK_POLICY_AUTOMATIC, GTK_POLICY_AUTOMATIC); + gtk_container_add(GTK_CONTAINER(scrolledwindow), GTK_WIDGET(solution_tree)); + + gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(scrolledwindow), TRUE, TRUE, 5); + + /* Add the solutions */ + sol = cdm->solutions; + do { + char tmp[32]; + gtk_list_store_append(solution_list_store, &iter); + snprintf(tmp, 31, "%06X", sol->phasing_code); + gtk_list_store_set(solution_list_store, &iter, SOLUTION_COLUMN_NUMBER, sol->n, SOLUTION_COLUMN_PHASECODE, tmp, + SOLUTION_COLUMN_UNFLATNESS, sol->unflatness, SOLUTION_COLUMN_ENTROPY, sol->entropy, -1); + sol = sol->next; + } while ( sol ); + g_signal_connect(G_OBJECT(gtk_tree_view_get_selection(GTK_TREE_VIEW(solution_tree))), "changed", + G_CALLBACK(cdm_solutionwindow_selectionchanged), cdm); + + + g_signal_connect(G_OBJECT(solution_window), "response", G_CALLBACK(cdm_solutionwindow_response), cdm); + gtk_widget_show_all(solution_window); + +} |