partialator: Add -w and --operator options

3 files changed, 112 insertions, 58 deletions

diff --git a/src/partialator.c b/src/partialator.c
index 5c751732..baab84c4 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -327,7 +327,9 @@ static void show_help(const char *s)
 "      --no-free              Disable cross-validation (testing only).\n"
 "      --custom-split         List of files for custom dataset splitting.\n"
 "      --max-rel-B            Maximum allowable relative |B| factor.\n"
-"      --no-logs              Do not write extensive log files.\n");
+"      --no-logs              Do not write extensive log files.\n"
+"  -w <pg>                    Apparent point group for resolving ambiguities.\n"
+"      --operator=<op>        Indexing ambiguity operator for resolving.\n");
 }
 
 
@@ -820,6 +822,8 @@ int main(int argc, char *argv[])
 	char *outfile = NULL;
 	char *sym_str = NULL;
 	SymOpList *sym;
+	SymOpList *amb;
+	SymOpList *w_sym;
 	int nthreads = 1;
 	int i;
 	int n_iter = 10;
@@ -852,6 +856,8 @@ int main(int argc, char *argv[])
 	char *rfile = NULL;
 	RefList *reference = NULL;
 	int no_logs = 0;
+	char *w_sym_str = NULL;
+	char *operator = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -875,6 +881,7 @@ int main(int argc, char *argv[])
 		{"max-rel-B",          1, NULL,                7},
 		{"max-rel-b",          1, NULL,                7}, /* compat */
 		{"reference",          1, NULL,                8}, /* ssshhh! */
+		{"operator",           1, NULL,                9},
 
 		{"no-scale",           0, &no_scale,           1},
 		{"no-pr",              0, &no_pr,              1},
