aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/likelihood.c6
-rw-r--r--src/likelihood.h4
-rw-r--r--src/process_hkl.c404
3 files changed, 130 insertions, 284 deletions
diff --git a/src/likelihood.c b/src/likelihood.c
index 417ebea3..b6a37994 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -17,12 +17,6 @@
#include "statistics.h"
#include "utils.h"
-void detwin_intensities(const double *model, double *new_pattern,
- const unsigned int *model_counts,
- ReflItemList *items)
-{
- /* Placeholder... */
-}
void scale_intensities(const double *model, double *new_pattern,
const unsigned int *model_counts,
diff --git a/src/likelihood.h b/src/likelihood.h
index ad6585eb..08956fe3 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -20,10 +20,6 @@
#include "utils.h"
-extern void detwin_intensities(const double *model, double *new_pattern,
- const unsigned int *model_counts,
- ReflItemList *items);
-
extern void scale_intensities(const double *model, double *new_pattern,
const unsigned int *model_counts,
ReflItemList *items, double f0,
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 29288154..b5b6df9a 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -56,16 +56,7 @@ static void show_help(const char *s)
" --stop-after=<n> Stop after processing n patterns. Zero means\n"
" keep going until the end of the input, and is\n"
" the default.\n"
-" -c, --compare-with=<file> Compare with reflection intensities in this file\n"
"\n"
-" -e, --output-every=<n> Analyse figures of merit after every n patterns\n"
-" Default: 1000. A value of zero means to do the\n"
-" analysis only after reading all the patterns.\n"
-" -r, --rvsq Output lists of R vs |q| (\"Luzzatti plots\")\n"
-" when analysing figures of merit.\n"
-" --detwin Correlate each new pattern with the current\n"
-" model and choose the best fitting out of the\n"
-" allowable twins.\n"
" --scale Scale each pattern for best fit with the current\n"
" model.\n"
" -y, --symmetry=<sym> Merge according to point group <sym>.\n"
@@ -73,113 +64,6 @@ static void show_help(const char *s)
}
-static void write_RvsQ(const char *name, double *ref, double *trueref,
- unsigned int *counts, double scale, UnitCell *cell)
-{
- FILE *fh;
- double smax, sbracket;
- signed int h, k, l;
-
- fh = fopen(name, "w");
-
- smax = 0.0;
- for ( h=-INDMAX; h<INDMAX; h++ ) {
- for ( k=-INDMAX; k<INDMAX; k++ ) {
- for ( l=-INDMAX; l<INDMAX; l++ ) {
- double s = 2.0*resolution(cell, h, k, l);
- if ( (lookup_count(counts, h, k, l) > 0) && (s > smax) ) {
- smax = s;
- }
- }
- }
- }
-
- for ( sbracket=0.0; sbracket<smax; sbracket+=smax/RVQDV ) {
-
- double top = 0.0;
- double bot = 0.0;
- int nhits = 0;
- int nrefl = 0;
- double R;
- double hits_per_refl;
-
- for ( h=-INDMAX; h<INDMAX; h++ ) {
- for ( k=-INDMAX; k<INDMAX; k++ ) {
- for ( l=-INDMAX; l<INDMAX; l++ ) {
-
- double s;
- int c;
-
- if ( (h==0) && (k==0) && (l==0) ) continue;
-
- c = lookup_count(counts, h, k, l);
- s = 2.0*resolution(cell, h, k, l);
- if ((s>=sbracket) && (s<sbracket+smax/RVQDV) && (c>0)) {
-
- double obs, calc, obsi;
-
- obs = lookup_intensity(ref, h, k, l);
- calc = lookup_intensity(trueref, h, k, l);
-
- obsi = obs / (double)c;
- top += pow(obsi - scale*calc, 2.0);
- bot += pow(obsi, 2.0);
- nhits += c;
- nrefl++;
-
- }
-
- }
- }
- }
-
- R = sqrt(top/bot);
- hits_per_refl = nrefl ? (double)nhits/nrefl : 0;
-
- fprintf(fh, "%8.5f %8.5f %5.2f\n", sbracket+smax/(2.0*RVQDV),
- R, hits_per_refl);
-
- }
- fclose(fh);
-}
-
-
-static void process_reflections(double *ref, unsigned int *counts,
- double *trueref, unsigned int *truecounts,
- unsigned int n_patterns,
- UnitCell *cell, int do_rvsq)
-{
- int j;
- double mean_counts;
- int ctot = 0;
- int nmeas = 0;
- double R, scale;
- FILE *fh;
-
- for ( j=0; j<LIST_SIZE; j++ ) {
- ctot += counts[j];
- if ( counts[j] > 0 ) nmeas++;
- }
- mean_counts = (double)ctot/nmeas;
-
- R = stat_r2(ref, counts, trueref, truecounts, &scale);
- STATUS("%8u: R=%5.2f%% (sf=%7.4e) mean meas/refl=%5.2f,"
- " %i reflections measured\n",
- n_patterns, R*100.0, scale, mean_counts, nmeas);
-
- if ( do_rvsq ) {
- /* Record graph of R against q for this N */
- char name[64];
- snprintf(name, 63, "R_vs_q-%u.dat", n_patterns);
- write_RvsQ(name, ref, trueref, counts, scale, cell);
- }
-
- fh = fopen("results/convergence.dat", "a");
- fprintf(fh, "%u %5.2f %5.2f\n", n_patterns, R*100.0, mean_counts);
- fclose(fh);
-}
-
-
/* Note "holo" needn't actually be a holohedral point group, if you want to try
* something strange like resolving from a low-symmetry group into an even
* lower symmetry one.
@@ -321,34 +205,132 @@ static void merge_pattern(double *model, const double *new,
}
+static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
+ int config_maxonly, int config_scale, int config_sum,
+ int config_stopafter,
+ ReflItemList *twins, const char *holo, const char *mero,
+ int n_total_patterns)
+{
+ char *rval;
+ float f0;
+ int n_nof0 = 0;
+ int f0_valid = 0;
+ int n_patterns = 0;
+ double *new_pattern = new_list_intensity();
+ ReflItemList *items = new_items();
+
+ do {
+
+ char line[1024];
+ signed int h, k, l;
+ float intensity;
+ int r;
+
+ rval = fgets(line, 1023, fh);
+ if ( (strncmp(line, "Reflections from indexing", 25) == 0)
+ || (strncmp(line, "New pattern", 11) == 0) ) {
+
+ /* Start of first pattern? */
+ if ( n_patterns == 0 ) {
+ n_patterns++;
+ continue;
+ }
+
+ /* Assume a default I0 if we don't have one by now */
+ if ( config_scale && !f0_valid ) {
+ n_nof0++;
+ f0 = 1.0;
+ }
+
+ /* Scale if requested */
+ if ( config_scale ) {
+ scale_intensities(model, new_pattern,
+ model_counts, items, f0,
+ f0_valid);
+ }
+
+ /* Start of second or later pattern */
+ merge_pattern(model, new_pattern, model_counts,
+ items, config_maxonly, config_sum, twins,
+ holo, mero);
+
+ if ( n_patterns == config_stopafter ) break;
+
+ n_patterns++;
+ clear_items(items);
+
+ progress_bar(n_patterns, n_total_patterns, "Merging");
+
+ f0_valid = 0;
+
+ }
+
+ if ( strncmp(line, "f0 = ", 5) == 0 ) {
+ r = sscanf(line, "f0 = %f", &f0);
+ if ( r != 1 ) {
+ f0 = 1.0;
+ f0_valid = 0;
+ continue;
+ }
+ f0_valid = 1;
+ }
+
+ r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
+ if ( r != 4 ) continue;
+
+ if ( (h==0) && (k==0) && (l==0) ) continue;
+
+ if ( find_item(items, h, k, l) != 0 ) {
+ ERROR("More than one measurement for %i %i %i in"
+ " pattern number %i\n", h, k, l, n_patterns);
+ }
+ set_intensity(new_pattern, h, k, l, intensity);
+ add_item(items, h, k, l);
+
+ } while ( rval != NULL );
+
+ delete_items(items);
+ free(new_pattern);
+
+ STATUS("%i patterns had no f0 valid value.\n", n_nof0);
+}
+
+
+static int count_patterns(FILE *fh)
+{
+ char *rval;
+
+ int n_total_patterns = 0;
+ do {
+ char line[1024];
+
+ rval = fgets(line, 1023, fh);
+ if ( (strncmp(line, "Reflections from indexing", 25) == 0)
+ || (strncmp(line, "New pattern", 11) == 0) ) {
+ n_total_patterns++;
+ }
+ } while ( rval != NULL );
+
+ return n_total_patterns;
+}
+
+
int main(int argc, char *argv[])
{
int c;
char *filename = NULL;
char *output = NULL;
FILE *fh;
- unsigned int n_patterns;
- double *model, *trueref = NULL;
+ double *model;
unsigned int *model_counts;
- char *rval;
UnitCell *cell;
int config_maxonly = 0;
- int config_every = 1000;
- int config_rvsq = 0;
int config_stopafter = 0;
int config_sum = 0;
- int config_detwin = 0;
int config_scale = 0;
- char *intfile = NULL;
- double *new_pattern = NULL;
unsigned int n_total_patterns;
- unsigned int *truecounts = NULL;
char *sym = NULL;
char *pdb = NULL;
- float f0;
- int f0_valid;
- int n_nof0 = 0;
- ReflItemList *items;
ReflItemList *twins;
int i;
@@ -359,11 +341,8 @@ int main(int argc, char *argv[])
{"output", 1, NULL, 'o'},
{"max-only", 0, &config_maxonly, 1},
{"output-every", 1, NULL, 'e'},
- {"rvsq", 0, NULL, 'r'},
{"stop-after", 1, NULL, 's'},
- {"compare-with", 0, NULL, 'c'},
{"sum", 0, &config_sum, 1},
- {"detwin", 0, &config_detwin, 1},
{"scale", 0, &config_scale, 1},
{"symmetry", 1, NULL, 'y'},
{"pdb", 1, NULL, 'p'},
@@ -387,23 +366,10 @@ int main(int argc, char *argv[])
output = strdup(optarg);
break;
-
- case 'r' :
- config_rvsq = 1;
- break;
-
- case 'e' :
- config_every = atoi(optarg);
- break;
-
case 's' :
config_stopafter = atoi(optarg);
break;
- case 'c' :
- intfile = strdup(optarg);
- break;
-
case 'p' :
pdb = strdup(optarg);
break;
@@ -426,15 +392,6 @@ int main(int argc, char *argv[])
return 1;
}
- if ( intfile != NULL ) {
- truecounts = new_list_count();
- STATUS("Comparing against '%s'\n", intfile);
- trueref = read_reflections(intfile, truecounts, NULL);
- free(intfile);
- } else {
- trueref = NULL;
- }
-
if ( pdb == NULL ) {
pdb = strdup("molecule.pdb");
}
@@ -445,37 +402,10 @@ int main(int argc, char *argv[])
model = new_list_intensity();
model_counts = new_list_count();
- new_pattern = new_list_intensity();
- items = new_items();
cell = load_cell_from_pdb(pdb);
free(pdb);
- if ( strcmp(filename, "-") == 0 ) {
- fh = stdin;
- } else {
- fh = fopen(filename, "r");
- }
- free(filename);
- if ( fh == NULL ) {
- ERROR("Failed to open input file\n");
- return 1;
- }
-
- /* Count the number of patterns in the file */
- n_total_patterns = 0;
- do {
- char line[1024];
-
- rval = fgets(line, 1023, fh);
- if ( (strncmp(line, "Reflections from indexing", 25) == 0)
- || (strncmp(line, "New pattern", 11) == 0) ) {
- n_total_patterns++;
- }
- } while ( rval != NULL );
- rewind(fh);
- STATUS("There are %i patterns to process\n", n_total_patterns);
-
/* Show useful symmetry information */
const char *holo = get_holohedral(sym);
int np = num_general_equivs(holo) / num_general_equivs(sym);
@@ -496,112 +426,38 @@ int main(int argc, char *argv[])
twins = NULL;
}
- n_patterns = 0;
- f0_valid = 0;
- do {
-
- char line[1024];
- signed int h, k, l;
- float intensity;
- int r;
-
- rval = fgets(line, 1023, fh);
- if ( (strncmp(line, "Reflections from indexing", 25) == 0)
- || (strncmp(line, "New pattern", 11) == 0) ) {
-
- /* Start of first pattern? */
- if ( n_patterns == 0 ) {
- n_patterns++;
- continue;
- }
-
- /* Detwin before scaling */
- if ( config_detwin ) {
- detwin_intensities(model, new_pattern,
- model_counts, items);
- }
-
- /* Assume a default I0 if we don't have one by now */
- if ( config_scale && !f0_valid ) {
- n_nof0++;
- f0 = 1.0;
- }
-
- /* Scale if requested */
- if ( config_scale ) {
- scale_intensities(model, new_pattern,
- model_counts, items, f0,
- f0_valid);
- }
-
- /* Start of second or later pattern */
- merge_pattern(model, new_pattern, model_counts,
- items, config_maxonly, config_sum, twins,
- holo, sym);
-
- if ( (trueref != NULL) && config_every
- && (n_patterns % config_every == 0) ) {
- process_reflections(model, model_counts,
- trueref, truecounts,
- n_patterns, cell,
- config_rvsq);
- }
-
- if ( n_patterns == config_stopafter ) break;
-
- n_patterns++;
- clear_items(items);
-
- progress_bar(n_patterns, n_total_patterns, "Merging");
-
- f0_valid = 0;
-
- }
-
- if ( strncmp(line, "f0 = ", 5) == 0 ) {
- r = sscanf(line, "f0 = %f", &f0);
- if ( r != 1 ) {
- f0 = 1.0;
- f0_valid = 0;
- continue;
- }
- f0_valid = 1;
- }
-
- r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
- if ( r != 4 ) continue;
-
- if ( (h==0) && (k==0) && (l==0) ) continue;
+ /* Open the data stream */
+ if ( strcmp(filename, "-") == 0 ) {
+ fh = stdin;
+ } else {
+ fh = fopen(filename, "r");
+ }
+ free(filename);
+ if ( fh == NULL ) {
+ ERROR("Failed to open input file\n");
+ return 1;
+ }
- if ( find_item(items, h, k, l) != 0 ) {
- ERROR("More than one measurement for %i %i %i in"
- " pattern number %i\n", h, k, l, n_patterns);
- }
- set_intensity(new_pattern, h, k, l, intensity);
- add_item(items, h, k, l);
+ /* Count the number of patterns in the file */
+ n_total_patterns = count_patterns(fh);
+ STATUS("There are %i patterns to process\n", n_total_patterns);
+ rewind(fh);
- } while ( rval != NULL );
+ merge_all(fh, model, model_counts,
+ config_maxonly, config_scale, config_sum, config_stopafter,
+ twins, holo, sym, n_total_patterns);
+ rewind(fh);
fclose(fh);
- if ( trueref != NULL ) {
- process_reflections(model, model_counts, trueref, truecounts,
- n_patterns, cell, config_rvsq);
- }
-
if ( output != NULL ) {
write_reflections(output, model_counts, model, NULL,
0, cell, 1);
}
- STATUS("There were %u patterns.\n", n_patterns);
- STATUS("%i had no f0 valid value.\n", n_nof0);
-
- delete_items(items);
free(sym);
free(model);
free(model_counts);
- free(new_pattern);
free(output);
free(cell);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-17 01:50:26 +0000