diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/cl-utils.c | 354 | ||||
| -rw-r--r-- | src/cl-utils.h | 55 | ||||
| -rw-r--r-- | src/diffraction-gpu.c | 596 | ||||
| -rw-r--r-- | src/diffraction-gpu.h | 77 | ||||
| -rw-r--r-- | src/diffraction.c | 502 | ||||
| -rw-r--r-- | src/diffraction.cl.h | 513 | ||||
| -rw-r--r-- | src/diffraction.h | 66 | ||||
| -rw-r--r-- | src/list_tmp.h | 104 | ||||
| -rw-r--r-- | src/partial_sim.c | 1034 | ||||
| -rw-r--r-- | src/pattern_sim.c | 1183 | ||||
| -rw-r--r-- | src/pattern_sim.h | 67 |
11 files changed, 0 insertions, 4551 deletions
diff --git a/src/cl-utils.c b/src/cl-utils.c deleted file mode 100644 index 1b6667c4..00000000 --- a/src/cl-utils.c +++ /dev/null @@ -1,354 +0,0 @@ -/* - * cl-utils.c - * - * OpenCL utility functions - * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2010-2020 Thomas White <taw@physics.org> - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <stdlib.h> -#include <stdio.h> -#include <string.h> - -#define CL_TARGET_OPENCL_VERSION 220 -#ifdef HAVE_CL_CL_H -#include <CL/cl.h> -#else -#include <cl.h> -#endif - -#include "utils.h" - - -/* Return 1 if a GPU device is present, 0 if not, 2 on error. */ -int have_gpu_device() -{ - cl_uint nplat; - cl_platform_id platforms[8]; - cl_context_properties prop[3]; - cl_int err; - int i; - - err = clGetPlatformIDs(8, platforms, &nplat); - if ( err != CL_SUCCESS ) return 2; - if ( nplat == 0 ) return 0; - - /* Find a GPU platform in the list */ - for ( i=0; i<nplat; i++ ) { - - prop[0] = CL_CONTEXT_PLATFORM; - prop[1] = (cl_context_properties)platforms[i]; - prop[2] = 0; - - clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU, - NULL, NULL, &err); - - if ( err != CL_SUCCESS ) { - if ( err != CL_DEVICE_NOT_FOUND ) return 2; - } else { - return 1; - } - } - - return 0; -} - - -const char *clError(cl_int err) -{ - switch ( err ) { - - case CL_SUCCESS : - return "no error"; - - case CL_DEVICE_NOT_AVAILABLE : - return "device not available"; - - case CL_DEVICE_NOT_FOUND : - return "device not found"; - - case CL_INVALID_DEVICE_TYPE : - return "invalid device type"; - - case CL_INVALID_PLATFORM : - return "invalid platform"; - - case CL_INVALID_KERNEL : - return "invalid kernel"; - - case CL_INVALID_ARG_INDEX : - return "invalid argument index"; - - case CL_INVALID_ARG_VALUE : - return "invalid argument value"; - - case CL_INVALID_MEM_OBJECT : - return "invalid memory object"; - - case CL_INVALID_SAMPLER : - return "invalid sampler"; - - case CL_INVALID_ARG_SIZE : - return "invalid argument size"; - - case CL_INVALID_COMMAND_QUEUE : - return "invalid command queue"; - - case CL_INVALID_CONTEXT : - return "invalid context"; - - case CL_INVALID_VALUE : - return "invalid value"; - - case CL_INVALID_EVENT_WAIT_LIST : - return "invalid wait list"; - - case CL_MAP_FAILURE : - return "map failure"; - - case CL_MEM_OBJECT_ALLOCATION_FAILURE : - return "object allocation failure"; - - case CL_OUT_OF_HOST_MEMORY : - return "out of host memory"; - - case CL_OUT_OF_RESOURCES : - return "out of resources"; - - case CL_INVALID_KERNEL_NAME : - return "invalid kernel name"; - - case CL_INVALID_KERNEL_ARGS : - return "invalid kernel arguments"; - - case CL_INVALID_WORK_GROUP_SIZE : - return "invalid work group size"; - - case CL_IMAGE_FORMAT_NOT_SUPPORTED : - return "image format not supported"; - - case CL_INVALID_WORK_DIMENSION : - return "invalid work dimension"; - - default : - return "unknown error"; - } -} - - -static char *get_device_string(cl_device_id dev, cl_device_info info) -{ - int r; - size_t size; - char *val; - - r = clGetDeviceInfo(dev, info, 0, NULL, &size); - if ( r != CL_SUCCESS ) { - ERROR("Couldn't get device vendor size: %s\n", - clError(r)); - return NULL; - } - val = malloc(size); - r = clGetDeviceInfo(dev, info, size, val, NULL); - if ( r != CL_SUCCESS ) { - ERROR("Couldn't get dev vendor: %s\n", clError(r)); - return NULL; - } - - return val; -} - - -cl_device_id get_cl_dev(cl_context ctx, int n) -{ - cl_device_id *dev; - cl_int r; - size_t size; - int num_devs; - - /* Get the required size of the array */ - r = clGetContextInfo(ctx, CL_CONTEXT_DEVICES, 0, NULL, &size); - if ( r != CL_SUCCESS ) { - ERROR("Couldn't get array size for devices: %s\n", clError(r)); - return 0; - } - - dev = malloc(size); - r = clGetContextInfo(ctx, CL_CONTEXT_DEVICES, size, dev, NULL); - if ( r != CL_SUCCESS ) { - ERROR("Couldn't get device: %s\n", clError(r)); - return 0; - } - num_devs = size/sizeof(cl_device_id); - - if ( n >= num_devs ) { - ERROR("Device ID out of range\n"); - return 0; - } - - if ( n < 0 ) { - - int i; - - STATUS("Available devices:\n"); - for ( i=0; i<num_devs; i++ ) { - - char *vendor; - char *name; - - vendor = get_device_string(dev[i], CL_DEVICE_VENDOR); - name = get_device_string(dev[i], CL_DEVICE_NAME); - - STATUS("Device %i: %s %s\n", i, vendor, name); - - } - n = 0; - - STATUS("Using device 0. Use --gpu-dev to choose another.\n"); - - } else { - - char *vendor; - char *name; - - vendor = get_device_string(dev[n], CL_DEVICE_VENDOR); - name = get_device_string(dev[n], CL_DEVICE_NAME); - - STATUS("Using device %i: %s %s\n", n, vendor, name); - - } - - return dev[n]; -} - - -static void show_build_log(cl_program prog, cl_device_id dev) -{ - cl_int r; - char log[4096]; - size_t s; - - r = clGetProgramBuildInfo(prog, dev, CL_PROGRAM_BUILD_LOG, 4096, log, - &s); - - STATUS("Status: %i\n", r); - STATUS("%s\n", log); -} - - -cl_program load_program_from_string(const char *source_in, size_t len, - cl_context ctx, cl_device_id dev, - cl_int *err, const char *extra_cflags, - const char *insert_stuff) -{ - cl_program prog; - cl_int r; - char cflags[1024] = ""; - char *insert_pos; - char *source; - - /* Copy the code because we need to zero-terminate it */ - source = malloc(len+1); - if ( source == NULL ) return 0; - memcpy(source, source_in, len); - source[len] = '\0'; - - if ( insert_stuff != NULL ) { - insert_pos = strstr(source, "INSERT_HERE"); - - if ( insert_pos != NULL ) { - - char *source2; - size_t il; - - source2 = malloc(strlen(source)+strlen(insert_stuff)+1); - if ( source2 == NULL ) { - free(source); - return 0; - } - - il = insert_pos - source; - memcpy(source2, source, il); - memcpy(source2+il, insert_stuff, strlen(insert_stuff)+1); - memcpy(source2+il+strlen(insert_stuff), - source+il+11, strlen(source+il+11)+1); - free(source); - source = source2; - - } - } - - prog = clCreateProgramWithSource(ctx, 1, (const char **)&source, - NULL, err); - if ( *err != CL_SUCCESS ) { - ERROR("Couldn't load source\n"); - return 0; - } - - cflags[0] = '\0'; - strncat(cflags, "-cl-no-signed-zeros ", 1023-strlen(cflags)); - strncat(cflags, extra_cflags, 1023-strlen(cflags)); - - r = clBuildProgram(prog, 0, NULL, cflags, NULL, NULL); - if ( r != CL_SUCCESS ) { - ERROR("Couldn't build program\n"); - show_build_log(prog, dev); - *err = r; - return 0; - } - - free(source); - *err = CL_SUCCESS; - return prog; -} - - -cl_program load_program(const char *filename, cl_context ctx, - cl_device_id dev, cl_int *err, const char *extra_cflags, - const char *insert_stuff) -{ - FILE *fh; - char *source; - size_t len; - - fh = fopen(filename, "r"); - if ( fh == NULL ) { - ERROR("Couldn't open '%s'\n", filename); - *err = CL_INVALID_PROGRAM; - return 0; - } - source = malloc(16384); - if ( source == NULL ) { - fclose(fh); - return 0; - } - len = fread(source, 1, 16383, fh); - fclose(fh); - - return load_program_from_string(source, len, ctx, dev,err, extra_cflags, - insert_stuff); -} diff --git a/src/cl-utils.h b/src/cl-utils.h deleted file mode 100644 index d72ee6d5..00000000 --- a/src/cl-utils.h +++ /dev/null @@ -1,55 +0,0 @@ -/* - * cl-utils.h - * - * OpenCL utility functions - * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2010-2020 Thomas White <taw@physics.org> - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - -#ifndef CLUTILS_H -#define CLUTILS_H - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#define CL_TARGET_OPENCL_VERSION 220 -#ifdef HAVE_CL_CL_H -#include <CL/cl.h> -#else -#include <cl.h> -#endif - -extern int have_gpu_device(void); -extern const char *clError(cl_int err); -extern cl_device_id get_cl_dev(cl_context ctx, int n); -extern cl_program load_program_from_string(const char *source_in, size_t len, - cl_context ctx, cl_device_id dev, - cl_int *err, const char *extra_cflags, - const char *insert_stuff); -extern cl_program load_program(const char *filename, cl_context ctx, - cl_device_id dev, cl_int *err, - const char *extra_cflags, const char *insert_stuff); - - -#endif /* CLUTILS_H */ diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c deleted file mode 100644 index 97badc36..00000000 --- a/src/diffraction-gpu.c +++ /dev/null @@ -1,596 +0,0 @@ -/* - * diffraction-gpu.c - * - * Calculate diffraction patterns by Fourier methods (GPU version) - * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2009-2020 Thomas White <taw@physics.org> - * 2013 Alexandra Tolstikova - * 2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de> - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <stdlib.h> -#include <math.h> -#include <stdio.h> -#include <string.h> -#include <complex.h> - -#define CL_TARGET_OPENCL_VERSION 220 -#ifdef HAVE_CL_CL_H -#include <CL/cl.h> -#else -#include <cl.h> -#endif - -#include "image.h" -#include "utils.h" -#include "cell.h" -#include "diffraction.h" -#include "cl-utils.h" -#include "pattern_sim.h" - -#include "diffraction.cl.h" - -#define SINC_LUT_ELEMENTS (4096) - - -struct gpu_context -{ - cl_context ctx; - cl_command_queue cq; - cl_program prog; - cl_kernel kern; - cl_mem intensities; - cl_mem flags; - - /* Array of sinc LUTs */ - cl_mem *sinc_luts; - cl_float **sinc_lut_ptrs; - int max_sinc_lut; /* Number of LUTs, i.e. one greater than the maximum - * index. This equals the highest allowable "n". */ -}; - - -static void check_sinc_lut(struct gpu_context *gctx, int n, - int no_fringes, int flat) -{ - cl_int err; - cl_image_format fmt; - int i; - - if ( n > gctx->max_sinc_lut ) { - - gctx->sinc_luts = realloc(gctx->sinc_luts, - n*sizeof(*gctx->sinc_luts)); - gctx->sinc_lut_ptrs = realloc(gctx->sinc_lut_ptrs, - n*sizeof(*gctx->sinc_lut_ptrs)); - - for ( i=gctx->max_sinc_lut; i<n; i++ ) { - gctx->sinc_lut_ptrs[i] = NULL; - } - - gctx->max_sinc_lut = n; - } - - if ( gctx->sinc_lut_ptrs[n-1] != NULL ) return; - - /* Create a new sinc LUT */ - gctx->sinc_lut_ptrs[n-1] = malloc(SINC_LUT_ELEMENTS*sizeof(cl_float)); - gctx->sinc_lut_ptrs[n-1][0] = n; - if ( n == 1 ) { - for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) { - gctx->sinc_lut_ptrs[n-1][i] = 1.0; - } - } else { - for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) { - double x, val; - x = (double)i/SINC_LUT_ELEMENTS; - if ( (flat || no_fringes) && (x > 1.0/n) && (1.0-x > 1.0/n) ) { - val = 0.0; - } else if ( flat ) { - val = n; - } else { - val = fabs(sin(M_PI*n*x)/sin(M_PI*x)); - } - gctx->sinc_lut_ptrs[n-1][i] = val; - } - } - - fmt.image_channel_order = CL_INTENSITY; - fmt.image_channel_data_type = CL_FLOAT; - - gctx->sinc_luts[n-1] = clCreateImage2D(gctx->ctx, - CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, - &fmt, SINC_LUT_ELEMENTS, 1, 0, - gctx->sinc_lut_ptrs[n-1], &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't create LUT for %i\n", n); - return; - } -} - - -static int set_arg_float(struct gpu_context *gctx, int idx, float val) -{ - cl_int err; - err = clSetKernelArg(gctx->kern, idx, sizeof(cl_float), &val); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't set kernel argument %i: %s\n", - idx, clError(err)); - return 1; - } - - return 0; -} - - -static int set_arg_int(struct gpu_context *gctx, int idx, int val) -{ - cl_int err; - - err = clSetKernelArg(gctx->kern, idx, sizeof(cl_int), &val); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't set kernel argument %i: %s\n", - idx, clError(err)); - return 1; - } - - return 0; -} - - -static int set_arg_mem(struct gpu_context *gctx, int idx, cl_mem val) -{ - cl_int err; - - err = clSetKernelArg(gctx->kern, idx, sizeof(cl_mem), &val); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't set kernel argument %i: %s\n", - idx, clError(err)); - return 1; - } - - return 0; -} - - -static int do_panels(struct gpu_context *gctx, struct image *image, - double k, double weight, - int *n_inf, int *n_neg, int *n_nan) -{ - int i; - const int sampling = 4; /* This, squared, number of samples / pixel */ - - if ( set_arg_float(gctx, 1, k) ) return 1; - if ( set_arg_float(gctx, 2, weight) ) return 1; - - /* Iterate over panels */ - for ( i=0; i<image->detgeom->n_panels; i++ ) { - - size_t dims[2]; - size_t ldims[2]; - struct detgeom_panel *p; - cl_mem diff; - size_t diff_size; - float *diff_ptr; - int fs, ss; - cl_int err; - - p = &image->detgeom->panels[i]; - - /* Buffer for the results of this panel */ - diff_size = p->w * p->h * sizeof(cl_float); - diff = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY, - diff_size, NULL, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't allocate diffraction memory\n"); - return 1; - } - - if ( set_arg_mem(gctx, 0, diff) ) return 1; - - if ( set_arg_int(gctx, 3, p->w) ) return 1; - if ( set_arg_float(gctx, 4, p->cnx) ) return 1; - if ( set_arg_float(gctx, 5, p->cny) ) return 1; - if ( set_arg_float(gctx, 6, p->fsx) ) return 1; - if ( set_arg_float(gctx, 7, p->fsy) ) return 1; - if ( set_arg_float(gctx, 8, p->fsz) ) return 1; - if ( set_arg_float(gctx, 9, p->ssx) ) return 1; - if ( set_arg_float(gctx, 10, p->ssy) ) return 1; - if ( set_arg_float(gctx, 11, p->ssz) ) return 1; - if ( set_arg_float(gctx, 12, 1.0/p->pixel_pitch) ) return 1; - if ( set_arg_float(gctx, 13, p->cnz*p->pixel_pitch) ) return 1; - - dims[0] = p->w * sampling; - dims[1] = p->h * sampling; - - ldims[0] = sampling; - ldims[1] = sampling; - - err = clSetKernelArg(gctx->kern, 20, - sampling*sampling*sizeof(cl_float), NULL); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't set local memory: %s\n", clError(err)); - return 1; - } - - err = clEnqueueNDRangeKernel(gctx->cq, gctx->kern, 2, NULL, - dims, ldims, 0, NULL, NULL); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't enqueue diffraction kernel: %s\n", - clError(err)); - return 1; - } - - clFinish(gctx->cq); - - diff_ptr = clEnqueueMapBuffer(gctx->cq, diff, CL_TRUE, - CL_MAP_READ, 0, diff_size, - 0, NULL, NULL, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't map diffraction buffer: %s\n", - clError(err)); - return 1; - } - - for ( ss=0; ss<p->h; ss++ ) { - for ( fs=0; fs<p->w; fs++ ) { - - float val; - - val = diff_ptr[fs + p->w*ss]; - if ( isinf(val) ) (*n_inf)++; - if ( val < 0.0 ) (*n_neg)++; - if ( isnan(val) ) (*n_nan)++; - - image->dp[i][fs + p->w*ss] += val; - - } - } - - clEnqueueUnmapMemObject(gctx->cq, diff, diff_ptr, - 0, NULL, NULL); - - clReleaseMemObject(diff); - - } - - return 0; -} - - -int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, - int na, int nb, int nc, UnitCell *ucell, - int no_fringes, int flat, int n_samples) -{ - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - cl_float16 cell; - cl_int err; - int n_inf = 0; - int n_neg = 0; - int n_nan = 0; - int i; - double kmin, kmax, step, norm; - - if ( gctx == NULL ) { - ERROR("GPU setup failed.\n"); - return 1; - } - - /* Ensure all required LUTs are available */ - check_sinc_lut(gctx, na, no_fringes, flat); - check_sinc_lut(gctx, nb, no_fringes, flat); - check_sinc_lut(gctx, nc, no_fringes, flat); - - /* Unit cell */ - cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - cell.s[0] = ax; cell.s[1] = ay; cell.s[2] = az; - cell.s[3] = bx; cell.s[4] = by; cell.s[5] = bz; - cell.s[6] = cx; cell.s[7] = cy; cell.s[8] = cz; - - err = clSetKernelArg(gctx->kern, 14, sizeof(cl_float16), &cell); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't set unit cell: %s\n", clError(err)); - return 1; - } - - if ( set_arg_mem(gctx, 15, gctx->intensities) ) return 1; - if ( set_arg_mem(gctx, 16, gctx->flags) ) return 1; - if ( set_arg_mem(gctx, 17, gctx->sinc_luts[na-1]) ) return 1; - if ( set_arg_mem(gctx, 18, gctx->sinc_luts[nb-1]) ) return 1; - if ( set_arg_mem(gctx, 19, gctx->sinc_luts[nc-1]) ) return 1; - - /* Allocate memory for the result */ - image->dp = malloc(image->detgeom->n_panels * sizeof(float *)); - if ( image->dp == NULL ) { - ERROR("Couldn't allocate memory for result.\n"); - return 1; - } - for ( i=0; i<image->detgeom->n_panels; i++ ) { - struct detgeom_panel *p = &image->detgeom->panels[i]; - image->dp[i] = calloc(p->w * p->h, sizeof(float)); - if ( image->dp[i] == NULL ) { - ERROR("Couldn't allocate memory for panel %i\n", i); - return 1; - } - } - - spectrum_get_range(image->spectrum, &kmin, &kmax); - step = (kmax-kmin)/(n_samples+1); - - /* Determine normalisation factor such that weights add up to 1 after - * sampling (bins must have constant width) */ - norm = 0.0; - for ( i=1; i<=n_samples; i++ ) { - double k = kmin + i*step; - norm += spectrum_get_density_at_k(image->spectrum, k); - } - for ( i=1; i<=n_samples; i++ ) { - - double k = kmin + i*step; - double prob; - - /* Probability = p.d.f. times step width */ - prob = spectrum_get_density_at_k(image->spectrum, k)/norm; - - STATUS("Wavelength: %e m, weight = %.5f\n", 1.0/k, prob); - - err = do_panels(gctx, image, k, prob, &n_inf, &n_neg, &n_nan); - if ( err ) return 1; - - } - - if ( n_neg + n_inf + n_nan ) { - ERROR("WARNING: The GPU calculation produced %i negative" - " values, %i infinities and %i NaNs.\n", - n_neg, n_inf, n_nan); - } - - return 0; -} - - -/* Setup the OpenCL stuff, create buffers, load the structure factor table */ -struct gpu_context *setup_gpu(int no_sfac, - const double *intensities, unsigned char *flags, - const char *sym, int dev_num) -{ - struct gpu_context *gctx; - cl_uint nplat; - cl_platform_id platforms[8]; - cl_context_properties prop[3]; - cl_int err; - cl_device_id dev; - size_t intensities_size; - float *intensities_ptr; - size_t flags_size; - float *flags_ptr; - size_t maxwgsize; - int iplat; - int have_ctx = 0; - char cflags[512] = ""; - char *insert_stuff = NULL; - - STATUS("Setting up GPU...\n"); - - err = clGetPlatformIDs(8, platforms, &nplat); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't get platform IDs: %i\n", err); - return NULL; - } - if ( nplat == 0 ) { - ERROR("Couldn't find at least one platform!\n"); - return NULL; - } - - /* Find a GPU platform in the list */ - for ( iplat=0; iplat<nplat; iplat++ ) { - - prop[0] = CL_CONTEXT_PLATFORM; - prop[1] = (cl_context_properties)platforms[iplat]; - prop[2] = 0; - - gctx = malloc(sizeof(*gctx)); - gctx->ctx = clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU, - NULL, NULL, &err); - - if ( err != CL_SUCCESS ) { - if ( err == CL_DEVICE_NOT_FOUND ) { - /* No GPU device in this platform */ - continue; /* Try next platform */ - } else { - ERROR("Couldn't create OpenCL context: %s (%i)\n", - clError(err), err); - free(gctx); - return NULL; - } - } else { - STATUS("Using OpenCL platform %i (%i total)\n", - iplat, nplat); - have_ctx = 1; - break; - } - } - - if ( !have_ctx ) { - ERROR("Couldn't find a GPU device in any platform.\n"); - return NULL; - } - - dev = get_cl_dev(gctx->ctx, dev_num); - - gctx->cq = clCreateCommandQueue(gctx->ctx, dev, 0, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't create OpenCL command queue\n"); - free(gctx); - return NULL; - } - - /* Create a single-precision version of the scattering factors */ - intensities_size = IDIM*IDIM*IDIM*sizeof(cl_float); - intensities_ptr = malloc(intensities_size); - if ( intensities != NULL ) { - int i; - for ( i=0; i<IDIM*IDIM*IDIM; i++ ) { - intensities_ptr[i] = intensities[i]; - } - } else { - int i; - for ( i=0; i<IDIM*IDIM*IDIM; i++ ) { - intensities_ptr[i] = 100.0; /* Does nothing */ - } - strncat(cflags, "-DFLAT_INTENSITIES ", 511-strlen(cflags)); - } - gctx->intensities = clCreateBuffer(gctx->ctx, - CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, - intensities_size, intensities_ptr, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't allocate intensities memory\n"); - free(gctx); - return NULL; - } - free(intensities_ptr); - - if ( sym != NULL ) { - - int i, n; - SymOpList *pg; - size_t islen = 0; - - insert_stuff = malloc(16384); - if ( insert_stuff == NULL ) return NULL; - insert_stuff[0] = '\0'; - - pg = get_pointgroup(sym); - n = num_equivs(pg, NULL); - for ( i=0; i<n; i++ ) { - - IntegerMatrix *op = get_symop(pg, NULL, i); - char line[1024]; - - snprintf(line, 1023, - "val += lookup_flagged_intensity(intensities, " - "flags, %s, %s, %s);\n\t", - get_matrix_name(op, 0), - get_matrix_name(op, 1), - get_matrix_name(op, 2)); - - islen += strlen(line); - if ( islen > 16383 ) { - ERROR("Too many symmetry operators.\n"); - return NULL; - } - strcat(insert_stuff, line); - - } - - free_symoplist(pg); - - printf("Inserting --->%s<---\n", insert_stuff); - - } else { - if ( intensities != NULL ) { - ERROR("You gave me an intensities file but no point" - " group. I'm assuming '1'.\n"); - strncat(cflags, "-DPG1 ", 511-strlen(cflags)); - } - } - - /* Create a flag array */ - flags_size = IDIM*IDIM*IDIM*sizeof(cl_float); - flags_ptr = malloc(flags_size); - if ( flags != NULL ) { - int i; - for ( i=0; i<IDIM*IDIM*IDIM; i++ ) { - flags_ptr[i] = flags[i]; - } - } else { - int i; - for ( i=0; i<IDIM*IDIM*IDIM; i++ ) { - flags_ptr[i] = 1.0; - } - } - gctx->flags = clCreateBuffer(gctx->ctx, - CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, - flags_size, flags_ptr, &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't allocate flag buffer\n"); - free(gctx); - return NULL; - } - free(flags_ptr); - - gctx->prog = load_program_from_string((char *)data_diffraction_cl, - data_diffraction_cl_len, gctx->ctx, - dev, &err, cflags, insert_stuff); - if ( err != CL_SUCCESS ) { - free(gctx); - return NULL; - } - - gctx->kern = clCreateKernel(gctx->prog, "diffraction", &err); - if ( err != CL_SUCCESS ) { - ERROR("Couldn't create kernel\n"); - free(gctx); - return NULL; - } - - gctx->max_sinc_lut = 0; - gctx->sinc_lut_ptrs = NULL; - gctx->sinc_luts = NULL; - - clGetDeviceInfo(dev, CL_DEVICE_MAX_WORK_GROUP_SIZE, - sizeof(size_t), &maxwgsize, NULL); - STATUS("Maximum work group size = %lli\n", (long long int)maxwgsize); - - return gctx; -} - - -void cleanup_gpu(struct gpu_context *gctx) -{ - int i; - - clReleaseProgram(gctx->prog); - clReleaseMemObject(gctx->intensities); - - /* Release LUTs */ - for ( i=1; i<=gctx->max_sinc_lut; i++ ) { - if ( gctx->sinc_lut_ptrs[i-1] != NULL ) { - clReleaseMemObject(gctx->sinc_luts[i-1]); - free(gctx->sinc_lut_ptrs[i-1]); - } - } - - free(gctx->sinc_luts); - free(gctx->sinc_lut_ptrs); - - clReleaseCommandQueue(gctx->cq); - clReleaseContext(gctx->ctx); - free(gctx); -} diff --git a/src/diffraction-gpu.h b/src/diffraction-gpu.h deleted file mode 100644 index 1128c8b1..00000000 --- a/src/diffraction-gpu.h +++ /dev/null @@ -1,77 +0,0 @@ -/* - * diffraction-gpu.h - * - * Calculate diffraction patterns by Fourier methods (GPU version) - * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2010-2019 Thomas White <taw@physics.org> - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#ifndef DIFFRACTION_GPU_H -#define DIFFRACTION_GPU_H - -#include "image.h" -#include "cell.h" - -struct gpu_context; - -#ifdef HAVE_OPENCL - -extern int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, - int na, int nb, int nc, UnitCell *ucell, - int no_fringes, int flat, int n_samples); -extern struct gpu_context *setup_gpu(int no_sfac, - const double *intensities, - const unsigned char *flags, - const char *sym, int dev_num); -extern void cleanup_gpu(struct gpu_context *gctx); - -#else - -static int get_diffraction_gpu(struct gpu_context *gctx, struct image *image, - int na, int nb, int nc, UnitCell *ucell, - int no_fringes, int flat, int n_samples) -{ - /* Do nothing */ - ERROR("This copy of CrystFEL was not compiled with OpenCL support.\n"); - return 1; -} - -static struct gpu_context *setup_gpu(int no_sfac, - const double *intensities, - const unsigned char *flags, - const char *sym, int dev_num) -{ - return NULL; -} - -static void cleanup_gpu(struct gpu_context *gctx) -{ -} - -#endif - -#endif /* DIFFRACTION_GPU_H */ diff --git a/src/diffraction.c b/src/diffraction.c deleted file mode 100644 index c0dba8aa..00000000 --- a/src/diffraction.c +++ /dev/null @@ -1,502 +0,0 @@ -/* - * diffraction.c - * - * Calculate diffraction patterns by Fourier methods - * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2009-2020 Thomas White <taw@physics.org> - * 2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de> - * 2013 Alexandra Tolstikova - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - - -#include <stdlib.h> -#include <math.h> -#include <stdio.h> -#include <string.h> -#include <complex.h> -#include <assert.h> -#include <fenv.h> - -#include "image.h" -#include "utils.h" -#include "cell.h" -#include "diffraction.h" -#include "symmetry.h" -#include "pattern_sim.h" - - -#define SINC_LUT_ELEMENTS (4096) - - -static double *get_sinc_lut(int n, int no_fringes, int flat) -{ - int i; - double *lut; - - lut = malloc(SINC_LUT_ELEMENTS*sizeof(double)); - lut[0] = n; - if ( n == 1 ) { - for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) { - lut[i] = 1.0; - } - } else { - for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) { - double x, val; - x = (double)i/SINC_LUT_ELEMENTS; - if ( (flat || no_fringes) && (x > 1.0/n) && (1.0-x > 1.0/n) ) { - val = 0.0; - } else if ( flat ) { - val = n; - } else { - val = fabs(sin(M_PI*n*x)/sin(M_PI*x)); - } - lut[i] = val; - } - } - - return lut; -} - - -static double interpolate_lut(double *lut, double val) -{ - double i, pos, f; - unsigned int low, high; - - pos = SINC_LUT_ELEMENTS * modf(fabs(val), &i); - low = (int)pos; /* Discard fractional part */ - high = low + 1; - f = modf(pos, &i); /* Fraction */ - if ( high == SINC_LUT_ELEMENTS ) high = 0; - - return (1.0-f)*lut[low] + f*lut[high]; -} - - -static double lattice_factor(struct rvec q, double ax, double ay, double az, - double bx, double by, double bz, - double cx, double cy, double cz, - double *lut_a, double *lut_b, - double *lut_c) -{ - struct rvec Udotq; - double f1, f2, f3; - - Udotq.u = ax*q.u + ay*q.v + az*q.w; - Udotq.v = bx*q.u + by*q.v + bz*q.w; - Udotq.w = cx*q.u + cy*q.v + cz*q.w; - - f1 = interpolate_lut(lut_a, Udotq.u); - f2 = interpolate_lut(lut_b, Udotq.v); - f3 = interpolate_lut(lut_c, Udotq.w); - - return f1 * f2 * f3; -} - - -static double sym_lookup_intensity(const double *intensities, - const unsigned char *flags, - const SymOpList *sym, - signed int h, signed int k, signed int l) -{ - int i; - double ret = 0.0; - - for ( i=0; i<num_equivs(sym, NULL); i++ ) { - - signed int he; - signed int ke; - signed int le; - double f, val; - - get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le); - - f = (double)lookup_arr_flag(flags, he, ke, le); - val = lookup_arr_intensity(intensities, he, ke, le); - - ret += f*val; - - } - - return ret; -} - - -static double sym_lookup_phase(const double *phases, - const unsigned char *flags, const SymOpList *sym, - signed int h, signed int k, signed int l) -{ - int i; - - for ( i=0; i<num_equivs(sym, NULL); i++ ) { - - signed int he; - signed int ke; - signed int le; - int f; - - get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le); - - f = lookup_arr_flag(flags, he, ke, le); - - if ( f ) return lookup_arr_phase(phases, he, ke, le); - - } - - return 0.0; -} - - -static double interpolate_linear(const double *ref, const unsigned char *flags, - const SymOpList *sym, float hd, - signed int k, signed int l) -{ - signed int h; - double val1, val2; - float f; - - h = (signed int)hd; - if ( hd < 0.0 ) h -= 1; - f = hd - (float)h; - assert(f >= 0.0); - - val1 = sym_lookup_intensity(ref, flags, sym, h, k, l); - val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l); - - return (1.0-f)*val1 + f*val2; -} - - -static double interpolate_bilinear(const double *ref, - const unsigned char *flags, - const SymOpList *sym, - float hd, float kd, signed int l) -{ - signed int k; - double val1, val2; - float f; - - k = (signed int)kd; - if ( kd < 0.0 ) k -= 1; - f = kd - (float)k; - assert(f >= 0.0); - - 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 SymOpList *sym, - float hd, float kd, float ld) -{ - signed int l; - double val1, val2; - float f; - - l = (signed int)ld; - if ( ld < 0.0 ) l -= 1; - f = ld - (float)l; - assert(f >= 0.0); - - 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; -} - - -static double complex interpolate_phased_linear(const double *ref, - const double *phases, - const unsigned char *flags, - const SymOpList *sym, - float hd, - signed int k, signed int l) -{ - signed int h; - double val1, val2; - float f; - double ph1, ph2; - double re1, re2, im1, im2; - double re, im; - - h = (signed int)hd; - if ( hd < 0.0 ) h -= 1; - f = hd - (float)h; - assert(f >= 0.0); - - 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); - - /* Calculate real and imaginary parts */ - re1 = val1 * cos(ph1); - im1 = val1 * sin(ph1); - re2 = val2 * cos(ph2); - im2 = val2 * sin(ph2); - - re = (1.0-f)*re1 + f*re2; - im = (1.0-f)*im1 + f*im2; - - return re + im*I; -} - - -static double complex interpolate_phased_bilinear(const double *ref, - const double *phases, - const unsigned char *flags, - const SymOpList *sym, - float hd, float kd, - signed int l) -{ - signed int k; - double complex val1, val2; - float f; - - k = (signed int)kd; - if ( kd < 0.0 ) k -= 1; - f = kd - (float)k; - assert(f >= 0.0); - - 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; -} - - -static double interpolate_phased_intensity(const double *ref, - const double *phases, - const unsigned char *flags, - const SymOpList *sym, - float hd, float kd, float ld) -{ - signed int l; - double complex val1, val2; - float f; - - l = (signed int)ld; - if ( ld < 0.0 ) l -= 1; - f = ld - (float)l; - assert(f >= 0.0); - - 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, - 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, const SymOpList *sym) -{ - float hd, kd, ld; - signed int h, k, l; - double r; - - hd = q.u * ax + q.v * ay + q.w * az; - 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 100.0; - - switch ( m ) { - - case GRADIENT_MOSAIC : - fesetround(1); /* Round to nearest */ - h = (signed int)rint(hd); - k = (signed int)rint(kd); - l = (signed int)rint(ld); - 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 = sym_lookup_intensity(intensities, flags, sym, h, k, l); - break; - - case GRADIENT_INTERPOLATE : - r = interpolate_intensity(intensities, flags, sym, hd, kd, ld); - break; - - case GRADIENT_PHASED : - r = interpolate_phased_intensity(intensities, phases, flags, - sym, hd, kd, ld); - break; - - default: - ERROR("This gradient method not implemented yet.\n"); - exit(1); - } - - return r; -} - - -static void diffraction_panel(struct image *image, const double *intensities, - const double *phases, const unsigned char *flags, - UnitCell *cell, GradientMethod m, - const SymOpList *sym, double k, - double ax, double ay, double az, - double bx, double by, double bz, - double cx, double cy, double cz, - double *lut_a, double *lut_b, double *lut_c, - int pn, double weight) -{ - int fs, ss; - const int nxs = 4; - const int nys = 4; - struct detgeom_panel *p = &image->detgeom->panels[pn]; - - weight /= nxs*nys; - - for ( ss=0; ss<p->h; ss++ ) { - for ( fs=0; fs<p->w; fs++ ) { - - int idx; - double f_lattice, I_lattice; - double I_molecule; - int xs, ys; - float xo, yo; - - for ( xs=0; xs<nxs; xs++ ) { - for ( ys=0; ys<nys; ys++ ) { - - double qv[3]; - struct rvec q; - - xo = (1.0/nxs) * xs; - yo = (1.0/nys) * ys; - - detgeom_transform_coords(p, fs+xo, ss+yo, - 1.0/k, 0.0, 0.0, qv); - - q.u = qv[0]; q.v = qv[1]; q.w = qv[2]; - - f_lattice = lattice_factor(q, ax, ay, az, - bx, by, bz, - cx, cy, cz, - lut_a, lut_b, lut_c); - - I_molecule = molecule_factor(intensities, - phases, flags, q, - ax, ay, az, - bx, by, bz, - cx, cy, cz, - m, sym); - - I_lattice = pow(f_lattice, 2.0); - - idx = fs + p->w*ss; - image->dp[pn][idx] += I_lattice * I_molecule * weight; - - } - } - } - progress_bar(ss, p->h-1, "Calculating diffraction"); - } -} - - -static void diffraction_at_k(struct image *image, const double *intensities, - const double *phases, const unsigned char *flags, - UnitCell *cell, GradientMethod m, - const SymOpList *sym, double k, - double ax, double ay, double az, - double bx, double by, double bz, - double cx, double cy, double cz, - double *lut_a, double *lut_b, double *lut_c, - double weight) -{ - int i; - - for ( i=0; i<image->detgeom->n_panels; i++ ) { - diffraction_panel(image, intensities, phases, flags, cell, m, - sym, k, ax, ay, az, bx, by, bz, cx, cy, cz, - lut_a, lut_b, lut_c, i, weight); - } -} - - -void get_diffraction(struct image *image, int na, int nb, int nc, - const double *intensities, const double *phases, - const unsigned char *flags, UnitCell *cell, - GradientMethod m, const SymOpList *sym, - int no_fringes, int flat, int n_samples) -{ - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - double *lut_a; - double *lut_b; - double *lut_c; - int i; - double kmin, kmax, step; - double norm = 0.0; - - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - - lut_a = get_sinc_lut(na, no_fringes, flat); - lut_b = get_sinc_lut(nb, no_fringes, flat); - lut_c = get_sinc_lut(nc, no_fringes, flat); - - spectrum_get_range(image->spectrum, &kmin, &kmax); - step = (kmax-kmin)/(n_samples+1); - - /* Determine normalisation factor such that weights add up to 1 after - * sampling (bins must have constant width) */ - for ( i=1; i<=n_samples; i++ ) { - double k = kmin + i*step; - norm += spectrum_get_density_at_k(image->spectrum, k); - } - for ( i=1; i<=n_samples; i++ ) { - - double k = kmin + i*step; - double prob; - - /* Probability = p.d.f. times step width */ - prob = spectrum_get_density_at_k(image->spectrum, k)/norm; - - STATUS("Wavelength: %e m, weight = %.5f\n", 1.0/k, prob); - - diffraction_at_k(image, intensities, phases, - flags, cell, m, sym, k, - ax, ay, az, bx, by, bz, cx, cy, cz, - lut_a, lut_b, lut_c, prob); - - } - - free(lut_a); - free(lut_b); - free(lut_c); -} diff --git a/src/diffraction.cl.h b/src/diffraction.cl.h deleted file mode 100644 index 1f788b67..00000000 --- a/src/diffraction.cl.h +++ /dev/null @@ -1,513 +0,0 @@ -/* - * This file was generated from data/diffraction.cl - * using the following command: - * xxd -i data/diffraction.cl src/diffraction.cl.h - * - * If you have 'xxd' installed, you can run the script - * data/gen-resources to re-create this file. - */ - -unsigned char data_diffraction_cl[] = { - 0x2f, 0x2a, 0x0a, 0x20, 0x2a, 0x20, 0x64, 0x69, 0x66, 0x66, 0x72, 0x61, - 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x2e, 0x63, 0x6c, 0x0a, 0x20, 0x2a, 0x0a, - 0x20, 0x2a, 0x20, 0x47, 0x50, 0x55, 0x20, 0x63, 0x61, 0x6c, 0x63, 0x75, - 0x6c, 0x61, 0x74, 0x69, 0x6f, 0x6e, 0x20, 0x6b, 0x65, 0x72, 0x6e, 0x65, - 0x6c, 0x20, 0x66, 0x6f, 0x72, 0x20, 0x74, 0x72, 0x75, 0x6e, 0x63, 0x61, - 0x74, 0x65, 0x64, 0x20, 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x20, - 0x64, 0x69, 0x66, 0x66, 0x72, 0x61, 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x0a, - 0x20, 0x2a, 0x0a, 0x20, 0x2a, 0x20, 0x43, 0x6f, 0x70, 0x79, 0x72, 0x69, - 0x67, 0x68, 0x74, 0x20, 0xc2, 0xa9, 0x20, 0x32, 0x30, 0x31, 0x32, 0x2d, - 0x32, 0x30, 0x31, 0x34, 0x20, 0x44, 0x65, 0x75, 0x74, 0x73, 0x63, 0x68, - 0x65, 0x73, 0x20, 0x45, 0x6c, 0x65, 0x6b, 0x74, 0x72, 0x6f, 0x6e, 0x65, - 0x6e, 0x2d, 0x53, 0x79, 0x6e, 0x63, 0x68, 0x72, 0x6f, 0x74, 0x72, 0x6f, - 0x6e, 0x20, 0x44, 0x45, 0x53, 0x59, 0x2c, 0x0a, 0x20, 0x2a, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x61, 0x20, 0x72, - 0x65, 0x73, 0x65, 0x61, 0x72, 0x63, 0x68, 0x20, 0x63, 0x65, 0x6e, 0x74, - 0x72, 0x65, 0x20, 0x6f, 0x66, 0x20, 0x74, 0x68, 0x65, 0x20, 0x48, 0x65, - 0x6c, 0x6d, 0x68, 0x6f, 0x6c, 0x74, 0x7a, 0x20, 0x41, 0x73, 0x73, 0x6f, - 0x63, 0x69, 0x61, 0x74, 0x69, 0x6f, 0x6e, 0x2e, 0x0a, 0x20, 0x2a, 0x0a, - 0x20, 0x2a, 0x20, 0x41, 0x75, 0x74, 0x68, 0x6f, 0x72, 0x73, 0x3a, 0x0a, - 0x20, 0x2a, 0x20, 0x20, 0x20, 0x32, 0x30, 0x30, 0x39, 0x2d, 0x32, 0x30, - 0x31, 0x34, 0x20, 0x54, 0x68, 0x6f, 0x6d, 0x61, 0x73, 0x20, 0x57, 0x68, - 0x69, 0x74, 0x65, 0x20, 0x3c, 0x74, 0x61, 0x77, 0x40, 0x70, 0x68, 0x79, - 0x73, 0x69, 0x63, 0x73, 0x2e, 0x6f, 0x72, 0x67, 0x3e, 0x0a, 0x20, 0x2a, - 0x20, 0x20, 0x20, 0x32, 0x30, 0x31, 0x33, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x41, 0x6c, 0x65, 0x78, 0x61, 0x6e, 0x64, 0x72, 0x61, 0x20, 0x54, - 0x6f, 0x6c, 0x73, 0x74, 0x69, 0x6b, 0x6f, 0x76, 0x61, 0x0a, 0x20, 0x2a, - 0x0a, 0x20, 0x2a, 0x20, 0x54, 0x68, 0x69, 0x73, 0x20, 0x66, 0x69, 0x6c, - 0x65, 0x20, 0x69, 0x73, 0x20, 0x70, 0x61, 0x72, 0x74, 0x20, 0x6f, 0x66, - 0x20, 0x43, 0x72, 0x79, 0x73, 0x74, 0x46, 0x45, 0x4c, 0x2e, 0x0a, 0x20, - 0x2a, 0x0a, 0x20, 0x2a, 0x20, 0x43, 0x72, 0x79, 0x73, 0x74, 0x46, 0x45, - 0x4c, 0x20, 0x69, 0x73, 0x20, 0x66, 0x72, 0x65, 0x65, 0x20, 0x73, 0x6f, - 0x66, 0x74, 0x77, 0x61, 0x72, 0x65, 0x3a, 0x20, 0x79, 0x6f, 0x75, 0x20, - 0x63, 0x61, 0x6e, 0x20, 0x72, 0x65, 0x64, 0x69, 0x73, 0x74, 0x72, 0x69, - 0x62, 0x75, 0x74, 0x65, 0x20, 0x69, 0x74, 0x20, 0x61, 0x6e, 0x64, 0x2f, - 0x6f, 0x72, 0x20, 0x6d, 0x6f, 0x64, 0x69, 0x66, 0x79, 0x0a, 0x20, 0x2a, - 0x20, 0x69, 0x74, 0x20, 0x75, 0x6e, 0x64, 0x65, 0x72, 0x20, 0x74, 0x68, - 0x65, 0x20, 0x74, 0x65, 0x72, 0x6d, 0x73, 0x20, 0x6f, 0x66, 0x20, 0x74, - 0x68, 0x65, 0x20, 0x47, 0x4e, 0x55, 0x20, 0x47, 0x65, 0x6e, 0x65, 0x72, - 0x61, 0x6c, 0x20, 0x50, 0x75, 0x62, 0x6c, 0x69, 0x63, 0x20, 0x4c, 0x69, - 0x63, 0x65, 0x6e, 0x73, 0x65, 0x20, 0x61, 0x73, 0x20, 0x70, 0x75, 0x62, - 0x6c, 0x69, 0x73, 0x68, 0x65, 0x64, 0x20, 0x62, 0x79, 0x0a, 0x20, 0x2a, - 0x20, 0x74, 0x68, 0x65, 0x20, 0x46, 0x72, 0x65, 0x65, 0x20, 0x53, 0x6f, - 0x66, 0x74, 0x77, 0x61, 0x72, 0x65, 0x20, 0x46, 0x6f, 0x75, 0x6e, 0x64, - 0x61, 0x74, 0x69, 0x6f, 0x6e, 0x2c, 0x20, 0x65, 0x69, 0x74, 0x68, 0x65, - 0x72, 0x20, 0x76, 0x65, 0x72, 0x73, 0x69, 0x6f, 0x6e, 0x20, 0x33, 0x20, - 0x6f, 0x66, 0x20, 0x74, 0x68, 0x65, 0x20, 0x4c, 0x69, 0x63, 0x65, 0x6e, - 0x73, 0x65, 0x2c, 0x20, 0x6f, 0x72, 0x0a, 0x20, 0x2a, 0x20, 0x28, 0x61, - 0x74, 0x20, 0x79, 0x6f, 0x75, 0x72, 0x20, 0x6f, 0x70, 0x74, 0x69, 0x6f, - 0x6e, 0x29, 0x20, 0x61, 0x6e, 0x79, 0x20, 0x6c, 0x61, 0x74, 0x65, 0x72, - 0x20, 0x76, 0x65, 0x72, 0x73, 0x69, 0x6f, 0x6e, 0x2e, 0x0a, 0x20, 0x2a, - 0x0a, 0x20, 0x2a, 0x20, 0x43, 0x72, 0x79, 0x73, 0x74, 0x46, 0x45, 0x4c, - 0x20, 0x69, 0x73, 0x20, 0x64, 0x69, 0x73, 0x74, 0x72, 0x69, 0x62, 0x75, - 0x74, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x20, 0x74, 0x68, 0x65, 0x20, 0x68, - 0x6f, 0x70, 0x65, 0x20, 0x74, 0x68, 0x61, 0x74, 0x20, 0x69, 0x74, 0x20, - 0x77, 0x69, 0x6c, 0x6c, 0x20, 0x62, 0x65, 0x20, 0x75, 0x73, 0x65, 0x66, - 0x75, 0x6c, 0x2c, 0x0a, 0x20, 0x2a, 0x20, 0x62, 0x75, 0x74, 0x20, 0x57, - 0x49, 0x54, 0x48, 0x4f, 0x55, 0x54, 0x20, 0x41, 0x4e, 0x59, 0x20, 0x57, - 0x41, 0x52, 0x52, 0x41, 0x4e, 0x54, 0x59, 0x3b, 0x20, 0x77, 0x69, 0x74, - 0x68, 0x6f, 0x75, 0x74, 0x20, 0x65, 0x76, 0x65, 0x6e, 0x20, 0x74, 0x68, - 0x65, 0x20, 0x69, 0x6d, 0x70, 0x6c, 0x69, 0x65, 0x64, 0x20, 0x77, 0x61, - 0x72, 0x72, 0x61, 0x6e, 0x74, 0x79, 0x20, 0x6f, 0x66, 0x0a, 0x20, 0x2a, - 0x20, 0x4d, 0x45, 0x52, 0x43, 0x48, 0x41, 0x4e, 0x54, 0x41, 0x42, 0x49, - 0x4c, 0x49, 0x54, 0x59, 0x20, 0x6f, 0x72, 0x20, 0x46, 0x49, 0x54, 0x4e, - 0x45, 0x53, 0x53, 0x20, 0x46, 0x4f, 0x52, 0x20, 0x41, 0x20, 0x50, 0x41, - 0x52, 0x54, 0x49, 0x43, 0x55, 0x4c, 0x41, 0x52, 0x20, 0x50, 0x55, 0x52, - 0x50, 0x4f, 0x53, 0x45, 0x2e, 0x20, 0x20, 0x53, 0x65, 0x65, 0x20, 0x74, - 0x68, 0x65, 0x0a, 0x20, 0x2a, 0x20, 0x47, 0x4e, 0x55, 0x20, 0x47, 0x65, - 0x6e, 0x65, 0x72, 0x61, 0x6c, 0x20, 0x50, 0x75, 0x62, 0x6c, 0x69, 0x63, - 0x20, 0x4c, 0x69, 0x63, 0x65, 0x6e, 0x73, 0x65, 0x20, 0x66, 0x6f, 0x72, - 0x20, 0x6d, 0x6f, 0x72, 0x65, 0x20, 0x64, 0x65, 0x74, 0x61, 0x69, 0x6c, - 0x73, 0x2e, 0x0a, 0x20, 0x2a, 0x0a, 0x20, 0x2a, 0x20, 0x59, 0x6f, 0x75, - 0x20, 0x73, 0x68, 0x6f, 0x75, 0x6c, 0x64, 0x20, 0x68, 0x61, 0x76, 0x65, - 0x20, 0x72, 0x65, 0x63, 0x65, 0x69, 0x76, 0x65, 0x64, 0x20, 0x61, 0x20, - 0x63, 0x6f, 0x70, 0x79, 0x20, 0x6f, 0x66, 0x20, 0x74, 0x68, 0x65, 0x20, - 0x47, 0x4e, 0x55, 0x20, 0x47, 0x65, 0x6e, 0x65, 0x72, 0x61, 0x6c, 0x20, - 0x50, 0x75, 0x62, 0x6c, 0x69, 0x63, 0x20, 0x4c, 0x69, 0x63, 0x65, 0x6e, - 0x73, 0x65, 0x0a, 0x20, 0x2a, 0x20, 0x61, 0x6c, 0x6f, 0x6e, 0x67, 0x20, - 0x77, 0x69, 0x74, 0x68, 0x20, 0x43, 0x72, 0x79, 0x73, 0x74, 0x46, 0x45, - 0x4c, 0x2e, 0x20, 0x20, 0x49, 0x66, 0x20, 0x6e, 0x6f, 0x74, 0x2c, 0x20, - 0x73, 0x65, 0x65, 0x20, 0x3c, 0x68, 0x74, 0x74, 0x70, 0x3a, 0x2f, 0x2f, - 0x77, 0x77, 0x77, 0x2e, 0x67, 0x6e, 0x75, 0x2e, 0x6f, 0x72, 0x67, 0x2f, - 0x6c, 0x69, 0x63, 0x65, 0x6e, 0x73, 0x65, 0x73, 0x2f, 0x3e, 0x2e, 0x0a, - 0x20, 0x2a, 0x0a, 0x20, 0x2a, 0x2f, 0x0a, 0x0a, 0x0a, 0x2f, 0x2a, 0x20, - 0x4d, 0x61, 0x78, 0x6d, 0x69, 0x6d, 0x75, 0x6d, 0x20, 0x69, 0x6e, 0x64, - 0x65, 0x78, 0x20, 0x74, 0x6f, 0x20, 0x68, 0x6f, 0x6c, 0x64, 0x20, 0x76, - 0x61, 0x6c, 0x75, 0x65, 0x73, 0x20, 0x75, 0x70, 0x20, 0x74, 0x6f, 0x20, - 0x28, 0x63, 0x61, 0x6e, 0x20, 0x62, 0x65, 0x20, 0x69, 0x6e, 0x63, 0x72, - 0x65, 0x61, 0x73, 0x65, 0x64, 0x20, 0x69, 0x66, 0x20, 0x6e, 0x65, 0x63, - 0x65, 0x73, 0x73, 0x61, 0x72, 0x79, 0x29, 0x0a, 0x20, 0x2a, 0x20, 0x57, - 0x41, 0x52, 0x4e, 0x49, 0x4e, 0x47, 0x3a, 0x20, 0x41, 0x6c, 0x74, 0x65, - 0x72, 0x69, 0x6e, 0x67, 0x20, 0x74, 0x68, 0x69, 0x73, 0x20, 0x76, 0x61, - 0x6c, 0x75, 0x65, 0x20, 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x69, 0x74, 0x75, - 0x74, 0x65, 0x73, 0x20, 0x61, 0x6e, 0x20, 0x41, 0x42, 0x49, 0x20, 0x63, - 0x68, 0x61, 0x6e, 0x67, 0x65, 0x2c, 0x20, 0x61, 0x6e, 0x64, 0x20, 0x6d, - 0x65, 0x61, 0x6e, 0x73, 0x20, 0x79, 0x6f, 0x75, 0x20, 0x6d, 0x75, 0x73, - 0x74, 0x0a, 0x20, 0x2a, 0x20, 0x75, 0x70, 0x64, 0x61, 0x74, 0x65, 0x20, - 0x73, 0x72, 0x63, 0x2f, 0x70, 0x61, 0x74, 0x74, 0x65, 0x72, 0x6e, 0x5f, - 0x73, 0x69, 0x6d, 0x2e, 0x68, 0x20, 0x74, 0x68, 0x65, 0x6e, 0x20, 0x72, - 0x65, 0x63, 0x6f, 0x6d, 0x70, 0x69, 0x6c, 0x65, 0x20, 0x61, 0x6e, 0x64, - 0x20, 0x72, 0x65, 0x69, 0x6e, 0x73, 0x74, 0x61, 0x6c, 0x6c, 0x20, 0x65, - 0x76, 0x65, 0x72, 0x79, 0x74, 0x68, 0x69, 0x6e, 0x67, 0x2e, 0x20, 0x2a, - 0x2f, 0x0a, 0x23, 0x64, 0x65, 0x66, 0x69, 0x6e, 0x65, 0x20, 0x49, 0x4e, - 0x44, 0x4d, 0x41, 0x58, 0x20, 0x31, 0x33, 0x30, 0x0a, 0x23, 0x64, 0x65, - 0x66, 0x69, 0x6e, 0x65, 0x20, 0x49, 0x44, 0x49, 0x4d, 0x20, 0x28, 0x49, - 0x4e, 0x44, 0x4d, 0x41, 0x58, 0x2a, 0x32, 0x20, 0x2b, 0x31, 0x29, 0x0a, - 0x0a, 0x23, 0x69, 0x66, 0x6e, 0x64, 0x65, 0x66, 0x20, 0x4d, 0x5f, 0x50, - 0x49, 0x0a, 0x23, 0x64, 0x65, 0x66, 0x69, 0x6e, 0x65, 0x20, 0x4d, 0x5f, - 0x50, 0x49, 0x20, 0x28, 0x28, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x29, 0x28, - 0x33, 0x2e, 0x31, 0x34, 0x31, 0x35, 0x39, 0x32, 0x36, 0x35, 0x29, 0x29, - 0x0a, 0x23, 0x65, 0x6e, 0x64, 0x69, 0x66, 0x0a, 0x0a, 0x63, 0x6f, 0x6e, - 0x73, 0x74, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x74, - 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x61, 0x20, 0x3d, - 0x20, 0x43, 0x4c, 0x4b, 0x5f, 0x4e, 0x4f, 0x52, 0x4d, 0x41, 0x4c, 0x49, - 0x5a, 0x45, 0x44, 0x5f, 0x43, 0x4f, 0x4f, 0x52, 0x44, 0x53, 0x5f, 0x54, - 0x52, 0x55, 0x45, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, 0x43, - 0x4c, 0x4b, 0x5f, 0x41, 0x44, 0x44, 0x52, 0x45, 0x53, 0x53, 0x5f, 0x52, - 0x45, 0x50, 0x45, 0x41, 0x54, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, - 0x20, 0x43, 0x4c, 0x4b, 0x5f, 0x46, 0x49, 0x4c, 0x54, 0x45, 0x52, 0x5f, - 0x4c, 0x49, 0x4e, 0x45, 0x41, 0x52, 0x3b, 0x0a, 0x63, 0x6f, 0x6e, 0x73, - 0x74, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x74, 0x20, - 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x62, 0x20, 0x3d, 0x20, - 0x43, 0x4c, 0x4b, 0x5f, 0x4e, 0x4f, 0x52, 0x4d, 0x41, 0x4c, 0x49, 0x5a, - 0x45, 0x44, 0x5f, 0x43, 0x4f, 0x4f, 0x52, 0x44, 0x53, 0x5f, 0x54, 0x52, - 0x55, 0x45, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, 0x43, 0x4c, - 0x4b, 0x5f, 0x41, 0x44, 0x44, 0x52, 0x45, 0x53, 0x53, 0x5f, 0x52, 0x45, - 0x50, 0x45, 0x41, 0x54, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, - 0x43, 0x4c, 0x4b, 0x5f, 0x46, 0x49, 0x4c, 0x54, 0x45, 0x52, 0x5f, 0x4c, - 0x49, 0x4e, 0x45, 0x41, 0x52, 0x3b, 0x0a, 0x63, 0x6f, 0x6e, 0x73, 0x74, - 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x74, 0x20, 0x73, - 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x63, 0x20, 0x3d, 0x20, 0x43, - 0x4c, 0x4b, 0x5f, 0x4e, 0x4f, 0x52, 0x4d, 0x41, 0x4c, 0x49, 0x5a, 0x45, - 0x44, 0x5f, 0x43, 0x4f, 0x4f, 0x52, 0x44, 0x53, 0x5f, 0x54, 0x52, 0x55, - 0x45, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, 0x43, 0x4c, 0x4b, - 0x5f, 0x41, 0x44, 0x44, 0x52, 0x45, 0x53, 0x53, 0x5f, 0x52, 0x45, 0x50, - 0x45, 0x41, 0x54, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x7c, 0x20, 0x43, - 0x4c, 0x4b, 0x5f, 0x46, 0x49, 0x4c, 0x54, 0x45, 0x52, 0x5f, 0x4c, 0x49, - 0x4e, 0x45, 0x41, 0x52, 0x3b, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f, 0x61, - 0x74, 0x34, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x71, 0x28, 0x66, 0x6c, 0x6f, - 0x61, 0x74, 0x20, 0x66, 0x73, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, - 0x20, 0x73, 0x73, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x72, - 0x65, 0x73, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6c, - 0x65, 0x6e, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x6b, 0x2c, - 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6f, 0x72, 0x6e, - 0x65, 0x72, 0x5f, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, - 0x63, 0x6f, 0x72, 0x6e, 0x65, 0x72, 0x5f, 0x79, 0x2c, 0x0a, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x66, - 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x66, 0x73, 0x78, 0x2c, 0x20, 0x66, 0x6c, - 0x6f, 0x61, 0x74, 0x20, 0x66, 0x73, 0x79, 0x2c, 0x20, 0x66, 0x6c, 0x6f, - 0x61, 0x74, 0x20, 0x66, 0x73, 0x7a, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, - 0x74, 0x20, 0x73, 0x73, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, - 0x20, 0x73, 0x73, 0x79, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, - 0x73, 0x73, 0x7a, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, - 0x74, 0x20, 0x72, 0x78, 0x2c, 0x20, 0x72, 0x79, 0x2c, 0x20, 0x72, 0x7a, - 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x61, 0x7a, 0x2c, - 0x20, 0x74, 0x74, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, - 0x20, 0x71, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x78, - 0x73, 0x2c, 0x20, 0x79, 0x73, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, - 0x74, 0x20, 0x6b, 0x78, 0x2c, 0x20, 0x6b, 0x79, 0x2c, 0x20, 0x6b, 0x7a, - 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x74, 0x74, - 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x43, 0x61, 0x6c, 0x63, 0x75, - 0x6c, 0x61, 0x74, 0x65, 0x20, 0x33, 0x44, 0x20, 0x70, 0x6f, 0x73, 0x69, - 0x74, 0x69, 0x6f, 0x6e, 0x20, 0x6f, 0x66, 0x20, 0x67, 0x69, 0x76, 0x65, - 0x6e, 0x20, 0x70, 0x6f, 0x73, 0x69, 0x74, 0x69, 0x6f, 0x6e, 0x2c, 0x20, - 0x69, 0x6e, 0x20, 0x6d, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x72, 0x78, 0x20, - 0x3d, 0x20, 0x28, 0x63, 0x6f, 0x72, 0x6e, 0x65, 0x72, 0x5f, 0x78, 0x20, - 0x2b, 0x20, 0x66, 0x73, 0x2a, 0x66, 0x73, 0x78, 0x20, 0x2b, 0x20, 0x73, - 0x73, 0x2a, 0x73, 0x73, 0x78, 0x29, 0x20, 0x2f, 0x20, 0x72, 0x65, 0x73, - 0x3b, 0x0a, 0x09, 0x72, 0x79, 0x20, 0x3d, 0x20, 0x28, 0x63, 0x6f, 0x72, - 0x6e, 0x65, 0x72, 0x5f, 0x79, 0x20, 0x2b, 0x20, 0x66, 0x73, 0x2a, 0x66, - 0x73, 0x79, 0x20, 0x2b, 0x20, 0x73, 0x73, 0x2a, 0x73, 0x73, 0x79, 0x29, - 0x20, 0x2f, 0x20, 0x72, 0x65, 0x73, 0x3b, 0x0a, 0x09, 0x72, 0x7a, 0x20, - 0x3d, 0x20, 0x63, 0x6c, 0x65, 0x6e, 0x20, 0x2b, 0x20, 0x28, 0x66, 0x73, - 0x2a, 0x66, 0x73, 0x7a, 0x20, 0x2b, 0x20, 0x73, 0x73, 0x2a, 0x73, 0x73, - 0x7a, 0x29, 0x2f, 0x72, 0x65, 0x73, 0x3b, 0x0a, 0x0a, 0x09, 0x63, 0x74, - 0x74, 0x20, 0x3d, 0x20, 0x72, 0x7a, 0x20, 0x2f, 0x20, 0x73, 0x71, 0x72, - 0x74, 0x28, 0x72, 0x78, 0x2a, 0x72, 0x78, 0x20, 0x2b, 0x20, 0x72, 0x79, - 0x2a, 0x72, 0x79, 0x20, 0x2b, 0x20, 0x72, 0x7a, 0x2a, 0x72, 0x7a, 0x29, - 0x3b, 0x20, 0x20, 0x2f, 0x2a, 0x20, 0x63, 0x6f, 0x73, 0x28, 0x32, 0x74, - 0x68, 0x65, 0x74, 0x61, 0x29, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x74, 0x74, - 0x20, 0x3d, 0x20, 0x61, 0x63, 0x6f, 0x73, 0x28, 0x63, 0x74, 0x74, 0x29, - 0x3b, 0x0a, 0x09, 0x61, 0x7a, 0x20, 0x3d, 0x20, 0x61, 0x74, 0x61, 0x6e, - 0x32, 0x28, 0x72, 0x79, 0x2c, 0x20, 0x72, 0x78, 0x29, 0x3b, 0x0a, 0x0a, - 0x09, 0x6b, 0x78, 0x20, 0x3d, 0x20, 0x6b, 0x2a, 0x6e, 0x61, 0x74, 0x69, - 0x76, 0x65, 0x5f, 0x73, 0x69, 0x6e, 0x28, 0x74, 0x74, 0x29, 0x2a, 0x6e, - 0x61, 0x74, 0x69, 0x76, 0x65, 0x5f, 0x63, 0x6f, 0x73, 0x28, 0x61, 0x7a, - 0x29, 0x3b, 0x0a, 0x09, 0x6b, 0x79, 0x20, 0x3d, 0x20, 0x6b, 0x2a, 0x6e, - 0x61, 0x74, 0x69, 0x76, 0x65, 0x5f, 0x73, 0x69, 0x6e, 0x28, 0x74, 0x74, - 0x29, 0x2a, 0x6e, 0x61, 0x74, 0x69, 0x76, 0x65, 0x5f, 0x73, 0x69, 0x6e, - 0x28, 0x61, 0x7a, 0x29, 0x3b, 0x0a, 0x09, 0x6b, 0x7a, 0x20, 0x3d, 0x20, - 0x6b, 0x2a, 0x28, 0x63, 0x74, 0x74, 0x20, 0x2d, 0x20, 0x31, 0x2e, 0x30, - 0x29, 0x3b, 0x0a, 0x0a, 0x09, 0x71, 0x20, 0x3d, 0x20, 0x28, 0x66, 0x6c, - 0x6f, 0x61, 0x74, 0x34, 0x29, 0x28, 0x6b, 0x78, 0x2c, 0x20, 0x6b, 0x79, - 0x2c, 0x20, 0x6b, 0x7a, 0x2c, 0x20, 0x30, 0x2e, 0x30, 0x29, 0x3b, 0x0a, - 0x0a, 0x09, 0x72, 0x65, 0x74, 0x75, 0x72, 0x6e, 0x20, 0x71, 0x3b, 0x0a, - 0x7d, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x6c, 0x61, - 0x74, 0x74, 0x69, 0x63, 0x65, 0x5f, 0x66, 0x61, 0x63, 0x74, 0x6f, 0x72, - 0x28, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x31, 0x36, 0x20, 0x63, 0x65, 0x6c, - 0x6c, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, 0x20, 0x71, 0x2c, - 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x72, 0x65, - 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, 0x61, 0x67, - 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x61, - 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x72, - 0x65, 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, 0x61, - 0x67, 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, - 0x62, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x72, 0x65, 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, - 0x61, 0x67, 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, - 0x5f, 0x63, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, - 0x20, 0x66, 0x31, 0x2c, 0x20, 0x66, 0x32, 0x2c, 0x20, 0x66, 0x33, 0x2c, - 0x20, 0x76, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, 0x20, - 0x55, 0x64, 0x6f, 0x74, 0x71, 0x3b, 0x0a, 0x0a, 0x09, 0x55, 0x64, 0x6f, - 0x74, 0x71, 0x2e, 0x78, 0x20, 0x3d, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, - 0x73, 0x30, 0x2a, 0x71, 0x2e, 0x78, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, - 0x6c, 0x2e, 0x73, 0x31, 0x2a, 0x71, 0x2e, 0x79, 0x20, 0x2b, 0x20, 0x63, - 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x32, 0x2a, 0x71, 0x2e, 0x7a, 0x3b, 0x0a, - 0x09, 0x55, 0x64, 0x6f, 0x74, 0x71, 0x2e, 0x79, 0x20, 0x3d, 0x20, 0x63, - 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x33, 0x2a, 0x71, 0x2e, 0x78, 0x20, 0x2b, - 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x34, 0x2a, 0x71, 0x2e, 0x79, - 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x35, 0x2a, 0x71, - 0x2e, 0x7a, 0x3b, 0x0a, 0x09, 0x55, 0x64, 0x6f, 0x74, 0x71, 0x2e, 0x7a, - 0x20, 0x3d, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x36, 0x2a, 0x71, - 0x2e, 0x78, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x37, - 0x2a, 0x71, 0x2e, 0x79, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, - 0x73, 0x38, 0x2a, 0x71, 0x2e, 0x7a, 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, - 0x20, 0x4c, 0x6f, 0x6f, 0x6b, 0x20, 0x75, 0x70, 0x20, 0x76, 0x61, 0x6c, - 0x75, 0x65, 0x73, 0x20, 0x66, 0x72, 0x6f, 0x6d, 0x20, 0x70, 0x72, 0x65, - 0x63, 0x61, 0x6c, 0x63, 0x75, 0x6c, 0x61, 0x74, 0x65, 0x64, 0x20, 0x73, - 0x69, 0x6e, 0x63, 0x28, 0x29, 0x20, 0x74, 0x61, 0x62, 0x6c, 0x65, 0x20, - 0x2a, 0x2f, 0x0a, 0x09, 0x66, 0x31, 0x20, 0x3d, 0x20, 0x72, 0x65, 0x61, - 0x64, 0x5f, 0x69, 0x6d, 0x61, 0x67, 0x65, 0x66, 0x28, 0x66, 0x75, 0x6e, - 0x63, 0x5f, 0x61, 0x2c, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, - 0x5f, 0x61, 0x2c, 0x20, 0x28, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x32, 0x29, - 0x28, 0x55, 0x64, 0x6f, 0x74, 0x71, 0x2e, 0x78, 0x2c, 0x20, 0x30, 0x2e, - 0x30, 0x29, 0x29, 0x2e, 0x73, 0x30, 0x3b, 0x0a, 0x09, 0x66, 0x32, 0x20, - 0x3d, 0x20, 0x72, 0x65, 0x61, 0x64, 0x5f, 0x69, 0x6d, 0x61, 0x67, 0x65, - 0x66, 0x28, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x62, 0x2c, 0x20, 0x73, 0x61, - 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x62, 0x2c, 0x20, 0x28, 0x66, 0x6c, - 0x6f, 0x61, 0x74, 0x32, 0x29, 0x28, 0x55, 0x64, 0x6f, 0x74, 0x71, 0x2e, - 0x79, 0x2c, 0x20, 0x30, 0x2e, 0x30, 0x29, 0x29, 0x2e, 0x73, 0x30, 0x3b, - 0x0a, 0x09, 0x66, 0x33, 0x20, 0x3d, 0x20, 0x72, 0x65, 0x61, 0x64, 0x5f, - 0x69, 0x6d, 0x61, 0x67, 0x65, 0x66, 0x28, 0x66, 0x75, 0x6e, 0x63, 0x5f, - 0x63, 0x2c, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x72, 0x5f, 0x63, - 0x2c, 0x20, 0x28, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x32, 0x29, 0x28, 0x55, - 0x64, 0x6f, 0x74, 0x71, 0x2e, 0x7a, 0x2c, 0x20, 0x30, 0x2e, 0x30, 0x29, - 0x29, 0x2e, 0x73, 0x30, 0x3b, 0x0a, 0x0a, 0x09, 0x72, 0x65, 0x74, 0x75, - 0x72, 0x6e, 0x20, 0x66, 0x31, 0x20, 0x2a, 0x20, 0x66, 0x32, 0x20, 0x2a, - 0x20, 0x66, 0x33, 0x3b, 0x0a, 0x7d, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f, - 0x61, 0x74, 0x20, 0x6c, 0x6f, 0x6f, 0x6b, 0x75, 0x70, 0x5f, 0x69, 0x6e, - 0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x79, 0x28, 0x67, 0x6c, 0x6f, 0x62, - 0x61, 0x6c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x69, 0x6e, - 0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x69, 0x65, 0x73, 0x2c, 0x0a, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x73, 0x69, - 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x68, 0x2c, 0x20, - 0x73, 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x6b, - 0x2c, 0x20, 0x73, 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, - 0x20, 0x6c, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x69, 0x6e, 0x74, 0x20, 0x69, - 0x64, 0x78, 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x4f, 0x75, 0x74, - 0x20, 0x6f, 0x66, 0x20, 0x72, 0x61, 0x6e, 0x67, 0x65, 0x3f, 0x20, 0x2a, - 0x2f, 0x0a, 0x09, 0x69, 0x66, 0x20, 0x28, 0x20, 0x28, 0x61, 0x62, 0x73, - 0x28, 0x68, 0x29, 0x20, 0x3e, 0x20, 0x49, 0x4e, 0x44, 0x4d, 0x41, 0x58, - 0x29, 0x20, 0x7c, 0x7c, 0x20, 0x28, 0x61, 0x62, 0x73, 0x28, 0x6b, 0x29, - 0x20, 0x3e, 0x20, 0x49, 0x4e, 0x44, 0x4d, 0x41, 0x58, 0x29, 0x20, 0x7c, - 0x7c, 0x20, 0x28, 0x61, 0x62, 0x73, 0x28, 0x6c, 0x29, 0x20, 0x3e, 0x20, - 0x49, 0x4e, 0x44, 0x4d, 0x41, 0x58, 0x29, 0x20, 0x29, 0x20, 0x7b, 0x0a, - 0x09, 0x09, 0x72, 0x65, 0x74, 0x75, 0x72, 0x6e, 0x20, 0x30, 0x2e, 0x30, - 0x3b, 0x0a, 0x09, 0x7d, 0x0a, 0x0a, 0x09, 0x68, 0x20, 0x3d, 0x20, 0x28, - 0x68, 0x3e, 0x3d, 0x30, 0x29, 0x20, 0x3f, 0x20, 0x68, 0x20, 0x3a, 0x20, - 0x68, 0x2b, 0x49, 0x44, 0x49, 0x4d, 0x3b, 0x0a, 0x09, 0x6b, 0x20, 0x3d, - 0x20, 0x28, 0x6b, 0x3e, 0x3d, 0x30, 0x29, 0x20, 0x3f, 0x20, 0x6b, 0x20, - 0x3a, 0x20, 0x6b, 0x2b, 0x49, 0x44, 0x49, 0x4d, 0x3b, 0x0a, 0x09, 0x6c, - 0x20, 0x3d, 0x20, 0x28, 0x6c, 0x3e, 0x3d, 0x30, 0x29, 0x20, 0x3f, 0x20, - 0x6c, 0x20, 0x3a, 0x20, 0x6c, 0x2b, 0x49, 0x44, 0x49, 0x4d, 0x3b, 0x0a, - 0x0a, 0x09, 0x69, 0x64, 0x78, 0x20, 0x3d, 0x20, 0x68, 0x20, 0x2b, 0x20, - 0x28, 0x49, 0x44, 0x49, 0x4d, 0x2a, 0x6b, 0x29, 0x20, 0x2b, 0x20, 0x28, - 0x49, 0x44, 0x49, 0x4d, 0x2a, 0x49, 0x44, 0x49, 0x4d, 0x2a, 0x6c, 0x29, - 0x3b, 0x0a, 0x0a, 0x09, 0x72, 0x65, 0x74, 0x75, 0x72, 0x6e, 0x20, 0x69, - 0x6e, 0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x69, 0x65, 0x73, 0x5b, 0x69, - 0x64, 0x78, 0x5d, 0x3b, 0x0a, 0x7d, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f, - 0x61, 0x74, 0x20, 0x6c, 0x6f, 0x6f, 0x6b, 0x75, 0x70, 0x5f, 0x66, 0x6c, - 0x61, 0x67, 0x67, 0x65, 0x64, 0x5f, 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73, - 0x69, 0x74, 0x79, 0x28, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x20, 0x66, - 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73, - 0x69, 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x67, 0x6c, 0x6f, 0x62, 0x61, - 0x6c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x66, 0x6c, 0x61, - 0x67, 0x73, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x73, - 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x68, 0x2c, - 0x20, 0x73, 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, 0x74, 0x20, - 0x6b, 0x2c, 0x20, 0x73, 0x69, 0x67, 0x6e, 0x65, 0x64, 0x20, 0x69, 0x6e, - 0x74, 0x20, 0x6c, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x72, 0x65, 0x74, 0x75, - 0x72, 0x6e, 0x20, 0x6c, 0x6f, 0x6f, 0x6b, 0x75, 0x70, 0x5f, 0x69, 0x6e, - 0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x79, 0x28, 0x69, 0x6e, 0x74, 0x65, - 0x6e, 0x73, 0x69, 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x68, 0x2c, 0x20, - 0x6b, 0x2c, 0x20, 0x6c, 0x29, 0x0a, 0x09, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x2a, 0x20, 0x6c, 0x6f, 0x6f, 0x6b, 0x75, 0x70, 0x5f, - 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73, 0x69, 0x74, 0x79, 0x28, 0x66, 0x6c, - 0x61, 0x67, 0x73, 0x2c, 0x20, 0x68, 0x2c, 0x20, 0x6b, 0x2c, 0x20, 0x6c, - 0x29, 0x3b, 0x0a, 0x7d, 0x0a, 0x0a, 0x0a, 0x66, 0x6c, 0x6f, 0x61, 0x74, - 0x20, 0x6d, 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, 0x5f, 0x66, 0x61, - 0x63, 0x74, 0x6f, 0x72, 0x28, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x20, - 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x69, 0x6e, 0x74, 0x65, 0x6e, - 0x73, 0x69, 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x67, 0x6c, 0x6f, 0x62, - 0x61, 0x6c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x66, 0x6c, - 0x61, 0x67, 0x73, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x31, 0x36, 0x20, 0x63, - 0x65, 0x6c, 0x6c, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, 0x20, - 0x71, 0x29, 0x0a, 0x7b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, - 0x68, 0x66, 0x2c, 0x20, 0x6b, 0x66, 0x2c, 0x20, 0x6c, 0x66, 0x3b, 0x0a, - 0x09, 0x69, 0x6e, 0x74, 0x20, 0x68, 0x2c, 0x20, 0x6b, 0x2c, 0x20, 0x6c, - 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x76, 0x61, 0x6c, - 0x20, 0x3d, 0x20, 0x30, 0x2e, 0x30, 0x3b, 0x0a, 0x0a, 0x09, 0x23, 0x69, - 0x66, 0x64, 0x65, 0x66, 0x20, 0x46, 0x4c, 0x41, 0x54, 0x5f, 0x49, 0x4e, - 0x54, 0x45, 0x4e, 0x53, 0x49, 0x54, 0x49, 0x45, 0x53, 0x0a, 0x09, 0x72, - 0x65, 0x74, 0x75, 0x72, 0x6e, 0x20, 0x31, 0x30, 0x30, 0x2e, 0x30, 0x3b, - 0x0a, 0x09, 0x23, 0x65, 0x6c, 0x73, 0x65, 0x0a, 0x0a, 0x09, 0x68, 0x66, - 0x20, 0x3d, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x30, 0x2a, 0x71, - 0x2e, 0x78, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x31, - 0x2a, 0x71, 0x2e, 0x79, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, - 0x73, 0x32, 0x2a, 0x71, 0x2e, 0x7a, 0x3b, 0x20, 0x20, 0x2f, 0x2a, 0x20, - 0x68, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x6b, 0x66, 0x20, 0x3d, 0x20, 0x63, - 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x33, 0x2a, 0x71, 0x2e, 0x78, 0x20, 0x2b, - 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x34, 0x2a, 0x71, 0x2e, 0x79, - 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x35, 0x2a, 0x71, - 0x2e, 0x7a, 0x3b, 0x20, 0x20, 0x2f, 0x2a, 0x20, 0x6b, 0x20, 0x2a, 0x2f, - 0x0a, 0x09, 0x6c, 0x66, 0x20, 0x3d, 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2e, - 0x73, 0x36, 0x2a, 0x71, 0x2e, 0x78, 0x20, 0x2b, 0x20, 0x63, 0x65, 0x6c, - 0x6c, 0x2e, 0x73, 0x37, 0x2a, 0x71, 0x2e, 0x79, 0x20, 0x2b, 0x20, 0x63, - 0x65, 0x6c, 0x6c, 0x2e, 0x73, 0x38, 0x2a, 0x71, 0x2e, 0x7a, 0x3b, 0x20, - 0x20, 0x2f, 0x2a, 0x20, 0x6c, 0x20, 0x2a, 0x2f, 0x0a, 0x0a, 0x09, 0x68, - 0x20, 0x3d, 0x20, 0x72, 0x6f, 0x75, 0x6e, 0x64, 0x28, 0x68, 0x66, 0x29, - 0x3b, 0x0a, 0x09, 0x6b, 0x20, 0x3d, 0x20, 0x72, 0x6f, 0x75, 0x6e, 0x64, - 0x28, 0x6b, 0x66, 0x29, 0x3b, 0x0a, 0x09, 0x6c, 0x20, 0x3d, 0x20, 0x72, - 0x6f, 0x75, 0x6e, 0x64, 0x28, 0x6c, 0x66, 0x29, 0x3b, 0x0a, 0x0a, 0x09, - 0x2f, 0x2a, 0x20, 0x53, 0x79, 0x6d, 0x6d, 0x65, 0x74, 0x72, 0x79, 0x20, - 0x73, 0x74, 0x75, 0x66, 0x66, 0x20, 0x67, 0x6f, 0x65, 0x73, 0x20, 0x68, - 0x65, 0x72, 0x65, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x49, 0x4e, 0x53, 0x45, - 0x52, 0x54, 0x5f, 0x48, 0x45, 0x52, 0x45, 0x0a, 0x0a, 0x09, 0x72, 0x65, - 0x74, 0x75, 0x72, 0x6e, 0x20, 0x76, 0x61, 0x6c, 0x3b, 0x0a, 0x09, 0x23, - 0x65, 0x6e, 0x64, 0x69, 0x66, 0x20, 0x2f, 0x2a, 0x20, 0x46, 0x4c, 0x41, - 0x54, 0x5f, 0x49, 0x4e, 0x54, 0x45, 0x4e, 0x53, 0x49, 0x54, 0x49, 0x49, - 0x45, 0x53, 0x20, 0x2a, 0x2f, 0x0a, 0x7d, 0x0a, 0x0a, 0x0a, 0x6b, 0x65, - 0x72, 0x6e, 0x65, 0x6c, 0x20, 0x76, 0x6f, 0x69, 0x64, 0x20, 0x64, 0x69, - 0x66, 0x66, 0x72, 0x61, 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x28, 0x67, 0x6c, - 0x6f, 0x62, 0x61, 0x6c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, - 0x64, 0x69, 0x66, 0x66, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, - 0x6b, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x77, 0x65, 0x69, - 0x67, 0x68, 0x74, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x77, 0x2c, 0x20, - 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6f, 0x72, 0x6e, 0x65, 0x72, - 0x5f, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6f, - 0x72, 0x6e, 0x65, 0x72, 0x5f, 0x79, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, - 0x74, 0x20, 0x66, 0x73, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, - 0x20, 0x66, 0x73, 0x79, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, - 0x66, 0x73, 0x7a, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x73, - 0x73, 0x78, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x73, 0x73, - 0x79, 0x2c, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x73, 0x73, 0x7a, - 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x72, 0x65, 0x73, 0x2c, - 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x63, 0x6c, 0x65, 0x6e, 0x2c, - 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x31, 0x36, 0x20, 0x63, 0x65, 0x6c, - 0x6c, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x20, 0x66, 0x6c, - 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73, 0x69, - 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, - 0x20, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x2a, 0x66, 0x6c, 0x61, 0x67, - 0x73, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x72, 0x65, 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, - 0x20, 0x69, 0x6d, 0x61, 0x67, 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, - 0x75, 0x6e, 0x63, 0x5f, 0x61, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x72, 0x65, 0x61, 0x64, 0x5f, - 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, 0x61, 0x67, 0x65, 0x32, 0x64, - 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x62, 0x2c, 0x0a, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x72, - 0x65, 0x61, 0x64, 0x5f, 0x6f, 0x6e, 0x6c, 0x79, 0x20, 0x69, 0x6d, 0x61, - 0x67, 0x65, 0x32, 0x64, 0x5f, 0x74, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, - 0x63, 0x2c, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x6c, 0x6f, 0x63, 0x61, 0x6c, 0x20, 0x66, 0x6c, 0x6f, - 0x61, 0x74, 0x20, 0x2a, 0x74, 0x6d, 0x70, 0x29, 0x0a, 0x7b, 0x0a, 0x09, - 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x66, 0x73, 0x2c, 0x20, 0x73, 0x73, - 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x66, 0x5f, 0x6c, - 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x2c, 0x20, 0x49, 0x5f, 0x6c, 0x61, - 0x74, 0x74, 0x69, 0x63, 0x65, 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, - 0x74, 0x20, 0x49, 0x5f, 0x6d, 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, - 0x3b, 0x0a, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x34, 0x20, 0x71, 0x3b, - 0x0a, 0x09, 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20, - 0x6c, 0x73, 0x30, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x6c, 0x6f, - 0x63, 0x61, 0x6c, 0x5f, 0x73, 0x69, 0x7a, 0x65, 0x28, 0x30, 0x29, 0x3b, - 0x0a, 0x09, 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20, - 0x6c, 0x73, 0x31, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x6c, 0x6f, - 0x63, 0x61, 0x6c, 0x5f, 0x73, 0x69, 0x7a, 0x65, 0x28, 0x31, 0x29, 0x3b, - 0x0a, 0x09, 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20, - 0x6c, 0x69, 0x30, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x6c, 0x6f, - 0x63, 0x61, 0x6c, 0x5f, 0x69, 0x64, 0x28, 0x30, 0x29, 0x3b, 0x0a, 0x09, - 0x63, 0x6f, 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x6c, 0x69, - 0x31, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x6c, 0x6f, 0x63, 0x61, - 0x6c, 0x5f, 0x69, 0x64, 0x28, 0x31, 0x29, 0x3b, 0x0a, 0x09, 0x63, 0x6f, - 0x6e, 0x73, 0x74, 0x20, 0x69, 0x6e, 0x74, 0x20, 0x6c, 0x73, 0x20, 0x3d, - 0x20, 0x6c, 0x73, 0x30, 0x20, 0x2a, 0x20, 0x6c, 0x73, 0x31, 0x3b, 0x0a, - 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x43, 0x61, 0x6c, 0x63, 0x75, 0x6c, 0x61, - 0x74, 0x65, 0x20, 0x66, 0x72, 0x61, 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x61, - 0x6c, 0x20, 0x63, 0x6f, 0x6f, 0x72, 0x64, 0x69, 0x6e, 0x61, 0x74, 0x65, - 0x73, 0x20, 0x69, 0x6e, 0x20, 0x66, 0x73, 0x2f, 0x73, 0x73, 0x20, 0x2a, - 0x2f, 0x0a, 0x09, 0x66, 0x73, 0x20, 0x3d, 0x20, 0x63, 0x6f, 0x6e, 0x76, - 0x65, 0x72, 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x67, 0x65, - 0x74, 0x5f, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x5f, 0x69, 0x64, 0x28, - 0x30, 0x29, 0x29, 0x20, 0x2f, 0x20, 0x63, 0x6f, 0x6e, 0x76, 0x65, 0x72, - 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x6c, 0x73, 0x30, 0x29, - 0x3b, 0x0a, 0x09, 0x73, 0x73, 0x20, 0x3d, 0x20, 0x63, 0x6f, 0x6e, 0x76, - 0x65, 0x72, 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x67, 0x65, - 0x74, 0x5f, 0x67, 0x6c, 0x6f, 0x62, 0x61, 0x6c, 0x5f, 0x69, 0x64, 0x28, - 0x31, 0x29, 0x29, 0x20, 0x2f, 0x20, 0x63, 0x6f, 0x6e, 0x76, 0x65, 0x72, - 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x6c, 0x73, 0x31, 0x29, - 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x47, 0x65, 0x74, 0x20, 0x74, - 0x68, 0x65, 0x20, 0x73, 0x63, 0x61, 0x74, 0x74, 0x65, 0x72, 0x69, 0x6e, - 0x67, 0x20, 0x76, 0x65, 0x63, 0x74, 0x6f, 0x72, 0x20, 0x2a, 0x2f, 0x0a, - 0x09, 0x71, 0x20, 0x3d, 0x20, 0x67, 0x65, 0x74, 0x5f, 0x71, 0x28, 0x66, - 0x73, 0x2c, 0x20, 0x73, 0x73, 0x2c, 0x20, 0x72, 0x65, 0x73, 0x2c, 0x20, - 0x63, 0x6c, 0x65, 0x6e, 0x2c, 0x20, 0x6b, 0x2c, 0x0a, 0x09, 0x20, 0x20, - 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x20, 0x63, 0x6f, 0x72, 0x6e, - 0x65, 0x72, 0x5f, 0x78, 0x2c, 0x20, 0x63, 0x6f, 0x72, 0x6e, 0x65, 0x72, - 0x5f, 0x79, 0x2c, 0x20, 0x66, 0x73, 0x78, 0x2c, 0x20, 0x66, 0x73, 0x79, - 0x2c, 0x20, 0x66, 0x73, 0x7a, 0x2c, 0x20, 0x73, 0x73, 0x78, 0x2c, 0x20, - 0x73, 0x73, 0x79, 0x2c, 0x20, 0x73, 0x73, 0x7a, 0x29, 0x3b, 0x0a, 0x0a, - 0x09, 0x2f, 0x2a, 0x20, 0x43, 0x61, 0x6c, 0x63, 0x75, 0x6c, 0x61, 0x74, - 0x65, 0x20, 0x74, 0x68, 0x65, 0x20, 0x64, 0x69, 0x66, 0x66, 0x72, 0x61, - 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x66, 0x5f, - 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x20, 0x3d, 0x20, 0x6c, 0x61, - 0x74, 0x74, 0x69, 0x63, 0x65, 0x5f, 0x66, 0x61, 0x63, 0x74, 0x6f, 0x72, - 0x28, 0x63, 0x65, 0x6c, 0x6c, 0x2c, 0x20, 0x71, 0x2c, 0x20, 0x66, 0x75, - 0x6e, 0x63, 0x5f, 0x61, 0x2c, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x62, - 0x2c, 0x20, 0x66, 0x75, 0x6e, 0x63, 0x5f, 0x63, 0x29, 0x3b, 0x0a, 0x09, - 0x49, 0x5f, 0x6d, 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, 0x20, 0x3d, - 0x20, 0x6d, 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, 0x5f, 0x66, 0x61, - 0x63, 0x74, 0x6f, 0x72, 0x28, 0x69, 0x6e, 0x74, 0x65, 0x6e, 0x73, 0x69, - 0x74, 0x69, 0x65, 0x73, 0x2c, 0x20, 0x66, 0x6c, 0x61, 0x67, 0x73, 0x2c, - 0x20, 0x63, 0x65, 0x6c, 0x6c, 0x2c, 0x20, 0x71, 0x29, 0x3b, 0x0a, 0x09, - 0x49, 0x5f, 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x20, 0x3d, 0x20, - 0x70, 0x6f, 0x77, 0x28, 0x66, 0x5f, 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, - 0x65, 0x2c, 0x20, 0x32, 0x2e, 0x30, 0x66, 0x29, 0x3b, 0x0a, 0x0a, 0x09, - 0x74, 0x6d, 0x70, 0x5b, 0x6c, 0x69, 0x30, 0x20, 0x2b, 0x20, 0x6c, 0x73, - 0x30, 0x2a, 0x6c, 0x69, 0x31, 0x5d, 0x20, 0x3d, 0x20, 0x49, 0x5f, 0x6d, - 0x6f, 0x6c, 0x65, 0x63, 0x75, 0x6c, 0x65, 0x20, 0x2a, 0x20, 0x49, 0x5f, - 0x6c, 0x61, 0x74, 0x74, 0x69, 0x63, 0x65, 0x3b, 0x0a, 0x0a, 0x09, 0x62, - 0x61, 0x72, 0x72, 0x69, 0x65, 0x72, 0x28, 0x43, 0x4c, 0x4b, 0x5f, 0x4c, - 0x4f, 0x43, 0x41, 0x4c, 0x5f, 0x4d, 0x45, 0x4d, 0x5f, 0x46, 0x45, 0x4e, - 0x43, 0x45, 0x29, 0x3b, 0x0a, 0x0a, 0x09, 0x2f, 0x2a, 0x20, 0x46, 0x69, - 0x72, 0x73, 0x74, 0x20, 0x74, 0x68, 0x72, 0x65, 0x61, 0x64, 0x20, 0x69, - 0x6e, 0x20, 0x67, 0x72, 0x6f, 0x75, 0x70, 0x20, 0x73, 0x75, 0x6d, 0x73, - 0x20, 0x74, 0x68, 0x65, 0x20, 0x73, 0x61, 0x6d, 0x70, 0x6c, 0x65, 0x73, - 0x20, 0x2a, 0x2f, 0x0a, 0x09, 0x69, 0x66, 0x20, 0x28, 0x20, 0x6c, 0x69, - 0x30, 0x20, 0x2b, 0x20, 0x6c, 0x69, 0x31, 0x20, 0x3d, 0x3d, 0x20, 0x30, - 0x20, 0x29, 0x20, 0x7b, 0x0a, 0x0a, 0x09, 0x09, 0x69, 0x6e, 0x74, 0x20, - 0x69, 0x3b, 0x0a, 0x09, 0x09, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x73, - 0x75, 0x6d, 0x20, 0x3d, 0x20, 0x30, 0x2e, 0x30, 0x3b, 0x0a, 0x09, 0x09, - 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x20, 0x76, 0x61, 0x6c, 0x3b, 0x0a, 0x09, - 0x09, 0x69, 0x6e, 0x74, 0x20, 0x69, 0x64, 0x78, 0x3b, 0x0a, 0x0a, 0x09, - 0x09, 0x69, 0x64, 0x78, 0x20, 0x3d, 0x20, 0x63, 0x6f, 0x6e, 0x76, 0x65, - 0x72, 0x74, 0x5f, 0x69, 0x6e, 0x74, 0x5f, 0x72, 0x74, 0x7a, 0x28, 0x66, - 0x73, 0x29, 0x20, 0x2b, 0x20, 0x77, 0x2a, 0x63, 0x6f, 0x6e, 0x76, 0x65, - 0x72, 0x74, 0x5f, 0x69, 0x6e, 0x74, 0x5f, 0x72, 0x74, 0x7a, 0x28, 0x73, - 0x73, 0x29, 0x3b, 0x0a, 0x0a, 0x09, 0x09, 0x66, 0x6f, 0x72, 0x20, 0x28, - 0x20, 0x69, 0x3d, 0x30, 0x3b, 0x20, 0x69, 0x3c, 0x6c, 0x73, 0x3b, 0x20, - 0x69, 0x2b, 0x2b, 0x20, 0x29, 0x20, 0x73, 0x75, 0x6d, 0x20, 0x2b, 0x3d, - 0x20, 0x74, 0x6d, 0x70, 0x5b, 0x69, 0x5d, 0x3b, 0x0a, 0x0a, 0x09, 0x09, - 0x76, 0x61, 0x6c, 0x20, 0x3d, 0x20, 0x77, 0x65, 0x69, 0x67, 0x68, 0x74, - 0x20, 0x2a, 0x20, 0x73, 0x75, 0x6d, 0x20, 0x2f, 0x20, 0x63, 0x6f, 0x6e, - 0x76, 0x65, 0x72, 0x74, 0x5f, 0x66, 0x6c, 0x6f, 0x61, 0x74, 0x28, 0x6c, - 0x73, 0x29, 0x3b, 0x0a, 0x09, 0x09, 0x64, 0x69, 0x66, 0x66, 0x5b, 0x69, - 0x64, 0x78, 0x5d, 0x20, 0x3d, 0x20, 0x76, 0x61, 0x6c, 0x3b, 0x0a, 0x0a, - 0x09, 0x7d, 0x0a, 0x0a, 0x7d, 0x0a -}; -unsigned int data_diffraction_cl_len = 6006; diff --git a/src/diffraction.h b/src/diffraction.h deleted file mode 100644 index 63b022c9..00000000 --- a/src/diffraction.h +++ /dev/null @@ -1,66 +0,0 @@ -/* - * diffraction.h - * - * Calculate diffraction patterns by Fourier methods - * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2009-2019 Thomas White <taw@physics.org> - * 2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de> - * 2013 Alexandra Tolstikova - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#ifndef DIFFRACTION_H -#define DIFFRACTION_H - -#include <gsl/gsl_rng.h> - -#include "image.h" -#include "cell.h" -#include "symmetry.h" - - -typedef enum { - GRADIENT_MOSAIC, - GRADIENT_INTERPOLATE, - GRADIENT_PHASED -} GradientMethod; - -extern void get_diffraction(struct image *image, int na, int nb, int nc, - const double *intensities, const double *phases, - const unsigned char *flags, UnitCell *cell, - GradientMethod m, const SymOpList *sym, - int no_fringes, int flat, int n_samples); - -extern struct sample *generate_tophat(struct image *image); - -extern struct sample *generate_SASE(struct image *image, gsl_rng *rng); - -extern struct sample *generate_twocolour(struct image *image); - -extern struct sample *generate_spectrum_fromfile(struct image *image, - char *spectrum_fn); - -#endif /* DIFFRACTION_H */ diff --git a/src/list_tmp.h b/src/list_tmp.h deleted file mode 100644 index 0c1a84e4..00000000 --- a/src/list_tmp.h +++ /dev/null @@ -1,104 +0,0 @@ -/* - * Template for creating indexed 3D lists of a given type, usually indexed - * as signed h,k,l values where -INDMAX<={h,k,l}<=+INDMAX. - * - * These are used, for example, for: - * - a list of 'double complex' values for storing structure factors, - * - a list of 'double' values for storing reflection intensities, - * - a list of 'unsigned int' values for counts of some sort. - * - * When LABEL and TYPE are #defined appropriately, including this header - * defines functions such as: - * - new_list_<LABEL>(), for creating a new list of the given type, - * - set_<LABEL>(), for setting a value in a list, - * - lookup_<LABEL>(), for retrieving values from a list, - * ..and so on. - * - * See src/utils.h for more information. - * - */ - -/* - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2009-2012 Thomas White <taw@physics.org> - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - -#include <stdio.h> - -#define ERROR_T(...) fprintf(stderr, __VA_ARGS__) - - -static inline void LABEL(set_arr)(TYPE *ref, signed int h, - signed int k, signed int l, - TYPE i) -{ - int idx; - - if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) { - ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l); - ERROR_T("You need to re-configure INDMAX and re-run.\n"); - exit(1); - } - - if ( h < 0 ) h += IDIM; - if ( k < 0 ) k += IDIM; - if ( l < 0 ) l += IDIM; - - idx = h + (IDIM*k) + (IDIM*IDIM*l); - ref[idx] = i; -} - - -static inline TYPE LABEL(lookup_arr)(const TYPE *ref, signed int h, - signed int k, signed int l) -{ - int idx; - - if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) { - ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l); - ERROR_T("You need to re-configure INDMAX and re-run.\n"); - exit(1); - } - - if ( h < 0 ) h += IDIM; - if ( k < 0 ) k += IDIM; - if ( l < 0 ) l += IDIM; - - idx = h + (IDIM*k) + (IDIM*IDIM*l); - return ref[idx]; -} - - -static inline TYPE *LABEL(new_arr)(void) -{ - TYPE *r; - size_t r_size; - r_size = IDIM*IDIM*IDIM*sizeof(TYPE); - r = malloc(r_size); - memset(r, 0, r_size); - return r; -} - - -#undef LABEL -#undef TYPE -#undef ERROR_T diff --git a/src/partial_sim.c b/src/partial_sim.c deleted file mode 100644 index f5588ceb..00000000 --- a/src/partial_sim.c +++ /dev/null @@ -1,1034 +0,0 @@ -/* - * partial_sim.c - * - * Generate partials for testing scaling - * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2011-2020 Thomas White <taw@physics.org> - * 2014 Valerio Mariani - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <stdarg.h> -#include <stdlib.h> -#include <stdio.h> -#include <string.h> -#include <unistd.h> -#include <getopt.h> -#include <assert.h> -#include <pthread.h> -#include <gsl/gsl_rng.h> - -#include <image.h> -#include <utils.h> -#include <reflist-utils.h> -#include <symmetry.h> -#include <geometry.h> -#include <stream.h> -#include <thread-pool.h> -#include <cell-utils.h> - -#include "version.h" - - -/* Number of bins for partiality graph */ -#define NBINS 50 - - -static void mess_up_cell(Crystal *cr, double cnoise, gsl_rng *rng) -{ - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - UnitCell *cell = crystal_get_cell(cr); - - //STATUS("Real:\n"); - //cell_print(cell); - - cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - ax = flat_noise(rng, ax, cnoise*fabs(ax)/100.0); - ay = flat_noise(rng, ay, cnoise*fabs(ay)/100.0); - az = flat_noise(rng, az, cnoise*fabs(az)/100.0); - bx = flat_noise(rng, bx, cnoise*fabs(bx)/100.0); - by = flat_noise(rng, by, cnoise*fabs(by)/100.0); - bz = flat_noise(rng, bz, cnoise*fabs(bz)/100.0); - cx = flat_noise(rng, cx, cnoise*fabs(cx)/100.0); - cy = flat_noise(rng, cy, cnoise*fabs(cy)/100.0); - cz = flat_noise(rng, cz, cnoise*fabs(cz)/100.0); - cell_set_reciprocal(cell, ax, ay, az, bx, by, bz, cx, cy, cz); - - //STATUS("Changed:\n"); - //cell_print(cell); -} - - -/* For each reflection in "partial", fill in what the intensity would be - * according to "full" */ -static void calculate_partials(Crystal *cr, - RefList *full, const SymOpList *sym, - int random_intensities, - pthread_rwlock_t *full_lock, - unsigned long int *n_ref, double *p_hist, - double *p_max, double max_q, double full_stddev, - double noise_stddev, gsl_rng *rng, - UnitCell *template_cell, RefList *template_reflist) -{ - Reflection *refl; - RefListIterator *iter; - - for ( refl = first_refl(crystal_get_reflections(cr), &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int h, k, l; - Reflection *rfull; - double L, p, Ip, If; - int bin; - double res; - - get_indices(refl, &h, &k, &l); - get_asymm(sym, h, k, l, &h, &k, &l); - p = get_partiality(refl); - L = get_lorentz(refl); - - pthread_rwlock_rdlock(full_lock); - rfull = find_refl(full, h, k, l); - pthread_rwlock_unlock(full_lock); - - if ( rfull == NULL ) { - if ( random_intensities ) { - - pthread_rwlock_wrlock(full_lock); - - /* In the gap between the unlock and the wrlock, - * the reflection might have been created by - * another thread. So, we must check again */ - rfull = find_refl(full, h, k, l); - if ( rfull == NULL ) { - rfull = add_refl(full, h, k, l); - If = fabs(gaussian_noise(rng, 0.0, - full_stddev)); - set_intensity(rfull, If); - set_redundancy(rfull, 1); - } else { - If = get_intensity(rfull); - } - pthread_rwlock_unlock(full_lock); - - } else { - set_redundancy(refl, 0); - If = 0.0; - } - } else { - If = get_intensity(rfull); - if ( random_intensities ) { - lock_reflection(rfull); - int red = get_redundancy(rfull); - set_redundancy(rfull, red+1); - unlock_reflection(rfull); - } - } - - Ip = crystal_get_osf(cr) * L * p * If; - - res = resolution(crystal_get_cell(cr), h, k, l); - bin = NBINS*2.0*res/max_q; - if ( (bin < NBINS) && (bin>=0) ) { - p_hist[bin] += p; - n_ref[bin]++; - if ( p > p_max[bin] ) p_max[bin] = p; - } else { - STATUS("Reflection out of histogram range: %e %i %f\n", - res, bin, p); - } - - Ip = gaussian_noise(rng, Ip, noise_stddev); - - set_intensity(refl, Ip); - set_esd_intensity(refl, noise_stddev); - } -} - - -static void draw_and_write_image(struct image *image, - const DataTemplate *dtempl, - RefList *reflections, - gsl_rng *rng, double background) -{ - Reflection *refl; - RefListIterator *iter; - int i; - - image->dp = malloc(image->detgeom->n_panels*sizeof(float *)); - if ( image->dp == NULL ) { - ERROR("Failed to allocate data\n"); - return; - } - for ( i=0; i<image->detgeom->n_panels; i++ ) { - - int j; - struct detgeom_panel *p = &image->detgeom->panels[i]; - - image->dp[i] = calloc(p->w * p->h, sizeof(float)); - if ( image->dp[i] == NULL ) { - ERROR("Failed to allocate data\n"); - return; - } - for ( j=0; j<p->w*p->h; j++ ) { - image->dp[i][j] = poisson_noise(rng, background); - } - - } - - for ( refl = first_refl(reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double Ip; - double dfs, dss; - int fs, ss; - struct detgeom_panel *p; - signed int pn; - - Ip = get_intensity(refl); - - get_detector_pos(refl, &dfs, &dss); - pn = get_panel_number(refl); - assert(pn < image->detgeom->n_panels); - p = &image->detgeom->panels[pn]; - - /* Explicit rounding, downwards */ - fs = dfs; ss = dss; - assert(fs >= 0); - assert(ss >= 0); - assert(fs < p->w); - assert(ss < p->h); - - image->dp[pn][fs + p->w*ss] += Ip; - } - - image_write(image, dtempl, image->filename); - - for ( i=0; i<image->detgeom->n_panels; i++ ) { - free(image->dp[i]); - } - free(image->dp); -} - - -static void show_help(const char *s) -{ - printf("Syntax: %s [options]\n\n", s); - printf( -"Generate a stream containing partials from a reflection list.\n" -"\n" -" -h, --help Display this help message.\n" -" --version Print CrystFEL version number and exit.\n" -"\n" -"You need to provide the following basic options:\n" -" -i, --input=<file> Read reflections from <file>.\n" -" Default: generate random ones instead (see -r).\n" -" -o, --output=<file> Write partials in stream format to <file>.\n" -" --images=<prefix> Write images to <prefix>NNN.h5.\n" -" -g. --geometry=<file> Get detector geometry from file.\n" -" -p, --pdb=<file> PDB file from which to get the unit cell.\n" -"\n" -" -y, --symmetry=<sym> Symmetry of the input reflection list.\n" -" -n <n> Simulate <n> patterns. Default: 2.\n" -" -r, --save-random=<file> Save randomly generated intensities to file.\n" -" --pgraph=<file> Save a histogram of partiality values to file.\n" -" -c, --cnoise=<val> Amount of reciprocal space cell noise, in percent.\n" -" --osf-stddev=<val> Standard deviation of the scaling factors.\n" -" --full-stddev=<val> Standard deviation of the randomly\n" -" generated full intensities, if not using -i.\n" -" --noise-stddev=<val> Set the standard deviation of the noise.\n" -" --background=<val> Background level in photons. Default 3000.\n" -" --profile-radius Reciprocal space reflection profile radius in m^-1.\n" -" Default 0.001e9 m^-1\n" -" --really-random Be non-deterministic.\n" -"\n" -); -} - - -struct partial_sim_queue_args -{ - RefList *full; - pthread_rwlock_t full_lock; - const DataTemplate *dtempl; - - int n_done; - int n_started; - int n_to_do; - - SymOpList *sym; - int random_intensities; - UnitCell *cell; - double cnoise; - double osf_stddev; - double full_stddev; - double noise_stddev; - double background; - double profile_radius; - - double max_q; - - char *image_prefix; - - /* The overall histogram */ - double p_hist[NBINS]; - unsigned long int n_ref[NBINS]; - double p_max[NBINS]; - - Stream *stream; - gsl_rng **rngs; - Stream *template_stream; -}; - - -struct partial_sim_worker_args -{ - struct partial_sim_queue_args *qargs; - struct image *image; - - UnitCell *template_cell; - RefList *template_reflist; - - /* Histogram for this image */ - double p_hist[NBINS]; - unsigned long int n_ref[NBINS]; - double p_max[NBINS]; - - int n; -}; - - -static void *create_job(void *vqargs) -{ - struct partial_sim_worker_args *wargs; - struct partial_sim_queue_args *qargs = vqargs; - - /* All done already? */ - if ( qargs->n_started == qargs->n_to_do ) return NULL; - - wargs = malloc(sizeof(struct partial_sim_worker_args)); - - wargs->qargs = qargs; - - if ( qargs->template_stream != NULL ) { - - struct image *image; - - image = stream_read_chunk(qargs->template_stream, - STREAM_REFLECTIONS); - if ( image == NULL ) { - ERROR("Failed to read template chunk!\n"); - free(wargs); - return NULL; - } - if ( image->n_crystals != 1 ) { - ERROR("Template stream must have exactly one crystal " - "per frame.\n"); - free(wargs); - return NULL; - } - - wargs->template_cell = crystal_get_cell(image->crystals[0]); - wargs->template_reflist = crystal_get_reflections(image->crystals[0]); - - image_free(image); - - } else { - wargs->template_cell = NULL; - wargs->template_reflist = NULL; - } - - qargs->n_started++; - wargs->n = qargs->n_started; - - return wargs; -} - - -static void run_job(void *vwargs, int cookie) -{ - struct partial_sim_worker_args *wargs = vwargs; - struct partial_sim_queue_args *qargs = wargs->qargs; - int i; - Crystal *cr; - double osf; - struct image *image; - - image = image_create_for_simulation(qargs->dtempl); - if ( image == NULL ) { - ERROR("Failed to create image.\n"); - return; - } - - cr = crystal_new(); - if ( cr == NULL ) { - ERROR("Failed to create crystal.\n"); - return; - } - crystal_set_image(cr, image); - image_add_crystal(image, cr); - - do { - osf = gaussian_noise(qargs->rngs[cookie], 1.0, - qargs->osf_stddev); - } while ( osf <= 0.0 ); - crystal_set_osf(cr, osf); - crystal_set_mosaicity(cr, 0.0); - crystal_set_profile_radius(cr, qargs->profile_radius); - - if ( wargs->template_cell == NULL ) { - /* Set up a random orientation */ - struct quaternion orientation; - char tmp[128]; - orientation = random_quaternion(qargs->rngs[cookie]); - crystal_set_cell(cr, cell_rotate(qargs->cell, orientation)); - snprintf(tmp, 127, "quaternion = %f %f %f %f", - orientation.w, orientation.x, orientation.y, orientation.z); - crystal_add_notes(cr, tmp); - } else { - crystal_set_cell(cr, wargs->template_cell); - } - - image->filename = malloc(256); - if ( image->filename == NULL ) { - ERROR("Failed to allocate filename.\n"); - return; - } - if ( qargs->image_prefix != NULL ) { - snprintf(image->filename, 255, "%s%i.h5", - qargs->image_prefix, wargs->n); - } else { - snprintf(image->filename, 255, "dummy-%i.h5", wargs->n); - } - - if ( wargs->template_reflist == NULL ) { - RefList *reflections; - reflections = predict_to_res(cr, qargs->max_q); - crystal_set_reflections(cr, reflections); - calculate_partialities(cr, PMODEL_XSPHERE); - } else { - crystal_set_reflections(cr, wargs->template_reflist); - update_predictions(cr); - calculate_partialities(cr, PMODEL_XSPHERE); - } - - for ( i=0; i<NBINS; i++ ) { - wargs->n_ref[i] = 0; - wargs->p_hist[i] = 0.0; - wargs->p_max[i] = 0.0; - } - - calculate_partials(cr, qargs->full, - qargs->sym, qargs->random_intensities, - &qargs->full_lock, - wargs->n_ref, wargs->p_hist, wargs->p_max, - qargs->max_q, qargs->full_stddev, - qargs->noise_stddev, qargs->rngs[cookie], - wargs->template_cell, wargs->template_reflist); - - if ( qargs->image_prefix != NULL ) { - draw_and_write_image(image, qargs->dtempl, - crystal_get_reflections(cr), - qargs->rngs[cookie], - qargs->background); - } - - /* Give a slightly incorrect cell in the stream */ - mess_up_cell(cr, qargs->cnoise, qargs->rngs[cookie]); - - wargs->image = image; -} - - -static void finalise_job(void *vqargs, void *vwargs) -{ - struct partial_sim_worker_args *wargs = vwargs; - struct partial_sim_queue_args *qargs = vqargs; - int i; - int ret; - - ret = stream_write_chunk(qargs->stream, wargs->image, - STREAM_REFLECTIONS); - if ( ret != 0 ) { - ERROR("WARNING: error writing stream file.\n"); - } - - for ( i=0; i<NBINS; i++ ) { - qargs->n_ref[i] += wargs->n_ref[i]; - qargs->p_hist[i] += wargs->p_hist[i]; - if ( wargs->p_max[i] > qargs->p_max[i] ) { - qargs->p_max[i] = wargs->p_max[i]; - } - } - - qargs->n_done++; - progress_bar(qargs->n_done, qargs->n_to_do, "Simulating"); - - image_free(wargs->image); -} - - -int main(int argc, char *argv[]) -{ - int c; - char *input_file = NULL; - char *output_file = NULL; - char *geomfile = NULL; - char *cellfile = NULL; - DataTemplate *dtempl; - RefList *full = NULL; - char *sym_str = NULL; - SymOpList *sym; - UnitCell *cell = NULL; - Stream *stream; - int n = 2; - int random_intensities = 0; - char *save_file = NULL; - struct partial_sim_queue_args qargs; - int n_threads = 1; - char *rval; - int i; - char *phist_file = NULL; - int config_random = 0; - char *image_prefix = NULL; - Stream *template_stream = NULL; - char *template = NULL; - struct image *test_image; - - /* Default simulation parameters */ - double profile_radius = 0.001e9; - double osf_stddev = 2.0; - double full_stddev = 1000.0; - double noise_stddev = 20.0; - double background = 3000.0; - double cnoise = 0.0; - - /* Long options */ - const struct option longopts[] = { - {"help", 0, NULL, 'h'}, - {"version", 0, NULL, 'v'}, - {"beam", 1, NULL, 'b'}, - {"output", 1, NULL, 'o'}, - {"input", 1, NULL, 'i'}, - {"pdb", 1, NULL, 'p'}, - {"geometry", 1, NULL, 'g'}, - {"symmetry", 1, NULL, 'y'}, - {"save-random", 1, NULL, 'r'}, - {"cnoise", 1, NULL, 'c'}, - - {"pgraph", 1, NULL, 2}, - {"osf-stddev", 1, NULL, 3}, - {"full-stddev", 1, NULL, 4}, - {"noise-stddev", 1, NULL, 5}, - {"images", 1, NULL, 6}, - {"background", 1, NULL, 7}, - {"beam-divergence", 1, NULL, 8}, - {"beam-bandwidth", 1, NULL, 9}, - {"profile-radius", 1, NULL, 10}, - {"photon-energy", 1, NULL, 11}, - {"template-stream", 1, NULL, 12}, - - {"really-random", 0, &config_random, 1}, - - {0, 0, NULL, 0} - }; - - /* Short options */ - while ((c = getopt_long(argc, argv, "hi:o:p:g:y:n:r:j:c:vb:", - longopts, NULL)) != -1) - { - switch (c) { - - case 'h' : - show_help(argv[0]); - return 0; - - case 'v' : - printf("CrystFEL: %s\n", - crystfel_version_string()); - printf("%s\n", - crystfel_licence_string()); - return 0; - - case 'b' : - ERROR("WARNING: This version of CrystFEL no longer " - "uses beam files. Please remove the beam file " - "from your partial_sim command line.\n"); - return 1; - - case 'i' : - input_file = strdup(optarg); - break; - - case 'o' : - output_file = strdup(optarg); - break; - - case 'p' : - cellfile = strdup(optarg); - break; - - case 'g' : - geomfile = strdup(optarg); - break; - - case 'y' : - sym_str = strdup(optarg); - break; - - case 'n' : - n = atoi(optarg); - break; - - case 'r' : - save_file = strdup(optarg); - break; - - case 'j' : - n_threads = atoi(optarg); - break; - - case 'c' : - cnoise = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid cell noise value.\n"); - return 1; - } - break; - - case 2 : - phist_file = strdup(optarg); - break; - - case 3 : - osf_stddev = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid OSF standard deviation.\n"); - return 1; - } - if ( osf_stddev < 0.0 ) { - ERROR("Invalid OSF standard deviation."); - ERROR(" (must be positive).\n"); - return 1; - } - break; - - case 4 : - full_stddev = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid full standard deviation.\n"); - return 1; - } - if ( full_stddev < 0.0 ) { - ERROR("Invalid full standard deviation."); - ERROR(" (must be positive).\n"); - return 1; - } - break; - - case 5 : - noise_stddev = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid noise standard deviation.\n"); - return 1; - } - if ( noise_stddev < 0.0 ) { - ERROR("Invalid noise standard deviation."); - ERROR(" (must be positive).\n"); - return 1; - } - break; - - case 6 : - image_prefix = strdup(optarg); - break; - - case 7 : - background = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid background level.\n"); - return 1; - } - if ( background < 0.0 ) { - ERROR("Background level must be positive.\n"); - return 1; - } - break; - - case 8 : - ERROR("--beam-divergence is no longer used.\n"); - ERROR("The 'xsphere' partiality model does not take divergence into account.\n"); - return 1; - - case 9 : - ERROR("--beam-bandwidth is no longer used.\n"); - ERROR("Set the bandwidth in the geometry file instead.\n"); - return 1; - - case 10 : - profile_radius = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid profile radius.\n"); - return 1; - } - if ( profile_radius < 0.0 ) { - ERROR("Profile radius must be positive.\n"); - return 1; - } - break; - - case 11 : - ERROR("--photon-energy is no longer used.\n"); - ERROR("Set the photon energy in the geometry file instead.\n"); - return 1; - - case 12 : - template = strdup(optarg); - break; - - case 0 : - break; - - case '?' : - return 1; - - default : - ERROR("Unhandled option '%c'\n", c); - break; - - } - } - - if ( n_threads < 1 ) { - ERROR("Invalid number of threads.\n"); - return 1; - } - - if ( (n_threads > 1) && (image_prefix != NULL) ) { - ERROR("Option \"--images\" is incompatible with \"-j\".\n"); - return 1; - } - - /* Load cell */ - if ( cellfile == NULL ) { - ERROR("You need to give a PDB file with the unit cell.\n"); - return 1; - } - cell = load_cell_from_file(cellfile); - if ( cell == NULL ) { - ERROR("Failed to get cell from '%s'\n", cellfile); - return 1; - } - free(cellfile); - - if ( !cell_is_sensible(cell) ) { - ERROR("Invalid unit cell parameters:\n"); - cell_print(cell); - return 1; - } - - /* Load geometry */ - if ( geomfile == NULL ) { - ERROR("You need to give a geometry file.\n"); - return 1; - } - dtempl = data_template_new_from_file(geomfile); - if ( dtempl == NULL ) { - ERROR("Failed to read geometry from '%s'\n", geomfile); - return 1; - } - - if ( save_file == NULL ) save_file = strdup("partial_sim.hkl"); - - /* Load (full) reflections */ - if ( input_file != NULL ) { - - RefList *as; - char *sym_str_fromfile = NULL; - - full = read_reflections_2(input_file, &sym_str_fromfile); - if ( full == NULL ) { - ERROR("Failed to read reflections from '%s'\n", - input_file); - return 1; - } - - /* If we don't have a point group yet, and if the file provides - * one, use the one from the file */ - if ( (sym_str == NULL) && (sym_str_fromfile != NULL) ) { - sym_str = sym_str_fromfile; - STATUS("Using symmetry from reflection file: %s\n", - sym_str); - } - - /* If we still don't have a point group, use "1" */ - if ( sym_str == NULL ) sym_str = strdup("1"); - - pointgroup_warning(sym_str); - sym = get_pointgroup(sym_str); - - if ( check_list_symmetry(full, sym) ) { - ERROR("The input reflection list does not appear to" - " have symmetry %s\n", symmetry_name(sym)); - if ( cell_get_lattice_type(cell) == L_MONOCLINIC ) { - ERROR("You may need to specify the unique axis " - "in your point group. The default is " - "unique axis c.\n"); - ERROR("See 'man crystfel' for more details.\n"); - } - return 1; - } - - as = asymmetric_indices(full, sym); - reflist_free(full); - full = as; - - } else { - random_intensities = 1; - if ( sym_str == NULL ) sym_str = strdup("1"); - sym = get_pointgroup(sym_str); - } - - if ( n < 1 ) { - ERROR("Number of patterns must be at least 1.\n"); - return 1; - } - - if ( output_file == NULL ) { - ERROR("You must give a filename for the output.\n"); - return 1; - } - stream = stream_open_for_write(output_file, dtempl); - if ( stream == NULL ) { - ERROR("Couldn't open output file '%s'\n", output_file); - return 1; - } - free(output_file); - - stream_write_geometry_file(stream, geomfile); - stream_write_commandline_args(stream, argc, argv); - - if ( template != NULL ) { - template_stream = stream_open_for_read(template); - if ( template_stream == NULL ) { - ERROR("Couldn't open template stream '%s'\n", template); - return 1; - } - } - - test_image = image_create_for_simulation(dtempl); - if ( test_image == NULL ) { - ERROR("Could not create image structure.\n"); - ERROR("Does the geometry file contain references?\n"); - return 1; - } - - STATUS("Simulation parameters:\n"); - STATUS(" Wavelength: %.5f A (photon energy %.2f eV)\n", - test_image->lambda*1e10, - ph_lambda_to_eV(test_image->lambda)); - STATUS(" Beam divergence: 0 (not modelled)\n"); - STATUS(" Beam bandwidth: %.5f %%\n", - test_image->bw*100.0); - STATUS("Reciprocal space profile radius: %e m^-1\n", - profile_radius); - if ( image_prefix != NULL ) { - STATUS(" Background: %.2f detector units\n", - background); - } else { - STATUS(" Background: none (no image " - "output)\n"); - } - STATUS(" Partiality model: xsphere (hardcoded)\n"); - STATUS(" Noise standard deviation: %.2f detector units\n", - noise_stddev); - if ( random_intensities ) { - STATUS(" Full intensities: randomly generated: " - "abs(Gaussian(sigma=%.2f)), symmetry %s\n", - full_stddev, sym_str); - } else { - STATUS(" Full intensities: from %s (symmetry %s)\n", - input_file, sym_str); - } - STATUS(" Max error in cell components: %.2f %%\n", cnoise); - STATUS("Scale factor standard deviation: %.2f\n", osf_stddev); - if ( template_stream != NULL ) { - STATUS("Crystal orientations and reflections to use from %s\n", - template); - } - - if ( random_intensities ) { - full = reflist_new(); - } - - qargs.full = full; - pthread_rwlock_init(&qargs.full_lock, NULL); - qargs.n_to_do = n; - qargs.n_done = 0; - qargs.n_started = 0; - qargs.sym = sym; - qargs.dtempl = dtempl; - qargs.random_intensities = random_intensities; - qargs.cell = cell; - qargs.stream = stream; - qargs.cnoise = cnoise; - qargs.osf_stddev = osf_stddev; - qargs.full_stddev = full_stddev; - qargs.noise_stddev = noise_stddev; - qargs.background = background; - qargs.max_q = detgeom_max_resolution(test_image->detgeom, - test_image->lambda); - qargs.image_prefix = image_prefix; - qargs.profile_radius = profile_radius; - qargs.template_stream = template_stream; - - qargs.rngs = malloc(n_threads * sizeof(gsl_rng *)); - if ( qargs.rngs == NULL ) { - ERROR("Failed to allocate RNGs\n"); - return 1; - } - - if ( config_random ) { - - FILE *fh; - - fh = fopen("/dev/urandom", "r"); - if ( fh == NULL ) { - ERROR("Failed to open /dev/urandom. Try again without" - " --really-random.\n"); - free(qargs.rngs); - return 1; - } - - for ( i=0; i<n_threads; i++ ) { - - unsigned long int seed; - qargs.rngs[i] = gsl_rng_alloc(gsl_rng_mt19937); - - if ( fread(&seed, sizeof(seed), 1, fh) == 1 ) { - gsl_rng_set(qargs.rngs[i], seed); - } else { - ERROR("Failed to seed RNG %i\n", i); - } - - } - - fclose(fh); - - } else { - gsl_rng *rng_for_seeds; - rng_for_seeds = gsl_rng_alloc(gsl_rng_mt19937); - for ( i=0; i<n_threads; i++ ) { - qargs.rngs[i] = gsl_rng_alloc(gsl_rng_mt19937); - gsl_rng_set(qargs.rngs[i], gsl_rng_get(rng_for_seeds)); - } - gsl_rng_free(rng_for_seeds); - } - - for ( i=0; i<NBINS; i++ ) { - qargs.n_ref[i] = 0; - qargs.p_hist[i] = 0.0; - qargs.p_max[i] = 0.0; - } - - run_threads(n_threads, run_job, create_job, finalise_job, - &qargs, n, 0, 0, 0); - - if ( random_intensities ) { - STATUS("Writing full intensities to %s\n", save_file); - write_reflist_2(save_file, full, sym); - } - - if ( phist_file != NULL ) { - - FILE *fh; - - fh = fopen(phist_file, "w"); - - if ( fh != NULL ) { - - double overall_max = 0.0; - double overall_mean = 0.0; - long long int overall_total = 0; - - for ( i=0; i<NBINS; i++ ) { - - double rcen; - - if ( qargs.p_max[i] > overall_max ) { - overall_max = qargs.p_max[i]; - } - - overall_mean += qargs.p_hist[i]; - overall_total += qargs.n_ref[i]; - - rcen = i/(double)NBINS*qargs.max_q - + qargs.max_q/(2.0*NBINS); - fprintf(fh, "%.2f %7lu %.3f %.3f\n", rcen/1.0e9, - qargs.n_ref[i], - qargs.p_hist[i]/qargs.n_ref[i], - qargs.p_max[i]); - - } - - fclose(fh); - - overall_mean /= overall_total; - - STATUS("Overall max partiality = %.2f\n", - overall_max); - STATUS("Overall mean partiality = %.2f\n", overall_mean); - STATUS("Total number of reflections = %lli\n", - overall_total); - - } else { - ERROR("Failed to open file '%s' for writing.\n", - phist_file); - } - - } - - for ( i=0; i<n_threads; i++ ) { - gsl_rng_free(qargs.rngs[i]); - } - free(qargs.rngs); - pthread_rwlock_destroy(&qargs.full_lock); - stream_close(stream); - cell_free(cell); - data_template_free(dtempl); - free_symoplist(sym); - reflist_free(full); - free(save_file); - free(geomfile); - free(input_file); - - return 0; -} diff --git a/src/pattern_sim.c b/src/pattern_sim.c deleted file mode 100644 index 9769a0d4..00000000 --- a/src/pattern_sim.c +++ /dev/null @@ -1,1183 +0,0 @@ -/* - * pattern_sim.c - * - * Simulate diffraction patterns from small crystals - * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2009-2020 Thomas White <taw@physics.org> - * 2013-2014 Chun Hong Yoon <chun.hong.yoon@desy.de> - * 2014 Valerio Mariani - * 2013 Alexandra Tolstikova - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <stdarg.h> -#include <stdlib.h> -#include <stdio.h> -#include <string.h> -#include <unistd.h> -#include <getopt.h> -#include <hdf5.h> - -#include <image.h> -#include <cell.h> -#include <cell-utils.h> -#include <utils.h> -#include <peaks.h> -#include <symmetry.h> -#include <reflist.h> -#include <reflist-utils.h> -#include <stream.h> -#include <datatemplate.h> - -#include "diffraction.h" -#include "diffraction-gpu.h" -#include "pattern_sim.h" -#include "version.h" - - -enum spectrum_type { - SPECTRUM_TOPHAT, /**< A top hat distribution */ - SPECTRUM_GAUSSIAN, /**< A Gaussian distribution */ - SPECTRUM_SASE, /**< A simulated SASE spectrum */ - SPECTRUM_TWOCOLOUR, /**< A spectrum containing two peaks */ - SPECTRUM_FROMFILE /**< An arbitrary spectrum read from a file */ -}; - -static void show_help(const char *s) -{ - printf("Syntax: %s [options]\n\n", s); - printf( -"Simulate diffraction patterns from small crystals probed with femtosecond\n" -"pulses of X-rays from a free electron laser.\n" -"\n" -" -h, --help Display this help message.\n" -" --version Print CrystFEL version number and exit.\n" -"\n" -" -p, --pdb=<file> File from which to get the unit cell.\n" -" (The actual Bragg intensities come from the\n" -" intensities file)\n" -" --gpu Use the GPU to speed up the calculation.\n" -" --gpu-dev=<n> Use GPU device <n>. Omit this option to see the\n" -" available devices.\n" -" -g, --geometry=<file> Get detector geometry from file.\n" -" -n, --number=<N> Generate N images. Default 1.\n" -" --no-images Do not output any HDF5 files.\n" -" -o, --output=<filename> Output HDF5 filename. Default: sim.h5.\n" -" -r, --random-orientation Use randomly generated orientations.\n" -" --powder=<file> Write a summed pattern of all images simulated by\n" -" 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> Select method for calculation of shape transforms\n" -" --really-random Seed the random number generator with /dev/urandom.\n" -" --min-size=<s> Minimum crystal size in nm.\n" -" --max-size=<s> Maximum crystal size in nm.\n" -" --no-noise Do not calculate Poisson noise.\n" -" -s, --sample-spectrum=<N> Use N samples from spectrum. Default 3.\n" -" -x, --spectrum=<type> Type of spectrum to simulate.\n" -" --background=<N> Add N photons of Poisson background (default 0).\n" -" --template=<file> Take orientations from stream <file>.\n" -" --no-fringes Exclude the side maxima of Bragg peaks.\n" -" --flat Make Bragg peaks flat.\n" -" --sase-spike-width SASE spike width (standard deviation in m^-1)\n" -" Default 2e6 m^-1\n" -" --twocol-separation Separation between peaks in two-colour mode in m^-1\n" -" Default 8e6 m^-1\n" -" --nphotons Number of photons per X-ray pulse. Default 1e12.\n" -" --beam-radius Radius of beam in metres (default 1e-6).\n" -); -} - - -static void record_panel(struct detgeom_panel *p, float *dp, - int do_poisson, - gsl_rng *rng, double ph_per_e, double background, - double lambda, - int *n_neg1, int *n_inf1, int *n_nan1, - int *n_neg2, int *n_inf2, int *n_nan2, - double *max_tt) -{ - int fs, ss; - - for ( ss=0; ss<p->h; ss++ ) { - for ( fs=0; fs<p->w; fs++ ) { - - double counts; - double intensity, sa; - double pix_area, Lsq; - double xs, ys, rx, ry; - double dsq, proj_area; - float dval; - double twotheta; - - intensity = (double)dp[fs + p->w*ss]; - if ( isinf(intensity) ) (*n_inf1)++; - if ( intensity < 0.0 ) (*n_neg1)++; - if ( isnan(intensity) ) (*n_nan1)++; - - /* Area of one pixel */ - pix_area = pow(p->pixel_pitch, 2.0); - Lsq = pow(p->cnz*p->pixel_pitch, 2.0); - - /* Calculate distance from crystal to pixel */ - xs = fs*p->fsx + ss*p->ssx; - ys = ss*p->fsy + ss*p->ssy; - rx = (xs + p->cnx) * p->pixel_pitch; - ry = (ys + p->cny) * p->pixel_pitch; - dsq = pow(rx, 2.0) + pow(ry, 2.0); - twotheta = atan2(sqrt(dsq), p->cnz*p->pixel_pitch); - - /* Area of pixel as seen from crystal (approximate) */ - proj_area = pix_area * cos(twotheta); - - /* Projected area of pixel divided by distance squared */ - sa = proj_area / (dsq + Lsq); - - if ( do_poisson ) { - counts = poisson_noise(rng, intensity * ph_per_e * sa); - } else { - counts = intensity * ph_per_e * sa; - } - - /* Number of photons in pixel */ - dval = counts + poisson_noise(rng, background); - - /* Convert to ADU */ - dval *= p->adu_per_photon; - - /* Saturation */ - if ( dval > p->max_adu ) dval = p->max_adu; - - dp[fs + p->w*ss] = dval; - - /* Sanity checks */ - if ( isinf(dp[fs + p->w*ss]) ) n_inf2++; - if ( isnan(dp[fs + p->w*ss]) ) n_nan2++; - if ( dp[fs + p->w*ss] < 0.0 ) n_neg2++; - - if ( twotheta > *max_tt ) *max_tt = twotheta; - - } - } -} - - -void record_image(struct image *image, int do_poisson, double background, - gsl_rng *rng, double beam_radius, double nphotons) -{ - double total_energy, energy_density; - double ph_per_e; - double area; - double max_tt = 0.0; - int n_inf1 = 0; - int n_neg1 = 0; - int n_nan1 = 0; - int n_inf2 = 0; - int n_neg2 = 0; - int n_nan2 = 0; - int pn; - - /* How many photons are scattered per electron? */ - area = M_PI*pow(beam_radius, 2.0); - total_energy = nphotons * ph_lambda_to_en(image->lambda); - energy_density = total_energy / area; - ph_per_e = (nphotons /area) * pow(THOMSON_LENGTH, 2.0); - STATUS("Fluence = %8.2e photons, " - "Energy density = %5.3f kJ/cm^2, " - "Total energy = %5.3f microJ\n", - nphotons, energy_density/1e7, total_energy*1e6); - - for ( pn=0; pn<image->detgeom->n_panels; pn++ ) { - - record_panel(&image->detgeom->panels[pn], image->dp[pn], - do_poisson, rng, ph_per_e, background, - image->lambda, - &n_neg1, &n_inf1, &n_nan1, - &n_neg2, &n_inf2, &n_nan2, &max_tt); - } - - STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n", - rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9); - - STATUS("Halve the d values to get the voxel size for a synthesis.\n"); - - if ( n_neg1 + n_inf1 + n_nan1 + n_neg2 + n_inf2 + n_nan2 ) { - ERROR("WARNING: The raw calculation produced %i negative" - " values, %i infinities and %i NaNs.\n", - n_neg1, n_inf1, n_nan1); - ERROR("WARNING: After processing, there were %i negative" - " values, %i infinities and %i NaNs.\n", - n_neg2, n_inf2, n_nan2); - } -} - - -static double *intensities_from_list(RefList *list, SymOpList *sym) -{ - Reflection *refl; - RefListIterator *iter; - double *out = new_arr_intensity(); - SymOpMask *m = new_symopmask(sym); - int neq = num_equivs(sym, NULL); - - for ( refl = first_refl(list, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) { - - signed int h, k, l; - int eps; - double intensity = get_intensity(refl); - - get_indices(refl, &h, &k, &l); - special_position(sym, m, h, k, l); - eps = neq / num_equivs(sym, m); - - set_arr_intensity(out, h, k, l, intensity / eps); - - } - - return out; -} - - -static double *phases_from_list(RefList *list) -{ - Reflection *refl; - RefListIterator *iter; - double *out = new_arr_phase(); - - for ( refl = first_refl(list, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) { - - signed int h, k, l; - double phase = get_phase(refl, NULL); - - get_indices(refl, &h, &k, &l); - - set_arr_phase(out, h, k, l, phase); - - } - - return out; - -} - - -static unsigned char *flags_from_list(RefList *list) -{ - Reflection *refl; - RefListIterator *iter; - unsigned char *out = new_arr_flag(); - - for ( refl = first_refl(list, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) { - - signed int h, k, l; - - get_indices(refl, &h, &k, &l); - - set_arr_flag(out, h, k, l, 1); - - } - - return out; - -} - - -static struct quaternion read_quaternion() -{ - do { - - int r; - float w, x, y, z; - char line[1024]; - char *rval; - - printf("Please input quaternion: w x y z\n"); - rval = fgets(line, 1023, stdin); - if ( rval == NULL ) return invalid_quaternion(); - chomp(line); - - r = sscanf(line, "%f %f %f %f", &w, &x, &y, &z); - if ( r == 4 ) { - - struct quaternion quat; - - quat.w = w; - quat.x = x; - quat.y = y; - quat.z = z; - - return quat; - - } else { - ERROR("Invalid rotation '%s'\n", line); - } - - } while ( 1 ); -} - - -static int random_ncells(double len, double min_m, double max_m) -{ - int min_cells, max_cells, wid; - - min_cells = min_m / len; - max_cells = max_m / len; - wid = max_cells - min_cells; - - return wid*random()/RAND_MAX + min_cells; -} - - -static void add_metadata(const char *filename, - struct quaternion q, - UnitCell *cell) -{ - hid_t fh, gh, sh, dh; - herr_t r; - hsize_t size[1]; - float data[4]; - double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; - - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - - fh = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT); - if ( fh < 0 ) { - ERROR("Failed to open file to add metadata.\n"); - return; - } - - gh = H5Gcreate2(fh, "pattern_sim", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - if ( gh < 0 ) { - ERROR("Failed to create metadata group.\n"); - H5Fclose(fh); - return; - } - - size[0] = 4; - sh = H5Screate_simple(1, size, NULL); - - dh = H5Dcreate2(gh, "quaternion", H5T_NATIVE_FLOAT, sh, - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - if ( dh < 0 ) { - ERROR("Couldn't create dataset\n"); - H5Fclose(fh); - return; - } - - data[0] = q.w; - data[1] = q.x; - data[2] = q.y; - data[3] = q.z; - r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, data); - if ( r < 0 ) { - ERROR("Failed to write quaternion to file\n"); - } - H5Dclose(dh); - - size[0] = 3; - sh = H5Screate_simple(1, size, NULL); - - dh = H5Dcreate2(gh, "astar", H5T_NATIVE_FLOAT, sh, - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - if ( dh < 0 ) { - ERROR("Couldn't create dataset\n"); - H5Fclose(fh); - return; - } - data[0] = asx; - data[1] = asy; - data[2] = asz; - r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, data); - if ( r < 0 ) { - ERROR("Failed to write astar to file.\n"); - } - H5Dclose(dh); - - dh = H5Dcreate2(gh, "bstar", H5T_NATIVE_FLOAT, sh, - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - if ( dh < 0 ) { - ERROR("Couldn't create dataset\n"); - H5Fclose(fh); - return; - } - data[0] = bsx; - data[1] = bsy; - data[2] = bsz; - r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, data); - if ( r < 0 ) { - ERROR("Failed to write bstar to file.\n"); - } - H5Dclose(dh); - - dh = H5Dcreate2(gh, "cstar", H5T_NATIVE_FLOAT, sh, - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - if ( dh < 0 ) { - ERROR("Couldn't create dataset\n"); - H5Fclose(fh); - return; - } - data[0] = csx; - data[1] = csy; - data[2] = csz; - r = H5Dwrite(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, data); - if ( r < 0 ) { - ERROR("Failed to write cstar to file.\n"); - } - H5Dclose(dh); - - H5Gclose(gh); - H5Fclose(fh); -} - - -int main(int argc, char *argv[]) -{ - int argn; - struct image *image; - DataTemplate *dtempl; - struct gpu_context *gctx = NULL; - struct image *powder; - char *intfile = NULL; - double *intensities; - char *rval; - double *phases; - unsigned char *flags; - int config_randomquat = 0; - int config_noimages = 0; - int config_nonoise = 0; - int config_nosfac = 0; - int config_gpu = 0; - int config_random = 0; - char *powder_fn = NULL; - char *cell_filename = NULL; - char *grad_str = NULL; - char *outfile = NULL; - char *geometry = NULL; - char *spectrum_str = NULL; - char *spectrum_fn = NULL; - GradientMethod grad; - enum spectrum_type spectrum_type; - int ndone = 0; /* Number of simulations done (images or not) */ - int number = 1; /* Number used for filename of image */ - int n_images = 1; /* Generate one image by default */ - int done = 0; - UnitCell *input_cell; - struct quaternion orientation; - int gpu_dev = -1; - int random_size = 0; - double min_size = 0.0; - double max_size = 0.0; - char *sym_str = NULL; - SymOpList *sym; - int nsamples = 3; - gsl_rng *rng; - double background = 0.0; - char *template_file = NULL; - Stream *st = NULL; - int no_fringes = 0; - int flat = 0; - double nphotons = 1e12; - double beam_radius = 1e-6; /* metres */ - double sase_spike_width = 2e6; /* m^-1 */ - double twocol_sep = 8e6; /* m^-1 */ - - /* Long options */ - const struct option longopts[] = { - {"help", 0, NULL, 'h'}, - {"version", 0, NULL, 'v'}, - {"gpu", 0, &config_gpu, 1}, - {"beam", 1, NULL, 'b'}, - {"random-orientation", 0, NULL, 'r'}, - {"number", 1, NULL, 'n'}, - {"no-images", 0, &config_noimages, 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'}, - {"output", 1, NULL, 'o'}, - {"geometry", 1, NULL, 'g'}, - {"sample-spectrum", 1, NULL, 's'}, - {"type-spectrum", 1, NULL, 'x'}, - {"spectrum", 1, NULL, 'x'}, - {"really-random", 0, &config_random, 1}, - {"no-fringes", 0, &no_fringes, 1}, - {"flat", 0, &flat, 1}, - - {"gpu-dev", 1, NULL, 2}, - {"min-size", 1, NULL, 3}, - {"max-size", 1, NULL, 4}, - {"background", 1, NULL, 5}, - {"template", 1, NULL, 6}, - {"beam-bandwidth", 1, NULL, 7}, - {"photon-energy", 1, NULL, 9}, - {"nphotons", 1, NULL, 10}, - {"beam-radius", 1, NULL, 11}, - {"spectrum-file", 1, NULL, 12}, - {"sase-spike-width", 1, NULL, 13}, - {"twocol-separation", 1, NULL, 14}, - - {0, 0, NULL, 0} - }; - - /* Short options */ - while ((argn = getopt_long(argc, argv, "hrn:i:t:p:o:g:y:s:x:vb:", - longopts, NULL)) != -1) { - - switch (argn) { - - case 'h' : - show_help(argv[0]); - return 0; - - case 'v' : - printf("CrystFEL: %s\n", - crystfel_version_string()); - printf("%s\n", - crystfel_licence_string()); - return 0; - - case 'b' : - ERROR("WARNING: This version of CrystFEL no longer " - "uses beam files. Please remove the beam file " - "from your pattern_sim command line.\n"); - return 1; - - case 'r' : - config_randomquat = 1; - break; - - case 'n' : - n_images = strtol(optarg, &rval, 10); - if ( *rval != '\0' ) { - ERROR("Invalid number of images.\n"); - return 1; - } - break; - - case 'i' : - intfile = strdup(optarg); - break; - - case 't' : - grad_str = strdup(optarg); - break; - - case 'p' : - cell_filename = strdup(optarg); - break; - - case 'o' : - outfile = strdup(optarg); - break; - - case 'w' : - powder_fn = strdup(optarg); - break; - - case 'g' : - geometry = strdup(optarg); - break; - - case 'y' : - sym_str = strdup(optarg); - break; - - case 's' : - nsamples = strtol(optarg, &rval, 10); - if ( *rval != '\0' ) { - ERROR("Invalid number of spectrum samples.\n"); - return 1; - } - break; - - case 'x' : - spectrum_str = strdup(optarg); - break; - - case 2 : - gpu_dev = atoi(optarg); - break; - - case 3 : - min_size = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid minimum size.\n"); - return 1; - } - min_size /= 1e9; - random_size++; - break; - - case 4 : - max_size = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid maximum size.\n"); - return 1; - } - max_size /= 1e9; - random_size++; - break; - - case 5 : - background = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid background level.\n"); - return 1; - } - break; - - case 6 : - template_file = strdup(optarg); - break; - - case 7 : - ERROR("--beam-bandwidth is no longer used.\n"); - ERROR("Set the bandwidth in the geometry file instead.\n"); - return 1; - - case 9 : - ERROR("--photon-energy is no longer used.\n"); - ERROR("Set the photon energy in the geometry file instead.\n"); - return 1; - - case 10 : - nphotons = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid number of photons.\n"); - return 1; - } - if ( nphotons < 0.0 ) { - ERROR("Number of photons must be positive.\n"); - return 1; - } - break; - - case 11 : - beam_radius = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid beam radius.\n"); - return 1; - } - if ( beam_radius < 0.0 ) { - ERROR("Beam radius must be positive.\n"); - return 1; - } - break; - - case 12 : - spectrum_fn = strdup(optarg); - break; - - case 13 : - sase_spike_width = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid SASE spike width.\n"); - return 1; - } - if ( beam_radius < 0.0 ) { - ERROR("SASE spike width must be positive.\n"); - return 1; - } - break; - - case 14 : - twocol_sep = strtod(optarg, &rval); - if ( *rval != '\0' ) { - ERROR("Invalid two colour separation.\n"); - return 1; - } - if ( beam_radius < 0.0 ) { - ERROR("Two-colour separation must be positive.\n"); - return 1; - } - break; - - case 0 : - break; - - case '?' : - break; - - default : - ERROR("Unhandled option '%c'\n", argn); - break; - - } - - } - - if ( random_size == 1 ) { - ERROR("You must specify both --min-size and --max-size.\n"); - return 1; - } - - if ( cell_filename == NULL ) { - cell_filename = strdup("molecule.pdb"); - } - - if ( outfile == NULL ) { - if ( n_images == 1 ) { - outfile = strdup("sim.h5"); - } else { - outfile = strdup("sim"); - } - } - - if ( template_file != NULL ) { - if ( config_randomquat ) { - ERROR("You cannot use -r and --template together.\n"); - return 1; - } - st = stream_open_for_read(template_file); - if ( st == NULL ) { - ERROR("Failed to open stream.\n"); - return 1; - } - free(template_file); - } - - if ( grad_str == NULL ) { - STATUS("You didn't specify a gradient calculation method, so" - " I'm using the 'mosaic' method, which is fastest.\n"); - grad = GRADIENT_MOSAIC; - } else if ( strcmp(grad_str, "mosaic") == 0 ) { - grad = GRADIENT_MOSAIC; - } else if ( strcmp(grad_str, "interpolate") == 0) { - grad = GRADIENT_INTERPOLATE; - } else if ( strcmp(grad_str, "phased") == 0) { - grad = GRADIENT_PHASED; - } else { - ERROR("Unrecognised gradient method '%s'\n", grad_str); - return 1; - } - free(grad_str); - - if ( config_gpu && (grad != GRADIENT_MOSAIC) ) { - ERROR("Only the mosaic method can be used for gradients when" - "calculating on the GPU.\n"); - return 1; - } - - if ( geometry == NULL ) { - ERROR("You need to specify a geometry file with --geometry\n"); - return 1; - } - dtempl = data_template_new_from_file(geometry); - if ( dtempl == NULL ) { - ERROR("Failed to read detector geometry from '%s'\n", - geometry); - return 1; - } - free(geometry); - - if ( spectrum_str == NULL ) { - STATUS("You didn't specify a spectrum type, so" - " I'm using a 'tophat' spectrum.\n"); - spectrum_type = SPECTRUM_TOPHAT; - } else if ( strcasecmp(spectrum_str, "tophat") == 0) { - spectrum_type = SPECTRUM_TOPHAT; - } else if ( strcasecmp(spectrum_str, "gaussian") == 0) { - spectrum_type = SPECTRUM_GAUSSIAN; - } else if ( strcasecmp(spectrum_str, "sase") == 0) { - spectrum_type = SPECTRUM_SASE; - } else if ( strcasecmp(spectrum_str, "twocolour") == 0 || - strcasecmp(spectrum_str, "twocolor") == 0 || - strcasecmp(spectrum_str, "twocolours") == 0 || - strcasecmp(spectrum_str, "twocolors") == 0) { - spectrum_type = SPECTRUM_TWOCOLOUR; - } else if ( strcasecmp(spectrum_str, "fromfile") == 0) { - spectrum_type = SPECTRUM_FROMFILE; - } else { - ERROR("Unrecognised spectrum type '%s'\n", spectrum_str); - return 1; - } - free(spectrum_str); - - /* Load unit cell */ - input_cell = load_cell_from_file(cell_filename); - if ( input_cell == NULL ) { - exit(1); - } - - 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; - - if ( sym_str == NULL ) sym_str = strdup("1"); - pointgroup_warning(sym_str); - sym = get_pointgroup(sym_str); - - } else { - - RefList *reflections; - char *sym_str_fromfile = NULL; - - reflections = read_reflections_2(intfile, &sym_str_fromfile); - if ( reflections == NULL ) { - ERROR("Problem reading input file %s\n", intfile); - return 1; - } - - /* If we don't have a point group yet, and if the file provides - * one, use the one from the file */ - if ( (sym_str == NULL) && (sym_str_fromfile != NULL) ) { - sym_str = sym_str_fromfile; - STATUS("Using symmetry from reflection file: %s\n", - sym_str); - } - - /* If we still don't have a point group, use "1" */ - if ( sym_str == NULL ) sym_str = strdup("1"); - - pointgroup_warning(sym_str); - sym = get_pointgroup(sym_str); - - if ( grad == GRADIENT_PHASED ) { - phases = phases_from_list(reflections); - } else { - phases = NULL; - } - intensities = intensities_from_list(reflections, sym); - flags = flags_from_list(reflections); - - /* Check that the intensities have the correct symmetry */ - if ( check_list_symmetry(reflections, sym) ) { - ERROR("The input reflection list does not appear to" - " have symmetry %s\n", symmetry_name(sym)); - if ( cell_get_lattice_type(input_cell) == L_MONOCLINIC ) - { - ERROR("You may need to specify the unique axis " - "in your point group. The default is " - "unique axis c.\n"); - ERROR("See 'man crystfel' for more details.\n"); - } - return 1; - } - - reflist_free(reflections); - - } - - rng = gsl_rng_alloc(gsl_rng_mt19937); - if ( config_random ) { - FILE *fh; - unsigned long int seed; - fh = fopen("/dev/urandom", "r"); - if ( fread(&seed, sizeof(seed), 1, fh) == 1 ) { - gsl_rng_set(rng, seed); - } else { - ERROR("Failed to seed random number generator\n"); - } - fclose(fh); - } - - image = image_create_for_simulation(dtempl); - powder = image_create_for_simulation(dtempl); - - /* Splurge a few useful numbers */ - STATUS("Simulation parameters:\n"); - STATUS(" Wavelength: %.5f A (photon energy %.2f eV)\n", - image->lambda*1e10, - ph_lambda_to_eV(image->lambda)); - STATUS(" Number of photons: %.0f (%.2f mJ)\n", - nphotons, eV_to_J(ph_lambda_to_eV(image->lambda))*nphotons*1e3); - STATUS(" Beam divergence: not simulated\n"); - STATUS(" Beam radius: %.2f microns\n", - beam_radius*1e6); - STATUS(" Background: %.2f photons\n", background); - - switch ( spectrum_type ) { - - case SPECTRUM_TOPHAT: - STATUS(" X-ray spectrum: top hat, " - "width %.5f %%\n", image->bw*100.0); - break; - - case SPECTRUM_GAUSSIAN: - STATUS(" X-ray spectrum: Gaussian, " - "bandwidth %.5f %%\n", image->bw*100.0); - break; - - case SPECTRUM_SASE: - STATUS(" X-ray spectrum: SASE, " - "bandwidth %.5f %%\n", image->bw*100.0); - break; - - case SPECTRUM_TWOCOLOUR: - STATUS(" X-ray spectrum: two colour, " - "bandwidth %.5f %%, separation %.5e m^-1\n", - image->bw*100.0, twocol_sep); - break; - - case SPECTRUM_FROMFILE: - STATUS(" X-ray spectrum: from %s\n", - spectrum_fn); - break; - } - if ( random_size ) { - STATUS(" Crystal size: random, between " - "%.2f and %.2f nm along each of a, b and c\n", - min_size*1e9, max_size*1e9); - } else { - STATUS(" Crystal size: 8 unit cells along " - "each of a, b and c\n"); - } - if ( intfile == NULL ) { - STATUS(" Full intensities: all equal\n"); - } else { - STATUS(" Full intensities: from %s (symmetry %s)\n", - intfile, sym_str); - } - do { - - int na, nb, nc; - double a, b, c, dis; - UnitCell *cell; - int err = 0; - int pi; - - /* Reset image data to zero for new pattern */ - for ( pi=0; pi<image->detgeom->n_panels; pi++ ) { - long j; - long np = image->detgeom->panels[pi].w - * image->detgeom->panels[pi].h; - for ( j=0; j<np; j++ ) { - image->dp[pi][j] = 0.0; - } - } - - if ( random_size ) { - - double alen, blen, clen, discard; - - cell_get_parameters(input_cell, &alen, &blen, &clen, - &discard, &discard, &discard); - - na = random_ncells(alen, min_size, max_size); - nb = random_ncells(blen, min_size, max_size); - nc = random_ncells(clen, min_size, max_size); - - } else { - - na = 8; - nb = 8; - nc = 8; - - } - - if ( st == NULL ) { - - if ( config_randomquat ) { - orientation = random_quaternion(rng); - } else { - orientation = read_quaternion(); - } - - STATUS("Orientation is %5.3f %5.3f %5.3f %5.3f\n", - orientation.w, orientation.x, - orientation.y, orientation.z); - - if ( !quaternion_valid(orientation) ) { - ERROR("Orientation modulus is not zero!\n"); - return 1; - } - - cell = cell_rotate(input_cell, orientation); - - } else { - - struct image *templ_image; - Crystal *cr; - - /* Get data from next chunk */ - templ_image = stream_read_chunk(st, 0); - if ( templ_image == NULL ) break; - if ( templ_image->n_crystals == 0 ) continue; - - cr = templ_image->crystals[0]; - cell = cell_new_from_cell(crystal_get_cell(cr)); - - if ( templ_image->n_crystals > 1 ) { - ERROR("Using the first crystal only.\n"); - } - - image_free(templ_image); - - } - - switch ( spectrum_type ) { - - case SPECTRUM_TOPHAT : - image->spectrum = spectrum_generate_tophat(image->lambda, - image->bw); - break; - - case SPECTRUM_GAUSSIAN : - image->spectrum = spectrum_generate_gaussian(image->lambda, - image->bw); - break; - - - case SPECTRUM_SASE : - image->spectrum = spectrum_generate_sase(image->lambda, - image->bw, - sase_spike_width, - rng); - break; - - case SPECTRUM_TWOCOLOUR : - image->spectrum = spectrum_generate_twocolour(image->lambda, - image->bw, - twocol_sep); - break; - - case SPECTRUM_FROMFILE : - image->spectrum = spectrum_load(spectrum_fn); - break; - - } - - cell_get_parameters(cell, &a, &b, &c, &dis, &dis, &dis); - STATUS("Particle size = %i x %i x %i" - " ( = %5.2f x %5.2f x %5.2f nm)\n", - na, nb, nc, na*a/1.0e-9, nb*b/1.0e-9, nc*c/1.0e-9); - - if ( config_gpu ) { - - if ( gctx == NULL ) { - gctx = setup_gpu(config_nosfac, - intensities, flags, sym_str, - gpu_dev); - } - err = get_diffraction_gpu(gctx, image, na, nb, nc, - cell, no_fringes, flat, - nsamples); - - } else { - get_diffraction(image, na, nb, nc, intensities, phases, - flags, cell, grad, sym, no_fringes, flat, - nsamples); - } - if ( err ) { - ERROR("Diffraction calculation failed.\n"); - done = 1; - goto skip; - } - - record_image(image, !config_nonoise, background, rng, - beam_radius, nphotons); - - if ( powder_fn != NULL ) { - - int pn; - - for ( pn=0; pn<image->detgeom->n_panels; pn++ ) { - - size_t w, i; - - w = image->detgeom->panels[pn].w - * image->detgeom->panels[pn].h; - - for ( i=0; i<w; i++ ) { - powder->dp[pn][i] += (double)image->dp[pn][i]; - } - - } - - if ( !(ndone % 10) ) { - image_write(powder, dtempl, powder_fn); - } - } - - if ( !config_noimages ) { - - char h5filename[1024]; - - if ( n_images != 1 ) { - snprintf(h5filename, 1023, "%s-%i.h5", - outfile, number); - } else { - strncpy(h5filename, outfile, 1023); - } - - number++; - - /* Write the output file */ - image_write(image, dtempl, h5filename); - - /* Add some pattern_sim-specific metadata */ - add_metadata(h5filename, orientation, cell); - - } - - cell_free(cell); - -skip: - ndone++; - - if ( n_images && (ndone >= n_images) ) done = 1; - - } while ( !done ); - - if ( powder_fn != NULL ) { - image_write(powder, dtempl, powder_fn); - } - - if ( gctx != NULL ) { - cleanup_gpu(gctx); - } - - image_free(image); - image_free(powder); - free(intfile); - cell_free(input_cell); - free(intensities); - free(outfile); - free(cell_filename); - free(sym_str); - free_symoplist(sym); - - gsl_rng_free(rng); - - return 0; -} diff --git a/src/pattern_sim.h b/src/pattern_sim.h deleted file mode 100644 index f013c0fb..00000000 --- a/src/pattern_sim.h +++ /dev/null @@ -1,67 +0,0 @@ -/* - * pattern_sim.h - * - * Simulate diffraction patterns from small crystals - * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * - * Authors: - * 2011-2021 Thomas White <taw@physics.org> - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#ifndef PATTERN_SIM_H -#define PATTERN_SIM_H - - -/* Maxmimum index to hold values up to (can be increased if necessary) - * WARNING: Altering this value constitutes an ABI change, and means you must - * update data/diffraction.cl then recompile and reinstall everything. */ -#define INDMAX 130 - -/* Array size */ -#define IDIM (INDMAX*2 +1) -#define LIST_SIZE (IDIM*IDIM*IDIM) - -/* Create functions for storing reflection intensities indexed as h,k,l */ -#define LABEL(x) x##_intensity -#define TYPE double -#include "list_tmp.h" - -/* As above, but for phase values */ -#define LABEL(x) x##_phase -#define TYPE double -#include "list_tmp.h" - -/* As above, but for (unsigned) integer counts */ -#define LABEL(x) x##_count -#define TYPE unsigned int -#include "list_tmp.h" - -/* As above, but for simple flags */ -#define LABEL(x) x##_flag -#define TYPE unsigned char -#include "list_tmp.h" - - -#endif /* PATTERN_SIM_H */ |