Skip to content
GitLab
Commit b7620b55 authored by Thomas White's avatar Thomas White
Browse files

process_hkl: Generate intensity histograms

parent ab97b033
Hide whitespace changes
Inline Side-by-side
scripts/frequency
View file @ b7620b55
......@@ -2,168 +2,37 @@
use strict;
my $bins = 200;
my $nplots = 32;
my $sym = "6/mmm";
open(FH, $ARGV[0]);
my $stream = $ARGV[0];
my $line;
my %counts;
# --- Step 1: Count reflections and find the most popular ones
while ( $line = <FH> ) {
my $h;
my $k;
my $l;
chomp($line);
if ( $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)/ ) {
my $key;
$h = $1;
$k = $2;
$l = $3;
$key = $h.",".$k.",".$l;
$counts{$key}++;
}
}
my $ref;
my $i;
$i = 0;
foreach $ref ( sort {-($counts{$a} <=> $counts{$b})} keys %counts )
{
$i++;
if ( $i > $nplots ) {
delete($counts{$ref});
} else {
printf("%s %i\n", $ref, $counts{$ref});
}
}
# --- Step 2: Work out histogram bins for each reflection
seek(FH, 0, 0);
my %mins;
my %maxs;
#system("~/crystfel/src/process_hkl -i ".$stream." -o test.hkl -y ".$sym);
open(FH, "test.hkl");
my $max = 0;
my $h;
my $k;
my $l;
my $line;
while ( $line = <FH> ) {
my $key;
my $intensity;
chomp($line);
if ( $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)\s+([\d\.]+)/ ) {
my $h;
my $k;
my $l;
chomp $line;
$h = $1;
$k = $2;
$l = $3;
$intensity = $4;
$key = $h.",".$k.",".$l;
next unless exists($counts{$key});
unless ( exists($mins{$key}) ) {
$mins{$key} = $intensity;
}
unless ( exists($maxs{$key}) ) {
$maxs{$key} = $intensity;
}
if ( $intensity < $mins{$key} ) {
$mins{$key} = $intensity;
}
if ( $intensity > $maxs{$key} ) {
$maxs{$key} = $intensity;
}
}
}
# --- Step 3: Bin the counts and plot
open(GP, "| gnuplot");
printf(GP "set term postscript enhanced font \"Helvetica,20\"\n");
printf(GP "set output \"histo.ps\"\n");
printf(GP "set xtics nomirror out rotate by -60\n");
#printf(GP "set logscale x\n");
foreach $ref ( sort {-($counts{$a} <=> $counts{$b})} keys %counts )
{
my $max = $maxs{$ref};
my $min = $mins{$ref};
my $binsize = ($max+1-$min)/$bins;
my @hist;
printf("%s %f -> %f\n", $ref, $mins{$ref}, $maxs{$ref});
for ( $i=0; $i<$bins; $i++ ) {
$hist[$i] = 0;
}
seek(FH, 0, 0);
while ( $line = <FH> ) {
my $bin;
my $h;
my $k;
my $l;
my $intensity;
chomp($line);
if ( $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)\s+([\d\.]+)/ ) {
my $key;
if ( $line =~ /\s+(\d+)$/ ) {
my $n = $1;
if ( $n > $max ) {
$line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)\s+/;
$h = $1;
$k = $2;
$l = $3;
$intensity = $4;
$key = $h.",".$k.",".$l;
next unless ( $ref eq $key );
$bin = int(($intensity-$min)/$binsize);
$hist[$bin]++;
$max = $n;
}
}
open(OFH, "> ".$ref.".dat");
for ( $i=0; $i<$bins; $i++ ) {
printf(OFH "%f\t%i\n", $min+(($i+0.5)*$binsize), $hist[$i]);
}
close(OFH);
my $nref = $ref;
$nref =~ s/\,/\ /g;
my $title = "t \"".$nref." reflection\"";
printf(GP "plot [] [0:20]\"".$ref.".dat\" u 1:2 w histeps ".$title."\n");
}
close(FH);
close(GP);
printf("%i %i %i = %i\n", $h, $k, $l, $max);
foreach $ref ( sort {-($counts{$a} <=> $counts{$b})} keys %counts ) {
unlink($ref.".dat");
}
system("ps2pdf histo.ps");
unlink("histo.ps");
exit 0;
src/process_hkl.c
View file @ b7620b55
......@@ -29,8 +29,8 @@
#include "symmetry.h"
/* Number of divisions for R vs |q| graphs */
#define RVQDV (20)
/* Number of divisions for intensity histograms */
#define NBINS (200)
static void show_help(const char *s)
......@@ -56,6 +56,8 @@ 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"
" -g, --histogram=<h,k,l> Calculate the histogram of measurements for this\n"
" reflection.\n"
"\n"
" --scale Scale each pattern for best fit with the current\n"
" model.\n"
......@@ -64,6 +66,47 @@ static void show_help(const char *s)
}
static void plot_histogram(double *vals, int n)
{
int i;
double max = -INFINITY;
double min = +INFINITY;
double step;
int histo[NBINS];
FILE *fh;
fh = fopen("histogram.dat", "w");
if ( fh == NULL ) {
ERROR("Couldn't open 'histogram.dat'\n");
return;
}
for ( i=0; i<n; i++ ) {
if ( vals[i] > max ) max = vals[i];
if ( vals[i] < min ) min = vals[i];
}
STATUS("%f %f\n", min, max);
for ( i=0; i<NBINS; i++ ) {
histo[i] = 0;
}
step = (max-min)/NBINS;
for ( i=0; i<n; i++ ) {
int bin;
bin = (vals[i]-min)/step;
histo[bin]++;
}
for ( i=0; i<NBINS; i++ ) {
fprintf(fh, "%f %i\n", min+step*i, histo[i]);
}
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.
......@@ -160,7 +203,9 @@ static void merge_pattern(double *model, ReflItemList *observed,
const double *new, ReflItemList *items,
unsigned int *model_counts, int mo,
ReflItemList *twins,
const char *holo, const char *mero)
const char *holo, const char *mero, double *hist_vals,
signed int hist_h, signed int hist_k,
signed int hist_l, int *hist_n)
{
int i;
int twin;
......@@ -216,6 +261,14 @@ static void merge_pattern(double *model, ReflItemList *observed,
/* Increase count count */
integrate_count(model_counts, h, k, l, 1);
if ( hist_vals != NULL ) {
int p = *hist_n;
if ( (h==hist_h) && (k==hist_k) && (l==hist_l) ) {
hist_vals[p] = intensity;
*hist_n = p+1;
}
}
}
/* Dump the reflections in this pattern into the overall list */
......@@ -230,7 +283,9 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
int config_maxonly, int config_scale, int config_sum,
int config_stopafter,
ReflItemList *twins, const char *holo, const char *mero,
int n_total_patterns)
int n_total_patterns, double *hist_vals,
signed int hist_h, signed int hist_k, signed int hist_l,
int *hist_i)
{
char *rval;
float f0;
......@@ -277,7 +332,9 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
/* Start of second or later pattern */
merge_pattern(model, observed, new_pattern, items,
counts, config_maxonly,
twins, holo, mero);
twins, holo, mero,
hist_vals, hist_h, hist_k, hist_l,
hist_i);
if ( n_patterns == config_stopafter ) break;
......@@ -380,6 +437,10 @@ int main(int argc, char *argv[])
ReflItemList *observed;
int i;
const char *holo = NULL;
char *histo = NULL;
signed int hist_h, hist_k, hist_l;
double *hist_vals = NULL;
int hist_i;
/* Long options */
const struct option longopts[] = {
......@@ -393,11 +454,12 @@ int main(int argc, char *argv[])
{"scale", 0, &config_scale, 1},
{"symmetry", 1, NULL, 'y'},
{"pdb", 1, NULL, 'p'},
{"histogram", 1, NULL, 'g'},
{0, 0, NULL, 0}
};
/* Short options */
while ((c = getopt_long(argc, argv, "hi:e:ro:p:y:",
while ((c = getopt_long(argc, argv, "hi:e:ro:p:y:g:",
longopts, NULL)) != -1) {
switch (c) {
......@@ -425,6 +487,10 @@ int main(int argc, char *argv[])
sym = strdup(optarg);
break;
case 'g' :
histo = strdup(optarg);
break;
case 0 :
break;
......@@ -475,6 +541,22 @@ int main(int argc, char *argv[])
holo = strdup("1");
}
if ( histo != NULL ) {
int r;
r = sscanf(histo, "%i,%i,%i", &hist_h, &hist_k, &hist_l);
if ( r != 3 ) {
ERROR("Invalid indices for '--histogram'\n");
return 1;
}
hist_vals = malloc(10*1024*sizeof(double));
free(histo);
STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l);
/* Put into the asymmetric cell for the target group */
get_asymm(hist_h, hist_k, hist_l,
&hist_h, &hist_k, &hist_l, sym);
STATUS("%i %i %i\n", hist_h, hist_k, hist_l);
}
/* Open the data stream */
if ( strcmp(filename, "-") == 0 ) {
fh = stdin;
......@@ -492,13 +574,21 @@ int main(int argc, char *argv[])
STATUS("There are %i patterns to process\n", n_total_patterns);
rewind(fh);
hist_i = 0;
merge_all(fh, &model, &observed, &counts,
config_maxonly, config_scale, config_sum, config_stopafter,
twins, holo, sym, n_total_patterns);
twins, holo, sym, n_total_patterns,
hist_vals, hist_h, hist_k, hist_l, &hist_i);
rewind(fh);
fclose(fh);
if ( hist_vals != NULL ) {
STATUS("%i %i %i was seen %i times.\n", hist_h, hist_k, hist_l,
hist_i);
plot_histogram(hist_vals, hist_i);
}
if ( output != NULL ) {
write_reflections(output, observed, model, NULL, counts, cell);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

imprint | privacy | Conditions of Use