/*
 * 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);

}