@@ -896,7 +903,7 @@ int main(int argc, char *argv[])
 	}
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:",
+	while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:w:",
 	                        longopts, NULL)) != -1)
 	{
 
@@ -955,6 +962,10 @@ int main(int argc, char *argv[])
 			pmodel_str = strdup(optarg);
 			break;
 
+			case 'w' :
+			w_sym_str = strdup(optarg);
+			break;
+
 			case 2 :
 			errno = 0;
 			min_measurements = strtod(optarg, &rval);
@@ -1005,6 +1016,10 @@ int main(int argc, char *argv[])
 			rfile = strdup(optarg);
 			break;
 
+			case 9 :
+			operator = strdup(optarg);
+			break;
+
 			case 0 :
 			break;
 
@@ -1044,6 +1059,44 @@ int main(int argc, char *argv[])
 	free(sym_str);
 	if ( sym == NULL ) return 1;
 
+	if ( (w_sym_str != NULL) && (operator != NULL) ) {
+		ERROR("Specify the apparent symmetry (-w) or the operator, "
+		      "not both.\n");
+		return 1;
+	}
+
+	if ( w_sym_str == NULL ) {
+		w_sym = NULL;
+		amb = NULL;
+	} else {
+		pointgroup_warning(w_sym_str);
+		w_sym = get_pointgroup(w_sym_str);
+		free(w_sym_str);
+		if ( w_sym == NULL ) return 1;
+		amb = get_ambiguities(w_sym, sym);
+		if ( amb == NULL ) {
+			ERROR("Couldn't find ambiguity operator.\n");
+			ERROR("Check that your values for -y and -w are "
+			      "correct.\n");
+			return 1;
+		}
+
+	}
+
+	if ( operator ) {
+		amb = parse_symmetry_operations(operator);
+		if ( amb == NULL ) return 1;
+		set_symmetry_name(amb, "Ambiguity");
+	}
+
+	if ( amb != NULL ) {
+		STATUS("Indexing ambiguity resolution enabled.  "
+		       "The ambiguity operation(s) are:\n");
+		describe_symmetry(amb);
+		/* In contrast to ambigator, partialator can deal with multiple
+		 * ambiguities at once */
+	}
+
 	if ( pmodel_str != NULL ) {
 		if ( strcmp(pmodel_str, "unity") == 0 ) {
 			pmodel = PMODEL_UNITY;
@@ -1264,7 +1317,7 @@ int main(int argc, char *argv[])
 
 		if ( !no_pr ) {
 			refine_all(crystals, n_crystals, full, nthreads, pmodel,
-			           0, i+1, no_logs);
+			           0, i+1, no_logs, sym, amb);
 		} else if ( !no_scale ) {
 			scale_all_to_reference(crystals, n_crystals, full, nthreads);
 		}
@@ -1357,7 +1410,7 @@ int main(int argc, char *argv[])
 		crystal_free(crystals[i]);
 	}
 	reflist_free(full);
-	free(sym);
+	free_symoplist(sym);
 	free(outfile);
 	free(crystals);
 	free(infile);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 1c33e211..7e5420d3 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -389,7 +389,7 @@ static int check_angle_shifts(gsl_vector *cur, const gsl_vector *initial,
 
 
 static RefList *reindex_reflections(RefList *input, SymOpList *amb,
-                                    SymOpList *sym)
+                                    SymOpList *sym, int idx)
 {
 	RefList *n;
 	Reflection *refl;
@@ -405,13 +405,13 @@ static RefList *reindex_reflections(RefList *input, SymOpList *amb,
 		Reflection *rn;
 
 		get_indices(refl, &h, &k, &l);
-		get_equiv(amb, NULL, 0, h, k, l, &h, &k, &l);
+		get_equiv(amb, NULL, idx, h, k, l, &h, &k, &l);
 		get_asymm(sym, h, k, l, &h, &k, &l);
 		rn = add_refl(n, h, k, l);
 		copy_data(rn, refl);
 
 		get_symmetric_indices(rn, &h, &k, &l);
-		get_equiv(amb, NULL, 0, h, k, l, &h, &k, &l);
+		get_equiv(amb, NULL, idx, h, k, l, &h, &k, &l);
 		set_symmetric_indices(rn, h, k, l);
 	}
 
@@ -419,7 +419,7 @@ static RefList *reindex_reflections(RefList *input, SymOpList *amb,
 }
 
 
-static void reindex_cell(UnitCell *cell, SymOpList *amb)
+static void reindex_cell(UnitCell *cell, SymOpList *amb, int idx)
 {
 	signed int h, k, l;
 	struct rvec na, nb, nc;
@@ -429,17 +429,17 @@ static void reindex_cell(UnitCell *cell, SymOpList *amb)
 	                          &bs.u, &bs.v, &bs.w,
 	                          &cs.u, &cs.v, &cs.w);
 
-	get_equiv(amb, NULL, 0, 1, 0, 0, &h, &k, &l);
+	get_equiv(amb, NULL, idx, 1, 0, 0, &h, &k, &l);
 	na.u = as.u*h + bs.u*k + cs.u*l;
 	na.v = as.v*h + bs.v*k + cs.v*l;
 	na.w = as.w*h + bs.w*k + cs.w*l;
 
-	get_equiv(amb, NULL, 0, 0, 1, 0, &h, &k, &l);
+	get_equiv(amb, NULL, idx, 0, 1, 0, &h, &k, &l);
 	nb.u = as.u*h + bs.u*k + cs.u*l;
 	nb.v = as.v*h + bs.v*k + cs.v*l;
 	nb.w = as.w*h + bs.w*k + cs.w*l;
 
-	get_equiv(amb, NULL, 0, 0, 0, 1, &h, &k, &l);
+	get_equiv(amb, NULL, idx, 0, 0, 1, &h, &k, &l);
 	nc.u = as.u*h + bs.u*k + cs.u*l;
 	nc.v = as.v*h + bs.v*k + cs.v*l;
 	nc.w = as.w*h + bs.w*k + cs.w*l;
@@ -450,45 +450,50 @@ static void reindex_cell(UnitCell *cell, SymOpList *amb)
 }
 
 
-void try_reindex(Crystal *crin, const RefList *full)
+void try_reindex(Crystal *crin, const RefList *full,
+                 SymOpList *sym, SymOpList *amb)
 {
 	RefList *list;
 	Crystal *cr;
 	UnitCell *cell;
-	SymOpList *amb;
-	SymOpList *sym;
 	double residual_original, residual_flipped;
+	int idx, n;
 
-	amb = parse_symmetry_operations("k,h,-l");
-	sym = get_pointgroup("m-3");
+	if ( sym == NULL || amb == NULL ) return;
 
 	residual_original = residual(crin, full, 0, NULL, NULL, 0);
 
 	cr = crystal_copy(crin);
 
-	cell = cell_new_from_cell(crystal_get_cell(crin));
-	if ( cell == NULL ) return;
-	reindex_cell(cell, amb);
-	crystal_set_cell(cr, cell);
+	n = num_equivs(amb, NULL);
 
-	list = reindex_reflections(crystal_get_reflections(crin), amb, sym);
-	crystal_set_reflections(cr, list);
+	for ( idx=0; idx<n; idx++ ) {
 
-	update_predictions(cr);
-	calculate_partialities(cr, PMODEL_XSPHERE);
+		cell = cell_new_from_cell(crystal_get_cell(crin));
+		if ( cell == NULL ) return;
+		reindex_cell(cell, amb, idx);
+		crystal_set_cell(cr, cell);
 
-	residual_flipped = residual(cr, full, 0, NULL, NULL, 1);
+		list = reindex_reflections(crystal_get_reflections(crin),
+		                           amb, sym, idx);
+		crystal_set_reflections(cr, list);
 
-	if ( residual_flipped < residual_original ) {
-		crystal_set_cell(crin, cell);
-		crystal_set_reflections(crin, list);
-	} else {
-		cell_free(crystal_get_cell(cr));
-		reflist_free(crystal_get_reflections(cr));
-		crystal_free(cr);
+		update_predictions(cr);
+		calculate_partialities(cr, PMODEL_XSPHERE);
+
+		residual_flipped = residual(cr, full, 0, NULL, NULL, 1);
+
+		if ( residual_flipped < residual_original ) {
+			crystal_set_cell(crin, cell);
+			crystal_set_reflections(crin, list);
+			residual_original = residual_flipped;
+		} else {
+			cell_free(crystal_get_cell(cr));
+			reflist_free(crystal_get_reflections(cr));
+		}
 	}
 
-	free_symoplist(amb);
+	crystal_free(cr);
 }
 
 
@@ -686,7 +691,8 @@ void write_gridscan(Crystal *cr, const RefList *full,
 
 static void do_pr_refine(Crystal *cr, const RefList *full,
                          PartialityModel pmodel, int verbose, int serial,
-                         int cycle, int write_logs)
+                         int cycle, int write_logs,
+                         SymOpList *sym, SymOpList *amb)
 {
 	gsl_multimin_fminimizer *min;
 	struct rf_priv priv;
@@ -697,7 +703,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
 	double residual_start, residual_free_start;
 	FILE *fh = NULL;
 
-	try_reindex(cr, full);
+	try_reindex(cr, full, sym, amb);
 
 	residual_start = residual(cr, full, 0, NULL, NULL, 0);
 	residual_free_start = residual(cr, full, 1, NULL, NULL, 0);
@@ -858,24 +864,6 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
 }
 
 
-static struct prdata pr_refine(Crystal *cr, const RefList *full,
-                               PartialityModel pmodel, int verbose, int serial,
-                               int cycle, int write_logs)
-{
-	struct prdata prdata;
-
-	prdata.refined = 0;
-
-	do_pr_refine(cr, full, pmodel, verbose, serial, cycle, write_logs);
-
-	if ( crystal_get_user_flag(cr) == 0 ) {
-		prdata.refined = 1;
-	}
-
-	return prdata;
-}
-
-
 struct refine_args
 {
 	RefList *full;
@@ -886,6 +874,8 @@ struct refine_args
 	int verbose;
 	int cycle;
 	int no_logs;
+	SymOpList *sym;
+	SymOpList *amb;
 };
 
 
@@ -906,9 +896,15 @@ static void refine_image(void *task, int id)
 	int write_logs = 0;
 
 	write_logs = !pargs->no_logs && (pargs->serial % 20 == 0);
-	pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel,
-	                          pargs->verbose, pargs->serial, pargs->cycle,
-	                          write_logs);
+	pargs->prdata.refined = 0;
+
+	do_pr_refine(cr, pargs->full, pargs->pmodel, pargs->verbose,
+	             pargs->serial, pargs->cycle, write_logs,
+	             pargs->sym, pargs->amb);
+
+	if ( crystal_get_user_flag(cr) == 0 ) {
+		pargs->prdata.refined = 1;
+	}
 }
 
 
