Commit 27d14f0b authored by Thomas White's avatar Thomas White
Browse files

Use symmetry when simulating (on the CPU only)

parent 796feb58
......@@ -11,7 +11,7 @@ AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\"
pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c cell.c \
hdf5-file.c detector.c sfac.c peaks.c reflections.c \
beam-parameters.c
beam-parameters.c symmetry.c
if HAVE_OPENCL
pattern_sim_SOURCES += diffraction-gpu.c cl-utils.c
endif
......
......@@ -119,12 +119,13 @@ indexamajig_OBJECTS = $(am_indexamajig_OBJECTS)
indexamajig_DEPENDENCIES =
am__pattern_sim_SOURCES_DIST = pattern_sim.c diffraction.c utils.c \
image.c cell.c hdf5-file.c detector.c sfac.c peaks.c \
reflections.c beam-parameters.c diffraction-gpu.c cl-utils.c
reflections.c beam-parameters.c symmetry.c diffraction-gpu.c \
cl-utils.c
am_pattern_sim_OBJECTS = pattern_sim.$(OBJEXT) diffraction.$(OBJEXT) \
utils.$(OBJEXT) image.$(OBJEXT) cell.$(OBJEXT) \
hdf5-file.$(OBJEXT) detector.$(OBJEXT) sfac.$(OBJEXT) \
peaks.$(OBJEXT) reflections.$(OBJEXT) \
beam-parameters.$(OBJEXT) $(am__objects_1)
beam-parameters.$(OBJEXT) symmetry.$(OBJEXT) $(am__objects_1)
pattern_sim_OBJECTS = $(am_pattern_sim_OBJECTS)
pattern_sim_DEPENDENCIES =
am_powder_plot_OBJECTS = powder_plot.$(OBJEXT) cell.$(OBJEXT) \
......@@ -284,7 +285,7 @@ AM_CFLAGS = -Wall
AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\"
pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c \
cell.c hdf5-file.c detector.c sfac.c peaks.c reflections.c \
beam-parameters.c $(am__append_2)
beam-parameters.c symmetry.c $(am__append_2)
pattern_sim_LDADD = @LIBS@
process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \
reflections.c symmetry.c stream.c beam-parameters.c
......
......@@ -330,7 +330,8 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
/* Setup the OpenCL stuff, create buffers, load the structure factor table */
struct gpu_context *setup_gpu(int no_sfac, struct image *image,
const double *intensities, int dev_num)
const double *intensities, unsigned char *flags,
int dev_num)
{
struct gpu_context *gctx;
cl_uint nplat;
......
......@@ -26,7 +26,8 @@ struct gpu_context;
extern void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
int na, int nb, int nc, UnitCell *ucell);
extern struct gpu_context *setup_gpu(int no_sfac, struct image *image,
const double *intensities, int dev_num);
const double *intensities,
const unsigned char *flags, int dev_num);
extern void cleanup_gpu(struct gpu_context *gctx);
#else
......@@ -39,7 +40,8 @@ static void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
}
static struct gpu_context *setup_gpu(int no_sfac, struct image *image,
const double *intensities, int dev_num)
const double *intensities,
const unsigned char *flags, int dev_num)
{
return NULL;
}
......
......@@ -23,6 +23,7 @@
#include "diffraction.h"
#include "sfac.h"
#include "beam-parameters.h"
#include "symmetry.h"
#define SAMPLING (4)
......@@ -69,8 +70,63 @@ static double lattice_factor(struct rvec q, double ax, double ay, double az,
}
static double interpolate_linear(const double *ref,
float hd, signed int k, signed int l)
static double sym_lookup_intensity(const double *intensities,
const unsigned char *flags, const char *sym,
signed int h, signed int k, signed int l)
{
int i;
double ret = 0.0;
for ( i=0; i<num_general_equivs(sym); i++ ) {
signed int he;
signed int ke;
signed int le;
double f, val;
get_general_equiv(h, k, l, &he, &ke, &le, sym, i);
f = (double)lookup_flag(flags, he, ke, le);
val = lookup_intensity(intensities, he, ke, le);
ret += f*val;
}
return ret;
}
static double sym_lookup_phase(const double *phases,
const unsigned char *flags, const char *sym,
signed int h, signed int k, signed int l)
{
int i;
double ret = 0.0;
for ( i=0; i<num_general_equivs(sym); i++ ) {
signed int he;
signed int ke;
signed int le;
double f, val;
get_general_equiv(h, k, l, &he, &ke, &le, sym, i);
f = (double)lookup_flag(flags, he, ke, le);
val = lookup_phase(phases, he, ke, le);
ret += f*val;
}
return ret;
}
static double interpolate_linear(const double *ref, const unsigned char *flags,
const char *sym, float hd,
signed int k, signed int l)
{
signed int h;
double val1, val2;
......@@ -81,8 +137,8 @@ static double interpolate_linear(const double *ref,
f = hd - (float)h;
assert(f >= 0.0);
val1 = lookup_intensity(ref, h, k, l);
val2 = lookup_intensity(ref, h+1, k, l);
val1 = sym_lookup_intensity(ref, flags, sym, h, k, l);
val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l);
val1 = val1;
val2 = val2;
......@@ -92,6 +148,7 @@ static double interpolate_linear(const double *ref,
static double interpolate_bilinear(const double *ref,
const unsigned char *flags, const char *sym,
float hd, float kd, signed int l)
{
signed int k;
......@@ -103,14 +160,15 @@ static double interpolate_bilinear(const double *ref,
f = kd - (float)k;
assert(f >= 0.0);
val1 = interpolate_linear(ref, hd, k, l);
val2 = interpolate_linear(ref, hd, k+1, l);
val1 = interpolate_linear(ref, flags, sym, hd, k, l);
val2 = interpolate_linear(ref, flags, sym, hd, k+1, l);
return (1.0-f)*val1 + f*val2;
}
static double interpolate_intensity(const double *ref,
const unsigned char *flags, const char *sym,
float hd, float kd, float ld)
{
signed int l;
......@@ -122,8 +180,8 @@ static double interpolate_intensity(const double *ref,
f = ld - (float)l;
assert(f >= 0.0);
val1 = interpolate_bilinear(ref, hd, kd, l);
val2 = interpolate_bilinear(ref, hd, kd, l+1);
val1 = interpolate_bilinear(ref, flags, sym, hd, kd, l);
val2 = interpolate_bilinear(ref, flags, sym, hd, kd, l+1);
return (1.0-f)*val1 + f*val2;
}
......@@ -131,6 +189,8 @@ static double interpolate_intensity(const double *ref,
static double complex interpolate_phased_linear(const double *ref,
const double *phases,
const unsigned char *flags,
const char *sym,
float hd,
signed int k, signed int l)
{
......@@ -146,10 +206,10 @@ static double complex interpolate_phased_linear(const double *ref,
f = hd - (float)h;
assert(f >= 0.0);
val1 = lookup_intensity(ref, h, k, l);
val2 = lookup_intensity(ref, h+1, k, l);
ph1 = lookup_phase(phases, h, k, l);
ph2 = lookup_phase(phases, h+1, k, l);
val1 = sym_lookup_intensity(ref, flags, sym, h, k, l);
val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l);
ph1 = sym_lookup_phase(phases, flags, sym, h, k, l);
ph2 = sym_lookup_phase(phases, flags, sym, h+1, k, l);
val1 = val1;
val2 = val2;
......@@ -169,6 +229,8 @@ static double complex interpolate_phased_linear(const double *ref,
static double complex interpolate_phased_bilinear(const double *ref,
const double *phases,
const unsigned char *flags,
const char *sym,
float hd, float kd,
signed int l)
{
......@@ -181,8 +243,8 @@ static double complex interpolate_phased_bilinear(const double *ref,
f = kd - (float)k;
assert(f >= 0.0);
val1 = interpolate_phased_linear(ref, phases, hd, k, l);
val2 = interpolate_phased_linear(ref, phases, hd, k+1, l);
val1 = interpolate_phased_linear(ref, phases, flags, sym, hd, k, l);
val2 = interpolate_phased_linear(ref, phases, flags, sym, hd, k+1, l);
return (1.0-f)*val1 + f*val2;
}
......@@ -190,6 +252,8 @@ static double complex interpolate_phased_bilinear(const double *ref,
static double interpolate_phased_intensity(const double *ref,
const double *phases,
const unsigned char *flags,
const char *sym,
float hd, float kd, float ld)
{
signed int l;
......@@ -201,20 +265,22 @@ static double interpolate_phased_intensity(const double *ref,
f = ld - (float)l;
assert(f >= 0.0);
val1 = interpolate_phased_bilinear(ref, phases, hd, kd, l);
val2 = interpolate_phased_bilinear(ref, phases, hd, kd, l+1);
val1 = interpolate_phased_bilinear(ref, phases, flags, sym,
hd, kd, l);
val2 = interpolate_phased_bilinear(ref, phases, flags, sym,
hd, kd, l+1);
return cabs((1.0-f)*val1 + f*val2);
}
/* Look up the structure factor for the nearest Bragg condition */
static double molecule_factor(const double *intensities,const double *phases,
struct rvec q,
static double molecule_factor(const double *intensities, const double *phases,
const unsigned char *flags, struct rvec q,
double ax, double ay, double az,
double bx, double by, double bz,
double cx, double cy, double cz,
GradientMethod m)
GradientMethod m, const char *sym)
{
float hd, kd, ld;
signed int h, k, l;
......@@ -224,6 +290,9 @@ static double molecule_factor(const double *intensities,const double *phases,
kd = q.u * bx + q.v * by + q.w * bz;
ld = q.u * cx + q.v * cy + q.w * cz;
/* No flags -> flat intensity distribution */
if ( flags == NULL ) return 1.0e5;
switch ( m ) {
case GRADIENT_MOSAIC :
h = (signed int)rint(hd);
......@@ -232,14 +301,14 @@ static double molecule_factor(const double *intensities,const double *phases,
if ( abs(h) > INDMAX ) r = 0.0;
else if ( abs(k) > INDMAX ) r = 0.0;
else if ( abs(l) > INDMAX ) r = 0.0;
else r = lookup_intensity(intensities, h, k, l);
else r = sym_lookup_intensity(intensities, flags, sym, h, k, l);
break;
case GRADIENT_INTERPOLATE :
r = interpolate_intensity(intensities, hd, kd, ld);
r = interpolate_intensity(intensities, flags, sym, hd, kd, ld);
break;
case GRADIENT_PHASED :
r = interpolate_phased_intensity(intensities, phases,
hd, kd, ld);
r = interpolate_phased_intensity(intensities, phases, flags,
sym, hd, kd, ld);
break;
default:
ERROR("This gradient method not implemented yet.\n");
......@@ -298,7 +367,8 @@ double water_diffraction(struct rvec q, double en,
void get_diffraction(struct image *image, int na, int nb, int nc,
const double *intensities, const double *phases,
UnitCell *cell, int do_water, GradientMethod m)
const unsigned char *flags, UnitCell *cell, int do_water,
GradientMethod m, const char *sym)
{
unsigned int xs, ys;
double ax, ay, az;
......@@ -353,10 +423,10 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
I_molecule = 1.0e10;
} else {
I_molecule = molecule_factor(intensities,
phases, q,
phases, flags, q,
ax,ay,az,
bx,by,bz,cx,cy,cz,
m);
m, sym);
}
I_lattice = pow(f_lattice, 2.0);
......
......@@ -26,7 +26,8 @@ typedef enum {
} GradientMethod;
extern void get_diffraction(struct image *image, int na, int nb, int nc,
const double *intensities,const double *phases,
UnitCell *cell, int do_water, GradientMethod m);
const double *intensities, const double *phases,
const unsigned char *flags, UnitCell *cell,
int do_water, GradientMethod m, const char *sym);
#endif /* DIFFRACTION_H */
......@@ -70,6 +70,8 @@ struct static_index_args
IndexingMethod indm;
IndexingPrivate *ipriv;
const double *intensities;
const unsigned char *flags;
const char *sym; /* Symmetry of "intensities" and "flags" */
struct gpu_context *gctx;
int gpu_dev;
int peaks;
......@@ -191,6 +193,7 @@ static void show_help(const char *s)
"\nIf you used --simulate, you may also want:\n\n"
" --intensities=<file> Specify file containing reflection intensities\n"
" to use when simulating.\n"
" -y, --symmetry=<sym> The symmetry of the intensities file.\n"
"\n"
"\nOptions for greater performance or verbosity:\n\n"
" --verbose Be verbose about indexing.\n"
......@@ -275,22 +278,23 @@ static struct image *get_simage(struct image *template, int alternate)
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
const double *intensities, UnitCell *cell,
int gpu_dev)
const double *intensities,
const unsigned char *flags, UnitCell *cell,
int gpu_dev, const char *sym)
{
/* Set up GPU if necessary.
* Unfortunately, setup has to go here since until now we don't know
* enough about the situation. */
if ( (gctx != NULL) && (*gctx == NULL) ) {
*gctx = setup_gpu(0, simage, intensities, gpu_dev);
*gctx = setup_gpu(0, simage, intensities, flags, gpu_dev);
}
if ( (gctx != NULL) && (*gctx != NULL) ) {
get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
} else {
get_diffraction(simage, 24, 24, 40,
intensities, NULL, cell, 0,
GRADIENT_MOSAIC);
intensities, NULL, flags, cell, 0,
GRADIENT_MOSAIC, sym);
}
record_image(simage, 0);
......@@ -321,7 +325,9 @@ static void process_image(void *pp, int cookie)
int config_polar = pargs->static_args.config_polar;
IndexingMethod indm = pargs->static_args.indm;
const double *intensities = pargs->static_args.intensities;
const unsigned char *flags = pargs->static_args.flags;
struct gpu_context *gctx = pargs->static_args.gctx;
const char *sym = pargs->static_args.sym;
image.features = NULL;
image.data = NULL;
......@@ -428,13 +434,13 @@ static void process_image(void *pp, int cookie)
if ( config_simulate ) {
if ( config_gpu ) {
pthread_mutex_lock(pargs->static_args.gpu_mutex);
simulate_and_write(simage, &gctx, intensities,
simulate_and_write(simage, &gctx, intensities, flags,
image.indexed_cell,
pargs->static_args.gpu_dev);
pargs->static_args.gpu_dev, sym);
pthread_mutex_unlock(pargs->static_args.gpu_mutex);
} else {
simulate_and_write(simage, NULL, intensities,
image.indexed_cell, 0);
simulate_and_write(simage, NULL, intensities, flags,
image.indexed_cell, 0, sym);
}
}
......@@ -539,6 +545,7 @@ int main(int argc, char *argv[])
char *indm_str = NULL;
UnitCell *cell;
double *intensities = NULL;
unsigned char *flags;
char *intfile = NULL;
char *pdb = NULL;
char *prefix = NULL;
......@@ -557,6 +564,7 @@ int main(int argc, char *argv[])
struct beam_params *beam = NULL;
double nominal_photon_energy;
int gpu_dev = -1;
char *sym = NULL;
/* Long options */
const struct option longopts[] = {
......@@ -579,6 +587,7 @@ int main(int argc, char *argv[])
{"verbose", 0, &config_verbose, 1},
{"alternate", 0, &config_alternate, 1},
{"intensities", 1, NULL, 'q'},
{"symmetry", 1, NULL, 'y'},
{"pdb", 1, NULL, 'p'},
{"prefix", 1, NULL, 'x'},
{"unpolarized", 0, &config_polar, 0},
......@@ -595,7 +604,7 @@ int main(int argc, char *argv[])
};
/* Short options */
while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:t:o:b:",
while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:t:o:b:y:",
longopts, NULL)) != -1) {
switch (c) {
......@@ -639,6 +648,10 @@ int main(int argc, char *argv[])
threshold = strtof(optarg, NULL);
break;
case 'y' :
sym = strdup(optarg);
break;
case 'b' :
beam = get_beam_parameters(optarg);
if ( beam == NULL ) {
......@@ -701,6 +714,8 @@ int main(int argc, char *argv[])
}
free(outfile);
if ( sym == NULL ) sym = strdup("1");
if ( speaks == NULL ) {
speaks = strdup("zaef");
STATUS("You didn't specify a peak detection method.\n");
......@@ -717,12 +732,26 @@ int main(int argc, char *argv[])
free(speaks);
if ( intfile != NULL ) {
ReflItemList *items;
int i;
items = read_reflections(intfile, intensities,
NULL, NULL, NULL);
flags = new_list_flag();
for ( i=0; i<num_items(items); i++ ) {
struct refl_item *it = get_item(items, i);
set_flag(flags, it->h, it->k, it->l, 1);
}
delete_items(items);
} else {
intensities = NULL;
flags = NULL;
}
if ( pdb == NULL ) {
......@@ -859,6 +888,8 @@ int main(int argc, char *argv[])
qargs.static_args.indm = indm;
qargs.static_args.ipriv = ipriv;
qargs.static_args.intensities = intensities;
qargs.static_args.flags = flags;
qargs.static_args.sym = sym;
qargs.static_args.gctx = gctx;
qargs.static_args.gpu_dev = gpu_dev;
qargs.static_args.peaks = peaks;
......@@ -881,6 +912,7 @@ int main(int argc, char *argv[])
free(det);
cell_free(cell);
if ( fh != stdout ) fclose(fh);
free(sym);
STATUS("There were %i images. %i could be indexed, of which %i"
" looked sane.\n", n_images, qargs.n_indexable, qargs.n_sane);
......
......@@ -68,6 +68,7 @@ static void show_help(const char *s)
" this invocation as the given filename.\n"
" -i, --intensities=<file> Specify file containing reflection intensities\n"
" (and phases) to use.\n"
" -y, --symmetry=<sym> The symmetry of the intensities file.\n"
" -t, --gradients=<method> Use <method> for the calculation of shape\n"
" transform intensities. Choose from:\n"
" mosaic : Take the intensity of the nearest\n"
......@@ -201,6 +202,7 @@ int main(int argc, char *argv[])
double *intensities;
char *rval;
double *phases;
unsigned char *flags;
int config_simdetails = 0;
int config_nearbragg = 0;
int config_randomquat = 0;
......@@ -227,6 +229,7 @@ int main(int argc, char *argv[])
int random_size = 0;
double min_size = 0.0;
double max_size = 0.0;
char *sym = NULL;
/* Long options */
const struct option longopts[] = {
......@@ -240,6 +243,7 @@ int main(int argc, char *argv[])
{"no-water", 0, &config_nowater, 1},
{"no-noise", 0, &config_nonoise, 1},
{"intensities", 1, NULL, 'i'},
{"symmetry", 1, NULL, 'y'},
{"powder", 1, NULL, 'w'},
{"gradients", 1, NULL, 't'},
{"pdb", 1, NULL, 'p'},
......@@ -254,7 +258,7 @@ int main(int argc, char *argv[])
};
/* Short options */
while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:g:b:",
while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:g:b:y:",
longopts, NULL)) != -1) {
switch (c) {
......@@ -302,6 +306,10 @@ int main(int argc, char *argv[])
beamfile = strdup(optarg);
break;
case 'y' :
sym = strdup(optarg);
break;
case 2 :
gpu_dev = atoi(optarg);
break;
......@@ -359,6 +367,8 @@ int main(int argc, char *argv[])
}
}
if ( sym == NULL ) sym = strdup("1");
if ( config_simdetails ) {
show_details();
return 0;
......@@ -404,14 +414,20 @@ int main(int argc, char *argv[])
}
if ( intfile == NULL ) {
/* Gentle reminder */
STATUS("You didn't specify the file containing the ");
STATUS("reflection intensities (with --intensities).\n");
STATUS("I'll simulate a flat intensity distribution.\n");
intensities = NULL;
phases = NULL;
flags = NULL;
} else {
int i;
ReflItemList *items;
if ( grad == GRADIENT_PHASED ) {
phases = new_list_phase();
} else {
......@@ -419,9 +435,16 @@ int main(int argc, char *argv[])
}
intensities = new_list_intensity();
phases = new_list_phase();
flags = new_list_flag();
items = read_reflections(intfile, intensities, phases,
NULL, NULL);
free(intfile);
for ( i=0; i<num_items(items); i++ ) {
struct refl_item *it = get_item(items, i);
set_flag(flags, it->h, it->k, it->l, 1);
}
delete_items(items);
}
......@@ -516,12 +539,13 @@ int main(int argc, char *argv[])
if ( config_gpu ) {
if ( gctx == NULL ) {
gctx = setup_gpu(config_nosfac, &image,
intensities, gpu_dev);
intensities, flags, gpu_dev);
}
get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
} else {
get_diffraction(&image, na, nb, nc, intensities, phases,
cell, !config_nowater, grad);
flags, cell, !config_nowater, grad,
sym);
}
if ( image.data == NULL ) {
ERROR("Diffraction calculation failed.\n");
......@@ -599,6 +623,7 @@ skip:
free(intensities);
free(outfile);
free(filename);
free(sym);
return 0;
}
......@@ -180,6 +180,11 @@ static inline double angle_between(double x1, double y1, double z1,
#define TYPE double
#include "list_tmp.h"
/* As above, but for simple flags */
#define LABEL(x) x##_flag
#define TYPE unsigned char
#include "list_tmp.h"
/* ----------- Reflection lists indexed by sequence (not indices) ----------- */
......