@@ -942,7 +938,8 @@ static void done_image(void *vqargs, void *task)
 
 void refine_all(Crystal **crystals, int n_crystals,
                 RefList *full, int nthreads, PartialityModel pmodel,
-                int verbose, int cycle, int no_logs)
+                int verbose, int cycle, int no_logs,
+                SymOpList *sym, SymOpList *amb)
 {
 	struct refine_args task_defaults;
 	struct queue_args qargs;
@@ -954,6 +951,8 @@ void refine_all(Crystal **crystals, int n_crystals,
 	task_defaults.verbose = verbose;
 	task_defaults.cycle = cycle;
 	task_defaults.no_logs = no_logs;
+	task_defaults.sym = sym;
+	task_defaults.amb = amb;
 
 	qargs.task_defaults = task_defaults;
 	qargs.n_started = 0;
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 3c4ac7ef..a27127ff 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -41,6 +41,7 @@
 #include "utils.h"
 #include "crystal.h"
 #include "geometry.h"
+#include "symmetry.h"
 
 
 enum prflag
@@ -59,7 +60,8 @@ extern const char *str_prflag(enum prflag flag);
 
 extern void refine_all(Crystal **crystals, int n_crystals,
                        RefList *full, int nthreads, PartialityModel pmodel,
-                       int verbose, int cycle, int no_logs);
+                       int verbose, int cycle, int no_logs,
+                       SymOpList *sym, SymOpList *amb);
 
 extern void write_gridscan(Crystal *cr, const RefList *full,
                            int cycle, int serial);