diff options
54 files changed, 2576 insertions, 1861 deletions
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index b234324d..cfe2054d 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -24,7 +24,7 @@ build_dev_images: - > /kaniko/executor --context $CI_PROJECT_DIR - --dockerfile $CI_PROJECT_DIR/Dockerfile_$IMAGE + --dockerfile $CI_PROJECT_DIR/deployment/docker/Dockerfile_$IMAGE --destination $CI_REGISTRY_IMAGE/$IMAGE:$IMAGE_TAG_LATEST --destination $CI_REGISTRY_IMAGE/$IMAGE:$IMAGE_TAG_SHA --build-arg HELM_VERSION=$HELM_VERSION @@ -91,9 +91,9 @@ build-and-test: build-and-test-minimal: image: $IMAGE:latest script: - - dnf update -y && dnf group install -y 'Development Tools' && - dnf install -y gsl-devel flex bison cmake meson - - meson setup build -Dhdf5=disabled --wrap-mode=nofallback && ninja -C build && ninja -C build test + - dnf update -y && dnf group install -y development-tools && + dnf install -y gsl-devel flex bison cmake meson + - meson setup build -Dhdf5=disabled --wrap-mode=nofallback && ninja -C build && ninja -C build test artifacts: when: always @@ -183,7 +183,7 @@ build_container_image: - > /kaniko/executor --context $CI_PROJECT_DIR - --dockerfile $CI_PROJECT_DIR/Dockerfile + --dockerfile $CI_PROJECT_DIR/deployment/docker/Dockerfile --destination $CI_REGISTRY_IMAGE/$IMAGE:$IMAGE_TAG_LATEST --destination $CI_REGISTRY_IMAGE/$IMAGE:$IMAGE_TAG_SHA --build-arg HELM_VERSION=$HELM_VERSION @@ -211,7 +211,7 @@ build_container_image_tag: - > /kaniko/executor --context $CI_PROJECT_DIR - --dockerfile $CI_PROJECT_DIR/Dockerfile + --dockerfile $CI_PROJECT_DIR/deployment/docker/Dockerfile --destination $CI_REGISTRY_IMAGE/$IMAGE:$IMAGE_TAG_VER --build-arg HELM_VERSION=$HELM_VERSION --build-arg RANCHER_CLI_VERSION=$RANCHER_CLI_VERSION diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 66d1ea62..7cda548a 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -14,7 +14,7 @@ The easiest way to ensure the future of CrystFEL is to use it and cite it in you You can also cite CrystFEL directly using one of the following DOIs: * All versions: [10.5281/zenodo.8183384](https://doi.org/10.5281/zenodo.8183384). -* This version (0.11.0): [10.5281/zenodo.10984700](https://doi.org/zenodo.10984700). +* This version (0.11.1): [10.5281/zenodo.13904047](https://doi.org/10.5281/zenodo.13904047). Getting started @@ -1,5 +1,5 @@ -Development version -------------------- +CrystFEL version 0.11.1, 8 October 2024 +--------------------------------------- - Introduction of Julia bindings for libcrystfel - General cleanup of libcrystfel API @@ -10,7 +10,7 @@ Development version - Catch NaNs in GGPM model [Gitlab #109] - GUI: Use SVG icons - indexamajig: Fix memory leak due to missing H5Sclose calls -- indexamajig: Fix buf with --peaks=hdf5 --no-revalidate +- indexamajig: Fix bug with --peaks=hdf5 --no-revalidate - indexamajig: Fix distorted cell distributions with xgandalf [GitLab #95] - indexamajig: Fix rare left-handed cells - indexamajig: Add ASAP::O acknowledgements @@ -25,10 +25,12 @@ Development version - GUI: Copy geometry file for detector refinement [Gitlab #101] - GUI: Fix unclosed bracket in script file [Gitlab #106] - Allow basic geometry refinement without hierarchy [Gitlab #96] +- Allow detector panel groups to appear at top of geometry file - Improve documentation about geometry hierarchy and refinement [GitHub #10] - Fix out-of-date documentation about geometry refinement [Gitlab #97] - Fix minor bug in unmerged sig(I) calculation at low photon counts - process_hkl: Fix ancient bug with --even-only multi-crystal indexing +- partialator: Add filename/frame ID to --unmerged-output CrystFEL version 0.11.0, 17 April 2024 @@ -246,7 +246,7 @@ Fedora 22 or later All dependencies can be taken from the Fedora repositories: ``` -$ sudo dnf group install 'Development Tools' +$ sudo dnf group install development-tools $ sudo dnf install hdf5-devel gsl-devel gtk3-devel cairo-devel pango-devel gdk-pixbuf2-devel meson gcc-c++ fftw-devel zeromq-devel msgpack-devel flex bison gcc-gfortran pandoc $ cd /home/user/downloads/crystfel $ meson build @@ -33,6 +33,7 @@ Documentation ------------- * [Basic Tutorial](doc/articles/tutorial.rst) +* [Which indexing method(s) should I use?](doc/articles/indexer-choice.rst) * [How to choose the right point group for merging](doc/articles/pointgroup.rst) * [How to increase data processing speed](doc/articles/speed.rst) * [Real-time data processing](doc/articles/online.rst) @@ -88,6 +89,26 @@ Citing CrystFEL Please see [CONTRIBUTING.md](CONTRIBUTING.md) for citation instructions. +Related software +---------------- + +[OnDA Monitor](https://www.ondamonitor.com/) for real-time monitoring of data +quality. Read [the paper](https://doi.org/10.1107/s1600576716007469). + +[Pixel Anomaly Detection Tool](https://github.com/gihankaushyal/PixelAnomalyDetectorTool) +uses machine learning techniques to find misbehaving detector pixels. +Read [the paper](https://doi.org/10.1107/s1600576724000116). + +[Careless](https://github.com/rs-station/careless) scales and merges data +using variational inference. Use its `stream2mtz` script to import data from +CrystFEL. Read [the paper](http://dx.doi.org/10.1038/s41467-022-35280-8). + +[DatView](https://github.com/nstander/DatView) helps with multivariate analysis +of large datasets. Read [the paper](https://doi.org/10.1107/s1600576719012044). +[More information here](https://sites.google.com/view/zatsepinlab/resources/datview) +including a tutorial and manual. + + Funding acknowledgements ------------------------ diff --git a/config.h.in b/config.h.in index c4ac7ebe..cf3599b9 100644 --- a/config.h.in +++ b/config.h.in @@ -7,3 +7,4 @@ #mesondefine HAVE_HDF5 #mesondefine HAVE_ASAPO #mesondefine HAVE_SCHED_SETAFFINITY +#mesondefine HAVE_FFTW diff --git a/crystfel.rb b/crystfel.rb index 3a132d39..59e233a6 100644 --- a/crystfel.rb +++ b/crystfel.rb @@ -6,8 +6,8 @@ class Crystfel < Formula homepage 'https://www.desy.de/~twhite/crystfel/index.html' stable do - url 'https://www.desy.de/~twhite/crystfel/crystfel-0.11.0.tar.gz' - sha256 'a203475987dbc51981e5374af0ad09a5c0a760e3b25daa501923e1f267e9f937' + url 'https://www.desy.de/~twhite/crystfel/crystfel-0.11.1.tar.gz' + sha256 'bd9b8df5d627a30c920bf397d04641d266e1a013b24623967d7f5deac97594a4' end head do diff --git a/Dockerfile b/deployment/docker/Dockerfile index 4dd9ab8a..4dd9ab8a 100644 --- a/Dockerfile +++ b/deployment/docker/Dockerfile diff --git a/Dockerfile_alma9 b/deployment/docker/Dockerfile_alma9 index 1d1e89c2..38df2b30 100644 --- a/Dockerfile_alma9 +++ b/deployment/docker/Dockerfile_alma9 @@ -59,7 +59,7 @@ RUN wget -nv https://github.com/msgpack/msgpack-c/releases/download/cpp-3.3.0 -DMSGPACK_BUILD_EXAMPLES=OFF \ -DMSGPACK_ENABLE_CXX=OFF \ -DCMAKE_BUILD_TYPE=Release ; \ - popd # build ; \ + popd ; \ make -C build -j 32 ; \ make -C build install ; \ popd # msgpack-3.3.0 diff --git a/Dockerfile_centos7 b/deployment/docker/Dockerfile_centos7 index 8c3cd696..5e591b04 100644 --- a/Dockerfile_centos7 +++ b/deployment/docker/Dockerfile_centos7 @@ -101,7 +101,7 @@ RUN wget -nv https://github.com/msgpack/msgpack-c/releases/download/cpp-3.3.0 -DMSGPACK_BUILD_EXAMPLES=OFF \ -DMSGPACK_ENABLE_CXX=OFF \ -DCMAKE_BUILD_TYPE=Release ; \ - popd # build ; \ + popd ; \ make -C build -j 32 ; \ make -C build install ; \ popd ; \ diff --git a/Dockerfile_fedora b/deployment/docker/Dockerfile_fedora index cfcb3c8c..52b9fb03 100644 --- a/Dockerfile_fedora +++ b/deployment/docker/Dockerfile_fedora @@ -8,7 +8,7 @@ ENV PATH=$CFPREFIX/bin:$PATH RUN dnf -y install epel-release ; \ dnf -y update ; \ - dnf -y group install 'Development Tools' ; \ + dnf -y group install development-tools; \ dnf -y install gsl-devel flex bison gcc-c++ gtk3-devel cairo-devel \ pango-devel gdk-pixbuf2-devel fftw-devel libpng-devel diffutils \ hdf5-devel zeromq-devel python3-devel which wget lz4 lz4-devel \ @@ -48,7 +48,7 @@ RUN wget -nv https://github.com/msgpack/msgpack-c/releases/download/cpp-3.3.0 -DMSGPACK_BUILD_EXAMPLES=OFF \ -DMSGPACK_ENABLE_CXX=OFF \ -DCMAKE_BUILD_TYPE=Release ; \ - popd # build ; \ + popd ; \ make -C build -j 32 ; \ make -C build install ; \ popd # msgpack-3.3.0 diff --git a/Dockerfile_ubuntu b/deployment/docker/Dockerfile_ubuntu index 161737fc..db2eafa3 100644 --- a/Dockerfile_ubuntu +++ b/deployment/docker/Dockerfile_ubuntu @@ -9,9 +9,9 @@ ENV PATH=$CFPREFIX/bin:/venv/bin:$PATH RUN export DEBIAN_FRONTEND=noninteractive; \ apt-get -y update ; \ - apt-get -q -y install build-essential flex bison cmake gfortran-14 libgsl-dev libgtk-3-dev \ + apt-get -q -y install build-essential flex bison cmake gfortran libgsl-dev libgtk-3-dev \ libcairo2-dev libpango1.0-dev libgdk-pixbuf2.0-dev libfftw3-dev libpng-dev libmsgpack-dev \ - libeigen3-dev libccp4-dev libhdf5-dev libzmq3-dev libpython3-dev python3-full python3-venv \ + libeigen3-dev libccp4-dev libhdf5-dev libzmq3-dev libpython3-dev python3-full python3-venv \ liblz4-dev libbz2-dev libcurl4-openssl-dev debianutils wget lz4 diffutils openssh-client \ krb5-user sssd-krb5 git curl unzip @@ -49,7 +49,7 @@ RUN wget -nv https://github.com/msgpack/msgpack-c/releases/download/cpp-3.3.0 -DMSGPACK_BUILD_EXAMPLES=OFF \ -DMSGPACK_ENABLE_CXX=OFF \ -DCMAKE_BUILD_TYPE=Release ; \ - popd # build ; \ + popd ; \ make -C build -j 32 ; \ make -C build install ; \ popd # msgpack-3.3.0 diff --git a/deployment/gitlab/config.toml b/deployment/gitlab/config.toml new file mode 100644 index 00000000..92e1ad3e --- /dev/null +++ b/deployment/gitlab/config.toml @@ -0,0 +1,16 @@ +[[runners]] + name = "macos_runner.desy.de" + url = "https://gitlab.desy.de/" + id = 9999 + token = "registration token" + token_obtained_at = 2023-05-10T11:19:56Z + token_expires_at = 0001-01-01T00:00:00Z + executor = "shell" + [runners.custom_build_dir] + [runners.cache] + MaxUploadedArchiveSize = 0 + [runners.cache.s3] + [runners.cache.gcs] + [runners.cache.azure] + [runners.feature_flags] + FF_ENABLE_JOB_CLEANUP = true diff --git a/doc/articles/indexer-choice.rst b/doc/articles/indexer-choice.rst new file mode 100644 index 00000000..57b1389f --- /dev/null +++ b/doc/articles/indexer-choice.rst @@ -0,0 +1,97 @@ +====================================== +Which indexing method(s) should I use? +====================================== + +The short answer +================ + +For most applications, the best choice of indexing method is ``xgandalf``. + +In most cases, it's not necessary to use several algorithms in sequence +(e.g. ``--indexing=xgandalf,mosflm,asdf``). + +Read on for recommendations in specific cases. + + +Unknown crystal lattice parameters +================================== + +Most of the indexing algorithms run fine without prior information about the +lattice parameters. The indexing algorithms will generate a set of estimated +lattice parameters for each crystal, and you can use ``cell_explorer`` to plot +histograms to determine the most "popular" parameters. + +Problems arise because there are an infinite number of ways to represent any +given lattice. Different indexing algorithms have different "opinions" +about which representation to use when not constrained by prior information. + +In this case, ``mosflm`` is a good choice of indexing method because it can +detect centered lattices instead of producing a primitive unit cell every time. +This makes it easier and quicker to see the correct symmetry (including indexing +ambiguities) and get to the crystallographically conventional representation of +the structure. + +Other indexing methods, particularly ``xgandalf``, are fine as long as you use +proper crystallographic knowledge to determine the symmetry for merging and +structure solution. + +It's particularly advisable to avoid using multiple indexing methods together +when the lattice parameters are unknown, because the "preferences" of different +algorithms will be mixed up together. + + +Fastest indexing (e.g. real-time analysis pipelines) +==================================================== + +Read the `document about data processing speed <speed.rst>`_. The best choice +is probably ``asdf`` or ``xgandalf`` with ``--xgandalf-fast``, but there are +many other considerations for increasing processing speed. + + +Wide-bandwidth X-ray diffraction (Laue/pink beam) +================================================= + +Use ``pinkindexer``, especially if the (FWHM) bandwidth is larger than about 3% +of the wavelength. For smaller bandwidths, try with ``xgandalf``. + +For larger bandwidths, it's also advisable to use options ``--no-refine`` and +``--no-check-peaks``, because these parts of the program assume that the +bandwidth is small. + + +Electron diffraction data +========================= + +Use ``pinkindexer``, but first read the `document about electrons <electrons.rst>`_. + + +Small unit cell +=============== + +In this version of CrystFEL, none of the algorithms are very well-suited for +small unit cells. Nevertheless, some researchers have had success using +``xgandalf``. It's also worth trying with ``taketwo`` because of the way its +algorithms are devised. + +A specific algorithm for small unit cells is expected in CrystFEL 0.12.0. + + +Multiple overlapping diffraction patterns +========================================= + +Try ``felix``, which is specifically designed for heavily overlapping +diffraction patterns. Note, however, that it's unfortunately difficult to get +hold of a copy of the Felix executable. + + +CrystFEL was installed minimally +================================ + +The ``taketwo`` algorithm has no dependencies, neither compile-time nor +run-time, so is always available. Note that ``taketwo`` always needs prior +lattice parameters. + +Even if you have no system-wide access to the computer system, you can still +make a user-local installation of ``dirax``, ``mosflm``, ``xds`` or ``felix``. +These methods have no compile-time requirements, and only need the +corresponding programs to be available at run-time. diff --git a/doc/articles/julia.rst b/doc/articles/julia.rst index 0b6cabeb..206848f8 100644 --- a/doc/articles/julia.rst +++ b/doc/articles/julia.rst @@ -71,6 +71,16 @@ or:: UnitCell(RhombohedralLattice, RhombohedralCell, NoUniqueAxis, 23.0 Å, 23.0 Å, 23.0 Å, 75.0°, 75.0°, 75.0°) +A ``UnitCell`` also keeps a record of the orientation of the lattice. You can +get the real and reciprocal space matrices of lattice vectors with +``cell.directcartesian`` and ``cell.reciprocalcartesian``, respectively:: + + julia> uc.directcartesian + 3×3 Matrix{Float64}: + 1.37e-8 8.38883e-25 8.38883e-25 + 0.0 1.37e-8 8.38883e-25 + 0.0 0.0 1.37e-8 + Reflection lists ================ @@ -212,12 +222,47 @@ same ``Image``). If necessary, you can make a copy using ``deepcopy``. Crystals ======== +A ``Crystal`` is made up of a ``UnitCell`` (which includes the orientation of +the crystal, as mentioned above) and a few other parameters:: + + julia> cr = Crystal(uc) + CrystFEL.Crystal(0x000000001ad0ada0): + + UnitCell(CubicLattice, FaceCenteredCell, NoUniqueAxis, + 137.000 Å, 137.000 Å, 137.000 Å, 90.000°, 90.000°, 90.000°) + + Linear scale factor: 1.0 + Debye-Walle factor: 0.0 + Mosaicity: 0.0 + Profile radius: 0.002 nm⁻¹ + Resolution limit: Inf + Flag: 0 Indexing ======== +Create an indexing engine like this:: + + indexer = Indexer("asdf", dtempl, cell) + +The first argument is handled exactly as by ``indexamajig --indexing=``. The +second argument is the detector geometry (``DataTemplate``), and you also need +to provide the target unit cell. + +Run the indexing engine on an image like this:: + + index(image, indexer) + Prediction ========== +Given a ``Crystal`` and and ``Image``, the function +``predictreflections(cr, image)`` will return a ``RefList{UnmergedReflection}`` +containing the predicted reflections. + +You can subsequently calculate partialities with:: + + calculatepartialities!(reflist, cr, image, model=XSphereModel) + diff --git a/doc/articles/tutorial.rst b/doc/articles/tutorial.rst index 9e02b2f3..4301b64b 100644 --- a/doc/articles/tutorial.rst +++ b/doc/articles/tutorial.rst @@ -775,7 +775,7 @@ refined geometry file. Choose any suitable new filename, then look at the bottom of the window. There is a drop-down menu where you should select the name of the indexing job from earlier in this step (using prior unit cell parameters). The geometry refinement calculations will be based on the results -of that indexing job. Leave **Heirarchy level** at zero and **Include +of that indexing job. Leave **Hierarchy level** at zero and **Include out-of-plane positions and tilts** unselected, to keep the refinement simple. This will only refine the overall beam center position on the detector. Higher refinement levels will do a more fine-grained refinement where the subunits diff --git a/doc/man/crystfel_geometry.5.md b/doc/man/crystfel_geometry.5.md index 640771ee..0c5dab1e 100644 --- a/doc/man/crystfel_geometry.5.md +++ b/doc/man/crystfel_geometry.5.md @@ -135,7 +135,7 @@ For each panel, the physical location is controlled by the **fs**, **ss**, **coffset** : The offset of the panel along the z-direction from the position given by -: **clen**. +: **clen**, in metres. The overall detector position in the z-direction is given by **clen**, which can only be specified once in the geometry file (not for each panel): @@ -535,9 +535,15 @@ above. The highest-level group should always be called **all**. All members of a group need to be defined before defining the group. This -means that the group definitions must come **after** the panel definitions, and -the groups should be defined from the bottom to top level of the hierarchy - -the **all** group coming last. +means that the group definitions shuld come **after** the panel definitions, +and the groups should be defined from the bottom to top level of the hierarchy +with the **all** group coming last. As a convenience, CrystFEL 0.11.1 and +later defer processing **group** lines until the rest of the geometry file +has been read. That means that you can put your **group** lines near the top +of the file, but the groups must still be defined from bottom to top level. +To maintain compatability with CrystFEL 0.11.0, tools that write geometry +files (e.g. **align_detector**) from later CrystFEL versions will put the group +directives at the end of the file. If you do not define any groups, CrystFEL will automatically create the **all** group for you, containing all panels in a flat hierarchy. This allows basic diff --git a/doc/man/indexamajig.1.md b/doc/man/indexamajig.1.md index 605782bd..d93df7ca 100644 --- a/doc/man/indexamajig.1.md +++ b/doc/man/indexamajig.1.md @@ -508,7 +508,8 @@ PEAK SEARCH OPTIONS **--max-res=px** : Only accept peaks if they lay at less than px pixels from the center of the -: detector when using **--peaks=peakfinder8**. The default is **--max-res=1200**. +: detector when using **--peaks=peakfinder8**. The default is **--max-res=0**, +: which disables the cutoff. **--no-use-saturated** : Normally, peaks which contain one or more pixels above max_adu (defined in the @@ -595,8 +596,9 @@ INDEXING OPTIONS **--mille** : Write detector calibration data in Millepede-II format. -**--mille-dir=dirname** -: Write the Millepede-II data into _dirname_. +**--mille-file=filename** +: Write the Millepede-II data into _filename_. The default is mille-data.bin, +: in the current directory. **--max-mille-level=n** : Write the Millepede-II data up to a maximum hierarchy depth _n_. If _n_ is diff --git a/libcrystfel/src/crystfel-mille.c b/libcrystfel/src/crystfel-mille.c index c2fe2a57..a62bf7b2 100644 --- a/libcrystfel/src/crystfel-mille.c +++ b/libcrystfel/src/crystfel-mille.c @@ -297,7 +297,7 @@ void write_mille(Mille *mille, int n, UnitCell *cell, } -Mille *crystfel_mille_new(const char *outFileName) +static Mille *mille_new() { Mille *m; @@ -311,6 +311,15 @@ Mille *crystfel_mille_new(const char *outFileName) m->have_local = NULL; m->n_local = 0; + return m; +} + + +Mille *crystfel_mille_new(const char *outFileName) +{ + Mille *m = mille_new(); + if ( m == NULL ) return NULL; + m->fh = fopen(outFileName, "wb"); if ( m->fh == NULL ) { ERROR("Failed to open Mille file '%s'\n", outFileName); @@ -318,6 +327,21 @@ Mille *crystfel_mille_new(const char *outFileName) return NULL; } + return m; +} + + +Mille *crystfel_mille_new_fd(int fd) +{ + Mille *m = mille_new(); + if ( m == NULL ) return NULL; + + m->fh = fdopen(fd, "wb"); + if ( m->fh == NULL ) { + ERROR("Failed to open Mille FD %i\n", fd); + cffree(m); + return NULL; + } return m; } diff --git a/libcrystfel/src/crystfel-mille.h b/libcrystfel/src/crystfel-mille.h index d6602ef0..54588e16 100644 --- a/libcrystfel/src/crystfel-mille.h +++ b/libcrystfel/src/crystfel-mille.h @@ -43,6 +43,7 @@ typedef struct mille Mille; */ extern Mille *crystfel_mille_new(const char *outFileName); +extern Mille *crystfel_mille_new_fd(int fd); extern void crystfel_mille_free(Mille *m); diff --git a/libcrystfel/src/datatemplate.c b/libcrystfel/src/datatemplate.c index c97bd6d7..a921fe37 100644 --- a/libcrystfel/src/datatemplate.c +++ b/libcrystfel/src/datatemplate.c @@ -1890,7 +1890,7 @@ static int all_panels_perpendicular_to_beam(const DataTemplate *dtempl) double z_diff; struct panel_template *p = &dtempl->panels[i]; z_diff = p->fsz*PANEL_WIDTH(p) + p->ssz*PANEL_HEIGHT(p); - if ( z_diff > 10.0*p->pixel_pitch ) return 0; + if ( z_diff > 10.0 ) return 0; } return 1; } diff --git a/libcrystfel/src/datatemplate_priv.h b/libcrystfel/src/datatemplate_priv.h index 26ee91ed..5fe104fe 100644 --- a/libcrystfel/src/datatemplate_priv.h +++ b/libcrystfel/src/datatemplate_priv.h @@ -41,7 +41,7 @@ #define MAX_PATH_PARTS (16) /* Maximum number of panel groups */ -#define MAX_PANEL_GROUPS (256) +#define MAX_PANEL_GROUPS (512) /* Maximum number of panel groups that can derive from one panel */ #define MAX_PANEL_GROUP_CHILDREN (64) @@ -128,7 +128,7 @@ struct panel_template double cny; /**@}*/ - /** The offset to be applied from clen */ + /** The offset to be applied from clen (in m) */ double cnz_offset; /** Mask definitions */ diff --git a/libcrystfel/src/image-seedee.c b/libcrystfel/src/image-seedee.c index 1b68aead..544efb00 100644 --- a/libcrystfel/src/image-seedee.c +++ b/libcrystfel/src/image-seedee.c @@ -98,6 +98,19 @@ static int load_seedee_data(struct panel_template *p, } } + } else if ( (array->datatype == 'u') + && (array->itemsize == 1)) + { + int fs, ss; + uint8_t *in_data = (uint8_t *)array->data; + + for ( ss=0; ss<PANEL_HEIGHT(p); ss++ ) { + for ( fs=0; fs<PANEL_WIDTH(p); fs++ ) { + size_t idx = fs+p->orig_min_fs + (ss+p->orig_min_ss)*data_size_fs; + data[fs+ss*PANEL_WIDTH(p)] = in_data[idx]; + } + } + } else if ( (array->datatype == 'i') && (array->itemsize == 4) && (array->byteorder == '<') ) diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c index 3b215c08..6b7b1696 100644 --- a/libcrystfel/src/image.c +++ b/libcrystfel/src/image.c @@ -1465,7 +1465,6 @@ struct image *image_new() image->ida = NULL; image->n_cached_headers = 0; - image->id = 0; image->serial = 0; image->spectrum = NULL; image->lambda = -1.0; diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 3db7de18..7a21aff2 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -165,9 +165,6 @@ struct image struct header_cache_entry *header_cache[HEADER_CACHE_SIZE]; int n_cached_headers; - /** ID number of the worker processing handling this image */ - int id; - /** Monotonically increasing serial number for this image */ int serial; diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 5b0f742f..66b6e5de 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -87,36 +87,6 @@ static const char *onoff(int a) return "off"; } -/* Definition and function definition duplicated here (from im-sandbox.{c,h}) - * because libcrystfel code cannot depend on core CrystFEL programs. - * - * Must match the value and definition in im-sandbox.h - */ -#define MAX_TASK_LEN (32) - -static void set_last_task(char *lt, const char *task) -{ - if ( lt == NULL ) return; - assert(strlen(task) < MAX_TASK_LEN-1); - strcpy(lt, task); -} - - -static void show_indexing_flags(IndexingFlags flags) -{ - STATUS("Indexing parameters:\n"); - STATUS(" Check unit cell parameters: %s\n", - onoff(flags & INDEXING_CHECK_CELL)); - STATUS(" Check peak alignment: %s\n", - onoff(flags & INDEXING_CHECK_PEAKS)); - STATUS(" Refine indexing solutions: %s\n", - onoff(flags & INDEXING_REFINE)); - STATUS(" Multi-lattice indexing (\"delete and retry\"): %s\n", - onoff(flags & INDEXING_MULTI)); - STATUS(" Retry indexing: %s\n", - onoff(flags & INDEXING_RETRY)); -} - char *base_indexer_str(IndexingMethod indm) { @@ -445,17 +415,34 @@ IndexingPrivate *setup_indexing(const char *method_list, } for ( i=0; i<6; i++ ) ipriv->tolerance[i] = tols[i]; + return ipriv; +} + + +void print_indexing_info(IndexingPrivate *ipriv) +{ + int i; + STATUS("List of indexing methods:\n"); - for ( i=0; i<n; i++ ) { - char *str = indexer_str(methods[i]); - char *tmp = friendly_indexer_name(methods[i]); + for ( i=0; i<ipriv->n_methods; i++ ) { + char *str = indexer_str(ipriv->methods[i]); + char *tmp = friendly_indexer_name(ipriv->methods[i]); STATUS(" %2i: %-25s (%s)\n", i, str, tmp); cffree(str); cffree(tmp); } - show_indexing_flags(flags); - return ipriv; + STATUS("Indexing parameters:\n"); + STATUS(" Check unit cell parameters: %s\n", + onoff(ipriv->flags & INDEXING_CHECK_CELL)); + STATUS(" Check peak alignment: %s\n", + onoff(ipriv->flags & INDEXING_CHECK_PEAKS)); + STATUS(" Refine indexing solutions: %s\n", + onoff(ipriv->flags & INDEXING_REFINE)); + STATUS(" Multi-lattice indexing (\"delete and retry\"): %s\n", + onoff(ipriv->flags & INDEXING_MULTI)); + STATUS(" Retry indexing: %s\n", + onoff(ipriv->flags & INDEXING_RETRY)); } @@ -579,94 +566,78 @@ static int check_cell(IndexingFlags flags, Crystal *cr, UnitCell *target, } -#ifdef MEASURE_INDEX_TIME -static float real_time() -{ - struct timespec tp; - clock_gettime(CLOCK_MONOTONIC_RAW, &tp); - return tp.tv_sec + tp.tv_nsec*1e-9; -} -#endif - - /* Return non-zero for "success" */ static int try_indexer(struct image *image, IndexingMethod indm, - IndexingPrivate *ipriv, void *mpriv, char *last_task, + IndexingPrivate *ipriv, void *mpriv, Mille *mille, int max_mille_level) { int i, r; int n_bad = 0; int n_before = image->n_crystals; - #ifdef MEASURE_INDEX_TIME - float time_start; - float time_end; - time_start = real_time(); - #endif - switch ( indm & INDEXING_METHOD_MASK ) { case INDEXING_NONE : return 0; case INDEXING_DIRAX : - set_last_task(last_task, "indexing:dirax"); + set_last_task("indexing:dirax"); profile_start("dirax"); r = run_dirax(image, mpriv); profile_end("dirax"); break; case INDEXING_ASDF : - set_last_task(last_task, "indexing:asdf"); + set_last_task("indexing:asdf"); profile_start("asdf"); r = run_asdf(image, mpriv); profile_end("asdf"); break; case INDEXING_MOSFLM : - set_last_task(last_task, "indexing:mosflm"); + set_last_task("indexing:mosflm"); profile_start("mosflm"); r = run_mosflm(image, mpriv); profile_end("mosflm"); break; case INDEXING_XDS : - set_last_task(last_task, "indexing:xds"); + set_last_task("indexing:xds"); profile_start("xds"); r = run_xds(image, mpriv); profile_end("xds"); break; case INDEXING_FILE : - set_last_task(last_task, "indexing:file"); + set_last_task("indexing:file"); profile_start("fromfile"); r = fromfile_index(image, mpriv); profile_end("fromfile"); break; case INDEXING_FELIX : - set_last_task(last_task, "indexing:felix"); + set_last_task("indexing:felix"); profile_start("felix"); r = felix_index(image, mpriv); profile_end("felix"); break; case INDEXING_TAKETWO : - set_last_task(last_task, "indexing:taketwo"); + set_last_task("indexing:taketwo"); profile_start("taketwo"); r = taketwo_index(image, mpriv); profile_end("taketwo"); break; case INDEXING_PINKINDEXER : - set_last_task(last_task, "indexing:pinkindexer"); + set_last_task("indexing:pinkindexer"); profile_start("pinkindexer"); r = run_pinkIndexer(image, mpriv, ipriv->n_threads); profile_end("pinkindexer"); break; case INDEXING_XGANDALF : - set_last_task(last_task, "indexing:xgandalf"); + set_last_task("indexing:xgandalf"); profile_start("xgandalf"); r = run_xgandalf(image, mpriv); profile_end("xgandalf"); @@ -678,11 +649,7 @@ static int try_indexer(struct image *image, IndexingMethod indm, } - set_last_task(last_task, "indexing:finalisation"); - - #ifdef MEASURE_INDEX_TIME - time_end = real_time(); - #endif + set_last_task("indexing:finalisation"); /* Stop a really difficult to debug situation in its tracks */ if ( image->n_crystals - n_before != r ) { @@ -786,16 +753,6 @@ static int try_indexer(struct image *image, IndexingMethod indm, n_bad = remove_flagged_crystals(image); assert(r >= n_bad); - #ifdef MEASURE_INDEX_TIME - printf("%s took %f s, %i crystals found of which %i accepted. %s %s\n", - indexer_str(indm & INDEXING_METHOD_MASK), - time_end - time_start, - r, r - n_bad, - image->filename, - image->ev); - fflush(stdout); - #endif - return r - n_bad; } @@ -998,13 +955,13 @@ void index_pattern_4(struct image *image, IndexingPrivate *ipriv, int *ping, r = try_indexer(image, ipriv->methods[n], ipriv, ipriv->engine_private[n], - last_task, mille, max_mille_level); + mille, max_mille_level); success += r; ntry++; done = finished_retry(ipriv->methods[n], ipriv->flags, r, image); if ( ntry > 5 ) done = 1; - if ( ping != NULL ) (*ping)++; + notify_alive(); } while ( !done ); @@ -1031,6 +988,13 @@ void index_pattern_4(struct image *image, IndexingPrivate *ipriv, int *ping, } +void index_pattern_5(struct image *image, IndexingPrivate *ipriv, + Mille *mille, int max_mille_level) +{ + index_pattern_4(image, ipriv, NULL, NULL, mille, max_mille_level); +} + + /* Set the indexer flags for "use no lattice type information" */ static IndexingMethod set_nolattice(IndexingMethod a) { diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index 10d07f08..c7e9efe3 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -240,6 +240,8 @@ extern IndexingPrivate *setup_indexing(const char *methods, struct fromfile_options *fromfile_opts, struct asdf_options *asdf_opts); +extern void print_indexing_info(IndexingPrivate *ipriv); + extern const IndexingMethod *indexing_methods(IndexingPrivate *p, int *n); extern char *detect_indexing_methods(UnitCell *cell); @@ -256,6 +258,9 @@ extern void index_pattern_4(struct image *image, IndexingPrivate *ipriv, int *ping, char *last_task, Mille *mille, int max_mille_level); +extern void index_pattern_5(struct image *image, IndexingPrivate *ipriv, + Mille *mille, int max_mille_level); + extern void cleanup_indexing(IndexingPrivate *ipriv); #ifdef __cplusplus diff --git a/libcrystfel/src/peakfinder8.c b/libcrystfel/src/peakfinder8.c index 28bee8d7..5822bba9 100644 --- a/libcrystfel/src/peakfinder8.c +++ b/libcrystfel/src/peakfinder8.c @@ -399,8 +399,9 @@ static struct peakfinder_mask *create_peakfinder_mask(const struct image *img, idx = ifs + iss*p.w; - if ( rmps->r_maps[i][idx] < max_res - && rmps->r_maps[i][idx] > min_res ) { + if ( ((max_res == 0) || (rmps->r_maps[i][idx] < max_res)) + && rmps->r_maps[i][idx] > min_res ) + { if (! ( ( img->bad != NULL ) && ( img->bad[i] != NULL ) diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c index f1f2173c..b64859bf 100644 --- a/libcrystfel/src/stream.c +++ b/libcrystfel/src/stream.c @@ -1174,9 +1174,9 @@ Stream *stream_open_for_write(const char *filename, } -int stream_get_fd(Stream *st) +FILE *stream_get_fh(Stream *st) { - return fileno(st->fh); + return st->fh; } diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h index beb620c7..b32f4927 100644 --- a/libcrystfel/src/stream.h +++ b/libcrystfel/src/stream.h @@ -36,6 +36,8 @@ * Stream functions (for indexing results) */ +#include <stdio.h> + struct image; #include "datatemplate.h" @@ -112,7 +114,7 @@ extern char *stream_audit_info(Stream *st); extern char *stream_geometry_file(Stream *st); /* Low-level stuff used for indexamajig sandbox */ -extern int stream_get_fd(Stream *st); +extern FILE *stream_get_fh(Stream *st); extern int stream_rewind(Stream *st); /* Random access */ diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c index a35ba272..96e2a371 100644 --- a/libcrystfel/src/utils.c +++ b/libcrystfel/src/utils.c @@ -466,7 +466,43 @@ char *safe_strdup(const char *in) } -/* ------------------------------ Useful functions ---------------------------- */ +/* -------------------------------- Debugging ------------------------------- */ + +static void set_last_task_dummy(const char *task) +{ + /* Do nothing */ +} + +static void notify_alive_dummy() +{ + /* Do nothing */ +} + +struct _debugconf { + void (*set_last_task)(const char *task); + void (*notify_alive)(void); +} debug_conf = { set_last_task_dummy, notify_alive_dummy }; + +int set_debug_funcs(void (*slt)(const char *), + void (*ping)(void)) +{ + debug_conf.set_last_task = slt; + debug_conf.notify_alive = ping; + return 0; +} + +void set_last_task(const char *task) +{ + debug_conf.set_last_task(task); +} + +void notify_alive() +{ + debug_conf.notify_alive(); +} + + +/* ----------------------------- Useful functions --------------------------- */ int convert_int(const char *str, int *pval) { diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index 08e059a7..e295c40b 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -254,6 +254,14 @@ extern int set_mm_funcs(void *(*cfmalloc)(size_t size), void *(*cfrealloc)(void *ptr, size_t size)); +/* -------------------------------- Debugging ------------------------------- */ + +extern void set_last_task(const char *task); +extern void notify_alive(void); +extern int set_debug_funcs(void (*slt)(const char *), + void (*ping)(void)); + + /* ------------------------------ File handling ----------------------------- */ extern char *check_prefix(char *prefix); diff --git a/meson.build b/meson.build index 6d41cccc..1cbf26fd 100644 --- a/meson.build +++ b/meson.build @@ -1,11 +1,11 @@ # Meson file for CrystFEL project('crystfel', 'c', - version: '0.11.0', + version: '0.11.1', license: 'GPL3+', meson_version: '>=0.55.0', default_options: ['buildtype=debugoptimized']) -libcrystfel_api_version = 17 +libcrystfel_api_version = 18 crystfel_rpath = '$ORIGIN/../lib64/:$ORIGIN/../lib:$ORIGIN/../../asapo-libraries/lib:$ORIGIN/../../asapo-libraries/lib64' add_project_arguments('-DHAVE_CONFIG_H', language: 'c') @@ -19,6 +19,7 @@ conf_data = configuration_data() # Mandatory dependencies cc = meson.get_compiler('c') mdep = cc.find_library('m', required: true) +rtdep = cc.find_library('rt', required: false) gsldep = dependency('gsl', required: true) pthreaddep = dependency('threads', required: true) @@ -176,7 +177,9 @@ adjust_detector = executable('adjust_detector', install_rpath: crystfel_rpath) # indexamajig -indexamajig_sources = ['src/indexamajig.c', 'src/im-sandbox.c', +indexamajig_sources = ['src/indexamajig.c', + 'src/im-sandbox.c', + 'src/im-argparse.c', 'src/process_image.c', versionc] if zmqdep.found() @@ -188,8 +191,9 @@ if asapodep.found() endif indexamajig = executable('indexamajig', indexamajig_sources, - dependencies: [mdep, libcrystfeldep, gsldep, - pthreaddep, zmqdep, asapodep, asapoproddep], + dependencies: [mdep, rtdep, libcrystfeldep, gsldep, + pthreaddep, zmqdep, asapodep, asapoproddep, + fftwdep], install: true, install_rpath: crystfel_rpath) diff --git a/relnotes-0.11.0.md b/relnotes-0.11.0.md deleted file mode 100644 index 655d23e7..00000000 --- a/relnotes-0.11.0.md +++ /dev/null @@ -1,172 +0,0 @@ -Release notes for CrystFEL version 0.11.0 -========================================= - -Copyright © 2012-2024 Deutsches Elektronen-Synchrotron DESY, - a research centre of the Helmholtz Association. - -See AUTHORS as well as the individual source code files for full contributor details. - -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/>. - - -Overview --------- - -This CrystFEL release adds a totally new system for detector geometry -refinement, as well as better support for compiling on Mac OS. The old CMake -build system has been removed in favour of Meson. The simulation tools -`pattern_sim` and `partial_sim`, which were not well maintained nor widely used, -have been removed. - -These changes are detailed below. In addition, there were many smaller fixes -and improvements. See the ChangeLog or the Git history for a comprehensive -list of all changes. - - -New system for detector geometry refinement -------------------------------------------- - -The new detector geometry system is based on the Millepede-II program, the same -software used for high energy physics experiments including the CMS. It -allows us to do a full refinement including all of the inter-dependencies -between crystal orientations, cell parameters and the detector geometry. -Previously, the refined updates to the detector positions would be biased -because they were based on biased crystal parameters, which were in turn biased -because of the incorrect detector geometry. To break this circular dependency, -we have to refine the individual parameters for each crystal at the same time -as we refine the geometry. Normally, doing this would result in an enormous -and computationally prohibitive least squares calculation, even with -specialised sparse matrix solvers. However, the Millepede algorithm reduces -the calculation to one with a sized based only on the number of geometrical -parameters for the detector. - -To get started with the new geometry refinement, read the manual for the new -`align_detector` program by running `man align_detector` or visiting -<https://gitlab.desy.de/thomas.white/crystfel/-/blob/master/doc/man/align_detector.1.md>. -An update of the geometry file will be necessary to add information about the -hierarchy of detector panels, but example files are available. The refinement -can also be performed using the CrystFEL GUI. - -A side-effect of this change is that the beam center is no longer refined -after indexing for each pattern individually. This should make the prediction -refinement more stable. - -Future CrystFEL versions will look at improving the stability and precision of -this method, in particular for three-dimensional and fine-grained refinement -tasks. - - -Homebrew formula and improved Mac OS support --------------------------------------------- - -This release adds a Homebrew formula for CrystFEL. If you're running on a Mac, -you can now install CrystFEL with a single `brew install` command. See -`INSTALL.md` for more details. Eventually we hope to submit this formula to -the main Homebrew repository to make it even easier to install. - -In addition, we are now regularly testing CrystFEL on Mac OS as part of our -continuous integration pipeline. - - -Removal of simulation tools and CMake build system --------------------------------------------------- - -Due mostly to a insufficient development resources, the simulation tools -`pattern_sim` and `partial_sim` have been removed from the CrystFEL suite -starting from this version. For more detailed discussion of the rationale for -this, see https://gitlab.desy.de/thomas.white/crystfel/-/issues/81 - -If you need quantitative image simulations, there are several better options -including [nanoBragg](https://bl831.als.lbl.gov/~jamesh/nanoBragg/). For -relative-scale geometry-only simulations of the kind done by `partial_sim`, -a better option is expected in a future CrystFEL version - watch this space! - -The CMake-based build system, which was barely maintained and lacked many -features compared to the Meson-based system, has now been removed altogether. - - -API changes ------------ - -Removed routines: -* `data_template_get_rigid_groups` -* `r_gradient` -* `x_gradient` -* `y_gradient` -* `image_hdf5_write` - -Added routines: -* `crystfel_mille_new` -* `crystfel_mille_free` -* `crystfel_mille_delete_last_record` -* `crystfel_mille_write_record` -* `mille_label` -* `mille_unlabel` -* `write_mille` -* `data_template_show_hierarchy` -* `data_template_translate_group_px` -* `data_template_translate_group_m` -* `data_template_rotate_group` -* `data_template_write_to_fh` -* `data_template_write_to_file` -* `data_template_group_info` -* `detgeom_find_panel` -* `detgeom_show_hierarchy` -* `detgeom_translate_detector_m` -* `detgeom_group_center` -* `make_panel_minvs` -* `is_cbf_file` -* `is_cbfgz_file` -* `is_hdf5_file` -* `image_data_arrays_new` -* `image_data_arrays_free` -* `image_create_dp_bad` -* `image_set_zero_data` -* `index_pattern_4` -* `r_dev` -* `fs_dev` -* `ss_dev` -* `r_gradient` -* `fs_ss_gradient` -* `read_reflections_3` -* `stream_write_data_template` -* `show_vector` -* `matrix_mult2` -* `matrix_mult3` -* `matrix_invert` -* `sq` -* `rotate2d` -* `is_dir` - -Added type definitions: -* `Mille` -* `struct dg_group_info` -* `struct detgeom_panel_group` -* `ImageDataArrays` -* `struct peak_params` -* `struct reflpeak` - -Removed type definitions: -* `struct rigid_group` -* `struct rg_collection` - -Changed structure definitions: -* `detgeom_panel` - addition of `group` -* `detgeom` - addition of `top_group` -* `image` - addition of `ida` -* `enum gparam` - changed values, moved from geometry.h to predict-refine.h - -Changed routine prototypes: -* `image_read` - addition of `ida` -* `image_read_data_block` - addition of `ida` -* `refine_prediction` - addition of `mille` and `max_mille_level` diff --git a/relnotes-0.11.1.md b/relnotes-0.11.1.md new file mode 100644 index 00000000..cc8e262f --- /dev/null +++ b/relnotes-0.11.1.md @@ -0,0 +1,127 @@ +Release notes for CrystFEL version 0.11.1 +========================================= + +Copyright © 2012-2024 Deutsches Elektronen-Synchrotron DESY, + a research centre of the Helmholtz Association. + +See AUTHORS as well as the individual source code files for full contributor details. + +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/>. + + +Overview +-------- + +This release is mainly a bug-fixing release to 0.11.0. + +These changes are detailed below. In addition, there were many smaller fixes +and improvements. See the ChangeLog or the Git history for a comprehensive +list of all changes. + +A screencast to accompany this release is available at the following locations: +<https://desy.de/~twhite/crystfel/presentations.html> +<https://vimeo.com/1017833941> + + +Improvements/bug fixes for detector geometry refinement +------------------------------------------------------- + +Improved detector geometry refinement, based on the "Millepede" algorithm +borrowed from particle physics, was introduced in the previous CrystFEL version +(0.11.0). Version 0.11.1 fixes several problems in the geometry files written +out by align_detector (and also adjust_detector), which prevented the refined +geometry files from being used in some circumstances. Corrupted "Mille" data +was also being written by indexamajig in some cases, which has been fixed. +A further problem with the graphical user interface (GUI), where the geometry +file was not correctly found by the geometry refinement task, was corrected. + +Two new options have been added to align_detector. `--out-of-plane-tilts` +allows panels to rotate outside the plane of the detector, and +`--camera-length` enables refinement of the overall camera length. +Previously, both of these motions were enabled with `--out-of-plane`. In the +new version. `--out-of-plane` only enables translations out of the detector +plane, not rotations, and does not allow the overall camera length to vary. + +If no panel groups are defined in the geometry file at all, CrystFEL will now +automatically create a single top-level group. This allows for basic (level 0) +geometry refinement without any extra work. + +As a convenience, CrystFEL now accepts panel group definitions appearing near +the start of the geometry file, before the panels have been defined. The +groups must still appear in correct hierarchy order (lowest levels first). +To maintain compatability with older CrystFEL versions, the panel groups will +be written at the end of the geometry file by align/adjust_detector. + + +GUI improvements +---------------- + +A new option "Request exclusive use of compute nodes" has been added to the +Slurm interface of the CrystFEL graphical user interface. This option should +allow for more efficient cluster usage. + +In addition, the "Jump to frame" dialogue box now gives two different ways to +specify which frame to jump to. + +The progress bar for indexing jobs run via Slurm has also been made more +realistic. + + +Julia bindings +-------------- + +An initial version of Julia bindings for "libcrystfel" has been introduced. +This will allow the various building blocks from CrystFEL (peak search, +indexing, spot prediction, reflection data management and so on) to be used +from a high-level programming language. To get started, see +<https://gitlab.desy.de/thomas.white/crystfel/-/blob/master/doc/articles/julia.rst> + + +API changes +----------- + +Removed routines: +* `crystal_copy_deep` +* `crystal_get_cell_const` +* `crystal_get_reflections` +* `crystal_get_image` +* `crystal_get_image_const` +* `crystal_set_reflections` +* `crystal_set_image` + +Added routines: +* `get_lattice_symmetry` +* `crystal_reliquish_cell` +* `image_add_crystal_refls` +* `free_reflistiterator` +* `name_equiv` +* `solve_inv` +* `set_mm_funcs` and `cfmalloc`/`cffree` etc. + +Added type definitions: +* `struct crystal_refls` + +Changed structure definitions: +* `image` - change of `crystals` to `struct crystal_refls` +* `image` - addition of `owns_peaklist` + +Changed routine prototypes: +* `update_predictions` - addition of `image` +* `predict_to_res` - addition of `image` +* `calculate_partialities` - addition of `image` +* `peakfinder8` - return type is now `ImageFeatureList` +* `search_peaks` - return type is now `ImageFeatureList` +* `search_peaks_peakfinder8` - return type is now `ImageFeatureList` +* `search_peaks_peakfinder9` - return type is now `ImageFeatureList` +* `indexing_peak_check` - addition of `peaks` +* `validate_peaks` - addition of `peaks` diff --git a/scripts/process_hkl.jl b/scripts/process_hkl.jl index 34e6104e..69e5cbc5 100644 --- a/scripts/process_hkl.jl +++ b/scripts/process_hkl.jl @@ -1,5 +1,4 @@ using CrystFEL -using LinearAlgebra let st = Stream("/home/twhite/experiments/cxidb-193/short.stream", "r"), merged = mergereflections(allcrystals(st), SymOpList("mmm")) diff --git a/src/crystfel_gui.c b/src/crystfel_gui.c index 760f8a91..53f8fb83 100644 --- a/src/crystfel_gui.c +++ b/src/crystfel_gui.c @@ -3,11 +3,11 @@ * * CrystFEL's main graphical user interface * - * Copyright © 2020-2023 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2020-2024 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2020-2023 Thomas White <taw@physics.org> + * 2020-2024 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -449,7 +449,7 @@ static gint goto_frame_sig(GtkWidget *widget, struct crystfelproject *proj) gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(stuff->framenumchk), FALSE, FALSE, 4.0); stuff->framenum = gtk_entry_new(); - snprintf(tmp, 64, "%i", proj->cur_frame); + snprintf(tmp, 64, "%i", proj->cur_frame+1); gtk_entry_set_text(GTK_ENTRY(stuff->framenum), tmp); gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(stuff->framenum), FALSE, FALSE, 4.0); @@ -685,11 +685,14 @@ static gint about_sig(GtkWidget *widget, struct crystfelproject *proj) "Richard Kirian <rkirian@asu.edu>", "Valerio Mariani <<alerio.mariani@desy.de>", "Andrew Martin <andrew.martin@desy.de>", + "Philipp Middendorf <philipp.middendorf@desy.de>", "Omri Mor <omor1@asu.edu>", "Takanori Nakane <nakane.t@gmail.com>", "Karol Nass <karol.nass@desy.de>", "Nicolas Riebesel <nicolas.riebesel@tuhh.de>", + "Silvan Schön <silvan.schoen@cfel.de>", "Mamoru Suzuki <mamoru.suzuki@protein.osaka-u.ac.jp>", + "Parthasarathy Tirumalai <parthasarathy.tirumalai@desy.de>", "Alexandra Tolstikova <alexandra.tolstikova@desy.de>", "Parker de Waal <Parker.deWaal@vai.org>", "Keitaro Yamashita <k.yamashita@spring8.or.jp>", @@ -724,7 +727,7 @@ static gint about_sig(GtkWidget *widget, struct crystfelproject *proj) gtk_about_dialog_set_version(GTK_ABOUT_DIALOG(window), crystfel_version_string()); gtk_about_dialog_set_copyright(GTK_ABOUT_DIALOG(window), - "© 2020-2023 Deutsches Elektronen-Synchrotron DESY, " + "© 2020-2024 Deutsches Elektronen-Synchrotron DESY, " "a research centre of the Helmholtz Association."); gtk_about_dialog_set_website(GTK_ABOUT_DIALOG(window), "https://www.desy.de/~twhite/crystfel"); diff --git a/src/gui_peaksearch.c b/src/gui_peaksearch.c index b67f09a8..0ab4172a 100644 --- a/src/gui_peaksearch.c +++ b/src/gui_peaksearch.c @@ -401,7 +401,7 @@ static void peaksearch_algo_changed(GtkWidget *combo, "--min-res"); add_int_param(proj->peak_params, "Maximum resolution (pixels):", &proj->peak_search_params.max_res, proj, - "--max-res"); + "--max-res (zero to disable)"); add_check_param(proj->peak_params, "Fast mode", &proj->peak_search_params.peakfinder8_fast, proj, "--peakfinder8-fast"); diff --git a/src/gui_project.c b/src/gui_project.c index 842b1b3b..a8df674d 100644 --- a/src/gui_project.c +++ b/src/gui_project.c @@ -1284,7 +1284,7 @@ int default_project(struct crystfelproject *proj) proj->peak_search_params.max_pix_count = 200; proj->peak_search_params.local_bg_radius = 3; proj->peak_search_params.min_res = 0; - proj->peak_search_params.max_res = 1200; + proj->peak_search_params.max_res = 0; proj->peak_search_params.min_snr_biggest_pix = 7.0; proj->peak_search_params.min_snr_peak_pix = 6.0; proj->peak_search_params.min_sig = 11.0; @@ -1298,6 +1298,7 @@ int default_project(struct crystfelproject *proj) proj->peak_search_params.median_filter = 0; proj->peak_search_params.check_hdf5_snr = 0; proj->peak_search_params.use_saturated = 1; + proj->peak_search_params.peakfinder8_fast = 0; proj->indexing_params.cell_file = NULL; proj->indexing_params.indexing_methods = NULL; diff --git a/src/im-argparse.c b/src/im-argparse.c new file mode 100644 index 00000000..da8321c7 --- /dev/null +++ b/src/im-argparse.c @@ -0,0 +1,1060 @@ +/* + * im-argparse.c + * + * Command line argument parsing for indexamajig + * + * Copyright © 2023-2024 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2023-2024 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 <stdarg.h> +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <unistd.h> +#include <argp.h> + +#include "version.h" +#include "im-argparse.h" + + +static void show_version(FILE *fh, struct argp_state *state) +{ + printf("CrystFEL: %s\n", crystfel_version_string()); + printf("%s\n", crystfel_licence_string()); +} + + +static DataSourceType parse_data_format(const char *str) +{ + if ( strcmp(str, "hdf5") == 0 ) return DATA_SOURCE_TYPE_HDF5; + if ( strcmp(str, "msgpack") == 0 ) return DATA_SOURCE_TYPE_MSGPACK; + if ( strcmp(str, "seedee") == 0 ) return DATA_SOURCE_TYPE_SEEDEE; + /* CBF and CBFGZ should be added here once image-cbf.c supports + * in-memory access */ + return DATA_SOURCE_TYPE_UNKNOWN; +} + + +static void add_copy_header(struct indexamajig_arguments *args, + const char *header) +{ + char **new_copy_headers = realloc(args->copy_headers, + (args->n_copy_headers+1)*sizeof(char *)); + if ( new_copy_headers == NULL ) { + ERROR("Failed to add copy header '%s'\n", header); + return; + } + + args->copy_headers = new_copy_headers; + args->copy_headers[args->n_copy_headers++] = strdup(header); +} + + +static error_t parse_arg(int key, char *arg, struct argp_state *state) +{ + float tmp; + int r; + struct indexamajig_arguments *args = state->input; + + switch ( key ) { + + case ARGP_KEY_INIT : + state->child_inputs[0] = args->taketwo_opts_ptr; + state->child_inputs[1] = args->felix_opts_ptr; + state->child_inputs[2] = args->xgandalf_opts_ptr; + state->child_inputs[3] = args->pinkindexer_opts_ptr; + state->child_inputs[4] = args->fromfile_opts_ptr; + state->child_inputs[5] = args->asdf_opts_ptr; + break; + + case 'h' : + argp_state_help(state, stdout, ARGP_HELP_STD_HELP); + break; /* argp_state_help doesn't return */ + + case 'v' : + show_version(stdout, state); + exit(0); + + case 'i' : + args->filename = strdup(arg); + break; + + case 'o' : + args->outfile = strdup(arg); + break; + + case 'x' : + free(args->prefix); + args->prefix = strdup(arg); + break; + + case 'j' : + args->n_proc = atoi(arg); + break; + + case 'g' : + args->geom_filename = arg; + break; + + case 201 : + args->basename = 1; + break; + + case 202 : + args->check_prefix = 0; + break; + + case 203 : + if ( sscanf(arg, "%f", &tmp) != 1 ) { + ERROR("Invalid value for --highres\n"); + return EINVAL; + } + args->iargs.highres = 1.0 / (tmp/1e10); /* A -> m^-1 */ + break; + + case 204 : + args->profile = 1; + break; + + case 205 : + args->temp_location = strdup(arg); + break; + + case 206 : + if (sscanf(arg, "%d", &args->iargs.wait_for_file) != 1) + { + ERROR("Invalid value for --wait-for-file\n"); + return EINVAL; + } + break; + + case 207 : + args->zmq_params.addr = strdup(arg); + break; + + case 208 : + args->iargs.no_image_data = 1; + break; + + case 209 : + ERROR("--spectrum-filename is no longer used.\n"); + return 1; + + case 210 : + args->iargs.no_mask_data = 1; + break; + + case 211 : + if ( args->zmq_params.n_subscriptions == 256 ) { + ERROR("Too many ZMQ subscriptions.\n"); + return 1; + } + args->zmq_params.subscriptions[args->zmq_params.n_subscriptions++] = strdup(arg); + break; + + case 212 : + args->zmq_params.request = strdup(arg); + break; + + case 213 : + args->asapo_params.endpoint = strdup(arg); + break; + + case 214 : + args->asapo_params.token = strdup(arg); + break; + + case 215 : + args->asapo_params.beamtime = strdup(arg); + break; + + case 217 : + args->asapo_params.group_id = strdup(arg); + break; + + case 218 : + args->asapo_params.source = strdup(arg); + break; + + case 219 : + args->iargs.data_format = parse_data_format(arg); + if ( args->iargs.data_format == DATA_SOURCE_TYPE_UNKNOWN ) { + ERROR("Unrecognised data format '%s'\n", arg); + return EINVAL; + } + break; + + case 220 : + args->asapo_params.stream = strdup(arg); + break; + + case 221 : + args->asapo_params.wait_for_stream = 1; + break; + + case 222 : + args->asapo_params.write_output_stream = 1; + break; + + case 223 : + args->cpu_pin = 1; + break; + + case 224 : + if (sscanf(arg, "%d", &args->no_data_timeout) != 1) + { + ERROR("Invalid value for --no-data-timeout\n"); + return EINVAL; + } + break; + + case 225 : + if (sscanf(arg, "%d", &args->asapo_params.consumer_timeout_ms) != 1) + { + ERROR("Invalid value for --asapo-consumer-timeout\n"); + return EINVAL; + } + break; + + case 226 : + args->asapo_params.use_ack = 1; + break; + + /* ---------- Peak search ---------- */ + + case 't' : + args->iargs.peak_search.threshold = strtof(arg, NULL); + break; + + case 301 : + args->iargs.peak_search.method = parse_peaksearch(arg); + if ( args->iargs.peak_search.method == PEAK_ERROR ) { + ERROR("Unrecognised peak detection method '%s'\n", arg); + return EINVAL; + } + break; + + case 302 : + r = sscanf(arg, "%f,%f,%f", &args->iargs.peak_search.pk_inn, + &args->iargs.peak_search.pk_mid, &args->iargs.peak_search.pk_out); + if ( (r != 3) || (args->iargs.peak_search.pk_inn < 0) ) { + ERROR("Invalid parameters for '--peak-radius'\n"); + return EINVAL; + } + break; + + case 303 : + if (sscanf(arg, "%d", &args->iargs.min_peaks) != 1) + { + ERROR("Invalid value for --min-peaks\n"); + return EINVAL; + } + break; + + case 304 : + ERROR("The option --hdf5-peaks is no longer used.\n"); + ERROR("Set the peak path in the geometry file.\n"); + break; + + case 305 : + if (sscanf(arg, "%d", &args->iargs.peak_search.median_filter) != 1) + { + ERROR("Invalid value for --median-filter\n"); + return EINVAL; + } + break; + + case 306 : + args->iargs.peak_search.noisefilter = 1; + break; + + case 307 : + if (sscanf(arg, "%f", &args->iargs.peak_search.min_sq_gradient) != 1) + { + ERROR("Invalid value for --min-squared-gradient\n"); + return EINVAL; + } + break; + + case 308 : + if (sscanf(arg, "%f", &args->iargs.peak_search.min_snr) != 1) + { + ERROR("Invalid value for --min-snr\n"); + return EINVAL; + } + break; + + case 309 : + if (sscanf(arg, "%d", &args->iargs.peak_search.min_pix_count) != 1) + { + ERROR("Invalid value for --min-pix-count\n"); + return EINVAL; + } + break; + + case 310 : + if (sscanf(arg, "%d", &args->iargs.peak_search.max_pix_count) != 1) + { + ERROR("Invalid value for --max-pix-count\n"); + return EINVAL; + } + break; + + case 311 : + if (sscanf(arg, "%d", &args->iargs.peak_search.local_bg_radius) != 1) + { + ERROR("Invalid value for --local-bg-radius\n"); + return EINVAL; + } + break; + + case 312 : + if (sscanf(arg, "%d", &args->iargs.peak_search.min_res) != 1) + { + ERROR("Invalid value for --min-res\n"); + return EINVAL; + } + break; + + case 313 : + if (sscanf(arg, "%d", &args->iargs.peak_search.max_res) != 1) + { + ERROR("Invalid value for --max-res\n"); + return EINVAL; + } + break; + + case 314 : + if (sscanf(arg, "%f", &args->iargs.peak_search.min_snr_biggest_pix) != 1) + { + ERROR("Invalid value for --max-snr-biggest-pix\n"); + return EINVAL; + } + break; + + case 315 : + if (sscanf(arg, "%f", &args->iargs.peak_search.min_snr_peak_pix) != 1) + { + ERROR("Invalid value for --max-snr-peak-pix\n"); + return EINVAL; + } + break; + + case 316 : + if (sscanf(arg, "%f", &args->iargs.peak_search.min_sig) != 1) + { + ERROR("Invalid value for --max-ssig\n"); + return EINVAL; + } + break; + + case 317 : + if (sscanf(arg, "%f", &args->iargs.peak_search.min_peak_over_neighbour) != 1) + { + ERROR("Invalid value for --max-peak-over-neighbour\n"); + return EINVAL; + } + break; + + case 318 : + args->iargs.peak_search.use_saturated = 0; + break; + + case 319 : + args->iargs.peak_search.revalidate = 0; + break; + + case 320 : + args->iargs.peak_search.half_pixel_shift = 0; + break; + + case 321 : + args->iargs.peak_search.check_hdf5_snr = 1; + break; + + case 322: + args->iargs.peak_search.peakfinder8_fast = 1; + break; + + /* ---------- Indexing ---------- */ + + case 400 : + case 'z' : + args->indm_str = strdup(arg); + break; + + case 'p' : + args->cellfile = strdup(arg); + break; + + case 401 : + /* Values in 'tols' are in frac (not %) and rad + * Conversion happens a few lines below */ + r = sscanf(arg, "%f,%f,%f,%f,%f,%f", + &args->iargs.tols[0], &args->iargs.tols[1], &args->iargs.tols[2], + &args->iargs.tols[3], &args->iargs.tols[4], &args->iargs.tols[5]); + if ( r != 6 ) { + /* Try old format */ + r = sscanf(arg, "%f,%f,%f,%f", + &args->iargs.tols[0], &args->iargs.tols[1], + &args->iargs.tols[2], &args->iargs.tols[3]); + if ( r != 4 ) { + ERROR("Invalid parameters for '--tolerance'\n"); + return EINVAL; + } + args->iargs.tols[4] = args->iargs.tols[3]; + args->iargs.tols[5] = args->iargs.tols[3]; + } + + /* Percent to fraction */ + args->iargs.tols[0] /= 100.0; + args->iargs.tols[1] /= 100.0; + args->iargs.tols[2] /= 100.0; + args->iargs.tols[3] = deg2rad(args->iargs.tols[3]); + args->iargs.tols[4] = deg2rad(args->iargs.tols[4]); + args->iargs.tols[5] = deg2rad(args->iargs.tols[5]); + break; + + case 402 : + args->if_checkcell = 0; + break; + + case 403 : + args->if_checkcell = 1; /* This is the default */ + break; + + case 404 : + args->if_multi = 1; + break; + + case 405 : + args->if_multi = 0; /* This is the default */ + break; + + case 406 : + args->if_retry = 0; + break; + + case 407 : + args->if_retry = 1; /* This is the default */ + break; + + case 408 : + args->if_refine = 0; + break; + + case 409 : + args->if_refine = 1; /* This is the default */ + break; + + case 410 : + args->if_peaks = 0; + break; + + case 411 : + args->if_peaks = 1; /* This is the default */ + break; + + case 412 : + ERROR("The option --no-cell-combinations is no longer used.\n"); + /* .. but we can still carry on. Results will probably be + * better than the user expected. */ + break; + + case 413 : + if (sscanf(arg, "%f", &args->iargs.wavelength_estimate) != 1) + { + ERROR("Invalid value for --wavelength-estimate\n"); + return EINVAL; + } + break; + + case 414 : + if (sscanf(arg, "%d", &args->iargs.n_threads) != 1) + { + ERROR("Invalid value for --max-indexer-threads\n"); + return EINVAL; + } + break; + + case 415 : + if (sscanf(arg, "%f", &args->iargs.clen_estimate) != 1) + { + ERROR("Invalid value for --camera-length-estimate\n"); + return EINVAL; + } + break; + + case 416 : + args->iargs.mille = 1; + break; + + case 417 : + free(args->milledir); + args->milledir = strdup(arg); + break; + + case 418 : + if (sscanf(arg, "%d", &args->iargs.max_mille_level) != 1) + { + ERROR("Invalid value for --max-mille-level\n"); + return EINVAL; + } + if ( args->iargs.max_mille_level < 0 ) { + ERROR("Invalid value for --max-mille-level\n"); + return EINVAL; + } + break; + + case 419 : + free(args->millefile); + args->millefile = strdup(arg); + break; + + /* ---------- Integration ---------- */ + + case 501 : + args->iargs.int_meth = integration_method(arg, &r); + if ( r ) { + ERROR("Invalid integration method '%s'\n", arg); + return EINVAL; + } + break; + + case 502 : + if ( sscanf(arg, "%f", &args->iargs.fix_profile_r) != 1 ) { + ERROR("Invalid value for --fix-profile-radius\n"); + return EINVAL; + } + break; + + case 503 : + ERROR("The option --fix-bandwidth is no longer used.\n"); + ERROR("Set the bandwidth in the geometry file instead.\n"); + break; + + case 504 : + if ( sscanf(arg, "%f", &args->iargs.fix_divergence) != 1 ) { + ERROR("Invalid value for --fix-divergence\n"); + return EINVAL; + } + break; + + case 505 : + r = sscanf(arg, "%f,%f,%f", &args->iargs.ir_inn, + &args->iargs.ir_mid, &args->iargs.ir_out); + if ( r != 3 ) { + ERROR("Invalid parameters for '--int-radius'\n"); + return EINVAL; + } + break; + + case 506 : + if ( strcmp(arg, "random") == 0 ) { + args->iargs.int_diag = INTDIAG_RANDOM; + } + + if ( strcmp(arg, "all") == 0 ) { + args->iargs.int_diag = INTDIAG_ALL; + } + + if ( strcmp(arg, "negative") == 0 ) { + args->iargs.int_diag = INTDIAG_NEGATIVE; + } + + if ( strcmp(arg, "implausible") == 0 ) { + args->iargs.int_diag = INTDIAG_IMPLAUSIBLE; + } + + if ( strcmp(arg, "strong") == 0 ) { + args->iargs.int_diag = INTDIAG_STRONG; + } + + r = sscanf(arg, "%i,%i,%i", &args->iargs.int_diag_h, + &args->iargs.int_diag_k, &args->iargs.int_diag_l); + if ( r == 3 ) { + args->iargs.int_diag = INTDIAG_INDICES; + } + + if ( (args->iargs.int_diag == INTDIAG_NONE) + && (strcmp(arg, "none") != 0) ) + { + ERROR("Invalid value for --int-diag.\n"); + return EINVAL; + } + + break; + + case 507 : + if ( sscanf(arg, "%f", &args->iargs.push_res) != 1 ) { + ERROR("Invalid value for --push-res\n"); + return EINVAL; + } + args->iargs.push_res *= 1e9; /* nm^-1 -> m^-1 */ + break; + + case 508 : + args->iargs.overpredict = 1; + break; + + case 509 : + args->iargs.cell_params_only = 1; + break; + + /* ---------- Output ---------- */ + + case 601 : + args->iargs.stream_nonhits = 0; + break; + + case 602 : + add_copy_header(args, arg); + break; + + case 603 : + args->iargs.stream_flags = CLEAR_BIT(args->iargs.stream_flags, + STREAM_PEAKS); + break; + + case 604 : + args->iargs.stream_flags = CLEAR_BIT(args->iargs.stream_flags, + STREAM_REFLECTIONS); + break; + + case 605 : + if ( sscanf(arg, "%d", &args->serial_start) != 1 ) { + ERROR("Invalid value for --serial-start\n"); + return EINVAL; + } + break; + + case 606 : + args->harvest_file = strdup(arg); + break; + + /* ---------- Secret muti-processing stuff ---------- */ + case 703 : + args->fd_stream = atoi(arg); + break; + + case 704 : + args->shm_name = strdup(arg); + break; + + case 705 : + args->fd_mille = atoi(arg); + break; + + case 706 : + args->worker_id = atoi(arg); + args->worker = 1; + break; + + case 707 : + args->worker_tmpdir = strdup(arg); + break; + + case 708 : + args->queue_sem = strdup(arg); + break; + + default : + return ARGP_ERR_UNKNOWN; + + } + + return 0; +} + +struct taketwo_options *taketwo_opts = NULL; +struct felix_options *felix_opts = NULL; +struct xgandalf_options *xgandalf_opts = NULL; +struct pinkindexer_options *pinkindexer_opts = NULL; +struct fromfile_options *fromfile_opts = NULL; +struct asdf_options *asdf_opts = NULL; + +struct indexamajig_arguments *parse_indexamajig_args(int argc, char *argv[]) +{ + struct indexamajig_arguments *args; + + args = malloc(sizeof(struct indexamajig_arguments)); + if ( args == NULL ) return NULL; + + /* Defaults for "top level" arguments */ + args->filename = NULL; + args->geom_filename = NULL; + args->outfile = NULL; + args->temp_location = strdup("."); + args->prefix = strdup(""); + args->check_prefix = 1; + args->n_proc = 1; + args->cellfile = NULL; + args->indm_str = NULL; + args->basename = 0; + args->zmq_params.addr = NULL; + args->zmq_params.request = NULL; + args->zmq_params.n_subscriptions = 0; + args->asapo_params.endpoint = NULL; + args->asapo_params.token = NULL; + args->asapo_params.beamtime = NULL; + args->asapo_params.group_id = NULL; + args->asapo_params.source = NULL; + args->asapo_params.stream = NULL; + args->asapo_params.wait_for_stream = 0; + args->asapo_params.write_output_stream = 0; + args->asapo_params.consumer_timeout_ms = 500; + args->asapo_params.use_ack = 0; + args->cpu_pin = 0; + args->serial_start = 1; + args->if_peaks = 1; + args->if_multi = 0; + args->if_retry = 1; + args->if_refine = 1; + args->if_checkcell = 1; + args->profile = 0; + args->no_data_timeout = 60; + args->copy_headers = NULL; + args->n_copy_headers = 0; + args->harvest_file = NULL; + args->taketwo_opts_ptr = &taketwo_opts; + args->felix_opts_ptr = &felix_opts; + args->xgandalf_opts_ptr = &xgandalf_opts; + args->pinkindexer_opts_ptr = &pinkindexer_opts; + args->fromfile_opts_ptr = &fromfile_opts; + args->asdf_opts_ptr = &asdf_opts; + args->worker = 0; + args->fd_stream = 0; + args->fd_mille = 0; + args->milledir = strdup("."); + args->millefile = strdup("mille-data.bin"); + args->worker_tmpdir = NULL; + args->queue_sem = NULL; + args->shm_name = NULL; + + /* Defaults for process_image arguments */ + args->iargs.cell = NULL; + args->iargs.peak_search.noisefilter = 0; + args->iargs.peak_search.median_filter = 0; + args->iargs.tols[0] = 0.05; /* frac (not %) */ + args->iargs.tols[1] = 0.05; /* frac (not %) */ + args->iargs.tols[2] = 0.05; /* frac (not %) */ + args->iargs.tols[3] = deg2rad(1.5); /* radians */ + args->iargs.tols[4] = deg2rad(1.5); /* radians */ + args->iargs.tols[5] = deg2rad(1.5); /* radians */ + args->iargs.peak_search.threshold = 800.0; + args->iargs.peak_search.min_sq_gradient = 100000.0; + args->iargs.peak_search.min_snr = 5.0; + args->iargs.peak_search.min_pix_count = 2; + args->iargs.peak_search.max_pix_count = 200; + args->iargs.peak_search.min_res = 0; + args->iargs.peak_search.max_res = 0; + args->iargs.peak_search.local_bg_radius = 3; + args->iargs.peak_search.min_snr_biggest_pix = 7.0; /* peak finder 9 */ + args->iargs.peak_search.min_snr_peak_pix = 6.0; + args->iargs.peak_search.min_sig = 11.0; + args->iargs.peak_search.min_peak_over_neighbour = -INFINITY; + args->iargs.peak_search.check_hdf5_snr = 0; + args->iargs.peak_search.peakfinder8_fast = 0; + args->iargs.pf_private = NULL; + args->iargs.dtempl = NULL; + args->iargs.peak_search.method = PEAK_ZAEF; + args->iargs.peak_search.half_pixel_shift = 1; + args->iargs.peak_search.pk_inn = -1.0; + args->iargs.peak_search.pk_mid = -1.0; + args->iargs.peak_search.pk_out = -1.0; + args->iargs.ir_inn = -1.0; + args->iargs.ir_mid = -1.0; + args->iargs.ir_out = -1.0; + args->iargs.peak_search.use_saturated = 1; + args->iargs.peak_search.revalidate = 1; + args->iargs.stream_flags = STREAM_PEAKS | STREAM_REFLECTIONS; + args->iargs.stream_nonhits = 1; + args->iargs.int_diag = INTDIAG_NONE; + args->iargs.min_peaks = 0; + args->iargs.overpredict = 0; + args->iargs.cell_params_only = 0; + args->iargs.wait_for_file = 0; + args->iargs.ipriv = NULL; /* No default */ + args->iargs.int_meth = integration_method("rings-nocen-nosat-nograd", NULL); + args->iargs.push_res = +INFINITY; + args->iargs.highres = +INFINITY; + args->iargs.fix_profile_r = -1.0; + args->iargs.fix_divergence = -1.0; + args->iargs.no_image_data = 0; + args->iargs.no_mask_data = 0; + args->iargs.wavelength_estimate = NAN; + args->iargs.clen_estimate = NAN; + args->iargs.n_threads = 1; + args->iargs.data_format = DATA_SOURCE_TYPE_UNKNOWN; + args->iargs.mille = 0; + args->iargs.max_mille_level = 99; + + argp_program_version_hook = show_version; + + static char doc[] = "Index and integrate snapshot diffraction images.\v" + "For more information including a tutorial, visit " + "https://www.desy.de/~twhite/crystfel"; + + static struct argp_option options[] = { + + {NULL, 0, 0, OPTION_DOC, "Basic options:", 2}, + + {NULL, 'h', NULL, OPTION_HIDDEN, NULL}, + {NULL, 'v', NULL, OPTION_HIDDEN, NULL}, + + {"input", 'i', "infile", 0, "List of input image filenames"}, + {"output", 'o', "filename.stream", 0, "Output stream filename"}, + {"geometry",'g', "experiment.geom", 0, "Detector geometry filename"}, + {"prefix", 'x', "/path/to/images/", OPTION_NO_USAGE, "Prefix filenames from input " + "file"}, + {NULL, 'j', "nproc", 0, "Run this many analyses in parallel, default 1"}, + {"basename", 201, NULL, OPTION_NO_USAGE, "Remove director parts from the " + "filenames"}, + {"no-check-prefix", 202, NULL, OPTION_NO_USAGE, "Don't attempt to correct the " + "--prefix"}, + {"highres", 203, "res", OPTION_NO_USAGE, "Absolute resolution cutoff in Angstroms"}, + {"profile", 204, NULL, OPTION_NO_USAGE, "Show timing data for performance " + "monitoring"}, + {"temp-dir", 205, "path", OPTION_NO_USAGE, "Location for temporary folder"}, + {"wait-for-file", 206, "seconds", OPTION_NO_USAGE, "Wait for each file before " + "processing"}, + {"zmq-input", 207, "addr", OPTION_NO_USAGE, "Receive data over ZeroMQ from " + "this location"}, + {"no-image-data", 208, NULL, OPTION_NO_USAGE, "Do not load image data"}, + {"spectrum-file", 209, "fn", OPTION_NO_USAGE | OPTION_HIDDEN, + "File containing radiation spectrum"}, + {"no-mask-data", 210, NULL, OPTION_NO_USAGE, "Do not load mask data"}, + {"zmq-subscribe", 211, "tag", OPTION_NO_USAGE, "Subscribe to ZMQ message" + "type"}, + {"zmq-request", 212, "str", OPTION_NO_USAGE, "Request messages using" + "this string."}, + {"asapo-endpoint", 213, "str", OPTION_NO_USAGE, "ASAP::O endpoint"}, + {"asapo-token", 214, "str", OPTION_NO_USAGE, "ASAP::O token"}, + {"asapo-beamtime", 215, "str", OPTION_NO_USAGE, "ASAP::O beamtime ID"}, + {"asapo-group", 217, "str", OPTION_NO_USAGE, "ASAP::O group ID"}, + {"asapo-source", 218, "str", OPTION_NO_USAGE, "ASAP::O data source"}, + {"data-format", 219, "str", OPTION_NO_USAGE, "Streamed data format"}, + {"asapo-stream", 220, "str", OPTION_NO_USAGE, "ASAP::O stream name"}, + {"asapo-wait-for-stream", 221, NULL, OPTION_NO_USAGE, + "Wait for ASAP::O stream to appear"}, + {"asapo-output-stream", 222, NULL, OPTION_NO_USAGE, + "Create an ASAP::O hits-only stream"}, + {"cpu-pin", 223, NULL, OPTION_NO_USAGE, "Pin worker processes to CPUs"}, + {"no-data-timeout", 224, "s", OPTION_NO_USAGE, + "Shut down after this many seconds without ASAP::O data"}, + {"asapo-consumer-timeout", 225, "ms", OPTION_NO_USAGE, + "ASAP::O get_next timeout for one frame (milliseconds)"}, + {"asapo-acks", 226, NULL, OPTION_NO_USAGE, "Use ASAP::O acknowledgements"}, + + {NULL, 0, 0, OPTION_DOC, "Peak search options:", 3}, + {"peaks", 301, "method", 0, "Peak search method. Default: zaef"}, + {"peak-radius", 302, "r1,r2,r3", OPTION_NO_USAGE, "Radii for peak search"}, + {"min-peaks", 303, "n", OPTION_NO_USAGE, "Minimum number of peaks for indexing"}, + {"hdf5-peaks", 304, "p", OPTION_HIDDEN, "Location of peak table in HDF5 file"}, + {"median-filter", 305, "n", OPTION_NO_USAGE, "Apply median filter to image data"}, + {"filter-noise", 306, NULL, OPTION_NO_USAGE, "Apply noise filter to image data"}, + {"threshold", 't', "adu", OPTION_NO_USAGE, "Threshold for peak detection " + "(zaef only, default 800)"}, + {"min-squared-gradient", 307, "n", OPTION_NO_USAGE, "Minimum squared gradient " + "(zaef only, default 100000)"}, + {"min-gradient", 307, "n", OPTION_ALIAS | OPTION_HIDDEN, NULL}, + {"min-snr", 308, "n", OPTION_NO_USAGE, "Minimum signal/noise ratio for peaks " + "(zaef,peakfinder8,peakfinder9 only, default 5)"}, + {"min-pix-count", 309, "n", OPTION_NO_USAGE, "Minimum number of pixels per peak " + "(peakfinder8 only, default 2)"}, + {"max-pix-count", 310, "n", OPTION_NO_USAGE, "Maximum number of pixels per peak " + "(peakfinder8 only, default 2)"}, + {"local-bg-radius", 311, "n", OPTION_NO_USAGE, "Radius (pixels) for local " + "background estimation (peakfinder8/9 only, default 3)"}, + {"min-res", 312, "n", OPTION_NO_USAGE, "Minimum resoultion (pixels) for peak " + "search (peakfinder8 only, default 0)"}, + {"max-res", 313, "n", OPTION_NO_USAGE, "Maximum resoultion (pixels) for peak " + "search (peakfinder8 only, default 1200)"}, + {"min-snr-biggest-pix", 314, "n", OPTION_NO_USAGE, "Minimum SNR of the biggest " + "pixel in the peak (peakfinder9 only)"}, + {"min-snr-peak-pix", 315, "n", OPTION_NO_USAGE, "Minimum SNR of peak pixel " + "(peakfinder9 only)"}, + {"min-sig", 316, "n", OPTION_NO_USAGE, "Minimum standard deviation of the " + "background (peakfinder9 only)"}, + {"min-peak-over-neighbour", 317, "n", OPTION_NO_USAGE, "Minimum difference between " + "highest pixel and neighbours (peakfinder9 only, just for speed)"}, + {"no-use-saturated", 318, NULL, OPTION_NO_USAGE, "Reject saturated peaks"}, + {"no-revalidate", 319, NULL, OPTION_NO_USAGE, "Don't re-integrate and check HDF5 " + "or MsgPack peaks"}, + {"no-half-pixel-shift", 320, NULL, OPTION_NO_USAGE, "Don't offset HDF5 peak " + "locations by 0.5 pixels"}, + {"check-hdf5-snr", 321, NULL, OPTION_NO_USAGE, "Check SNR for peaks from HDF5, " + "CXI or MsgPack (see --min-snr)"}, + {"peakfinder8-fast", 322, NULL, OPTION_NO_USAGE, "peakfinder8 fast execution"}, + + {NULL, 0, 0, OPTION_DOC, "Indexing options:", 4}, + {"indexing", 400, "method", 0, "List of indexing methods"}, + {NULL, 'z', "method", OPTION_HIDDEN | OPTION_ALIAS, NULL}, + {"pdb", 'p', "parameters.cell", 0, "PDB or CrystFEL Unit Cell File"}, + {"tolerance", 401, "a,b,c,al,be,ga", OPTION_NO_USAGE, "Tolerances for cell " + "comparison in percent and degrees, default 5,5,5,1.5,1.5,1.5"}, + {"no-check-cell", 402, NULL, OPTION_NO_USAGE, "Don't check cell parameters " + "against target cell"}, + {"check-cell", 403, NULL, OPTION_HIDDEN, NULL}, + {"multi", 404, NULL, OPTION_NO_USAGE, "Repeat indexing to index multiple hits"}, + {"no-multi", 405, NULL, OPTION_HIDDEN, NULL}, + {"no-retry", 406, NULL, OPTION_NO_USAGE, "Don't repeat indexing to increase " + "indexing rate"}, + {"retry", 407, NULL, OPTION_HIDDEN, NULL}, + {"no-refine", 408, NULL, OPTION_NO_USAGE, "Skip prediction refinement"}, + {"refine", 409, NULL, OPTION_HIDDEN, NULL}, + {"no-check-peaks", 410, NULL, OPTION_NO_USAGE, "Don't check that most peaks can be " + "accounted for by the indexing solution"}, + {"check-peaks", 411, NULL, OPTION_HIDDEN, NULL}, + {"no-cell-combinations", 412, NULL, OPTION_HIDDEN, NULL}, + {"wavelength-estimate", 413, "metres", 0, + "Estimate of the incident radiation wavelength, in metres."}, + {"max-indexer-threads", 414, "n", 0, + "Maximum number of threads allowed for indexing engines."}, + {"camera-length-estimate", 415, "metres", 0, + "Estimate of the camera length, in metres."}, + {"mille", 416, NULL, 0, + "Generate data for detector geometry refinement using Millepede"}, + {"mille-dir", 417, "dirname", OPTION_HIDDEN, "Save Millepede data in folder"}, + {"max-mille-level", 418, "n", 0, "Maximum geometry refinement level"}, + {"mille-file", 419, "filename", 0, "Filename for Millepede data (default mille-data.bin)"}, + + {NULL, 0, 0, OPTION_DOC, "Integration options:", 5}, + {"integration", 501, "method", OPTION_NO_USAGE, "Integration method"}, + {"fix-profile-radius", 502, "r", OPTION_NO_USAGE, "Fix profile radius for spot " + "prediction, instead of automatically determining"}, + {"fix-bandwidth", 503, "bw", OPTION_NO_USAGE, "Set the bandwidth for spot " + "prediction"}, + {"fix-divergence", 504, "deg", OPTION_NO_USAGE, "Set the divergence (full angle) " + "for spot prediction"}, + {"int-radius", 505, "r1,r2,r3", 0, "Set the integration radii (inner,mid,outer)"}, + {"int-diag", 506, "condition", 0, "Show debugging information about reflections"}, + {"push-res", 507, "dist", 0, "Integrate higher than apparent resolution cutoff (m^-1)"}, + {"overpredict", 508, NULL, 0, "Over-predict reflections"}, + {"cell-parameters-only", 509, NULL, 0, "Don't predict reflections at all"}, + + {NULL, 0, 0, OPTION_DOC, "Output options:", 6}, + {"no-non-hits-in-stream", 601, NULL, OPTION_NO_USAGE, "Don't include non-hits in " + "stream (see --min-peaks)"}, + {"copy-hdf5-field", 602, "f", OPTION_HIDDEN, NULL}, + {"copy-header", 602, "f", OPTION_NO_USAGE, "Put the value of this image header " + "field into the stream"}, + {"no-peaks-in-stream", 603, NULL, OPTION_NO_USAGE, "Don't put peak search results " + "in stream"}, + {"no-refls-in-stream", 604, NULL, OPTION_NO_USAGE, "Don't put integration results " + "in stream"}, + {"serial-start", 605, "n", OPTION_NO_USAGE, "Start the serial numbers in the stream " + "here"}, + {"harvest-file", 606, "filename", OPTION_NO_USAGE, "Write the actual parameters " + "used in JSON format"}, + + {"fd-stream", 703, "fd", OPTION_HIDDEN, "File descriptor for input"}, + {"shm-name", 704, "n", OPTION_HIDDEN, "SHM name for queue structure"}, + {"fd-mille", 705, "fd", OPTION_HIDDEN, "File descriptor for Mille data"}, + {"worker-id", 706, "id", OPTION_HIDDEN, "Worker process number"}, + {"worker-tmpdir", 707, "dir", OPTION_HIDDEN, "Worker temporary dir"}, + {"queue-semaphore", 708, "sem", OPTION_HIDDEN, "Queue semaphore name"}, + + {NULL, 0, 0, OPTION_DOC, "More information:", 99}, + + {0} + + }; + + static struct argp_child argp_children[] = { + {&taketwo_argp, 0, NULL, -2}, + {&felix_argp, 0, NULL, -2}, + {&xgandalf_argp, 0, NULL, -2}, + {&pinkIndexer_argp, 0, NULL, -2}, + {&fromfile_argp, 0, NULL, -2}, + {&asdf_argp, 0, NULL, -2}, + {0} + }; + + #if defined(__APPLE__) + int q; + for ( q=0; q<argc; q++ ) { + if ( strcmp(argv[q], "--help") == 0 ) { + fprintf(stderr, "\n" + "WARNING: 'indexamajig --help' crashes on some Mac OS versions.\n" + "This is a known problem, and is due to a bug in external library " + "code.\n\n"); + fflush(stderr); + break; + } + } + #endif + + static struct argp argp = { options, parse_arg, NULL, doc, + argp_children, NULL, NULL }; + if ( argp_parse(&argp, argc, argv, 0, NULL, args) ) return NULL; + + return args; +} + + +void cleanup_indexamajig_args(struct indexamajig_arguments *args) +{ + int i; + + free(args->milledir); + free(args->millefile); + free(args->prefix); + free(args->temp_location); + free(args->harvest_file); + free(*args->taketwo_opts_ptr); + free(*args->felix_opts_ptr); + free(*args->xgandalf_opts_ptr); + free(*args->pinkindexer_opts_ptr); + free(*args->fromfile_opts_ptr); + free(*args->asdf_opts_ptr); + free(args->filename); + free(args->outfile); + free(args->zmq_params.addr); + free(args->zmq_params.request); + for ( i=0; i<args->zmq_params.n_subscriptions; i++ ) { + free(args->zmq_params.subscriptions[i]); + } + free(args->asapo_params.endpoint); + free(args->asapo_params.token); + free(args->asapo_params.beamtime); + free(args->asapo_params.endpoint); + free(args->asapo_params.group_id); + free(args->asapo_params.source); + free(args->asapo_params.stream); + free(args->indm_str); + free(args->cellfile); + free(args->worker_tmpdir); + free(args->queue_sem); + free(args->shm_name); + for ( i=0; i<args->n_copy_headers; i++ ) { + free(args->copy_headers[i]); + } + free(args->copy_headers); + free(args); +} diff --git a/src/im-argparse.h b/src/im-argparse.h new file mode 100644 index 00000000..f6cbee34 --- /dev/null +++ b/src/im-argparse.h @@ -0,0 +1,86 @@ +/* + * im-argparse.h + * + * Command line argument parsing for indexamajig + * + * Copyright © 2023-2024 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2023-2024 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 IM_ARGPARSE_H +#define IM_ARGPARSE_H + +struct indexamajig_arguments; + +#include <index.h> + +#include "process_image.h" + +struct indexamajig_arguments +{ + struct index_args iargs; /* These are the options that will be + * given to process_image */ + char *filename; + char *geom_filename; + char *outfile; + char *prefix; + int check_prefix; + int n_proc; + char *cellfile; + char *indm_str; + int basename; + struct im_zmq_params zmq_params; + struct im_asapo_params asapo_params; + int serial_start; + char *temp_location; + int if_refine; + int if_checkcell; + int if_peaks; + int if_multi; + int if_retry; + int profile; /* Whether to do wall-clock time profiling */ + int no_data_timeout; + char **copy_headers; + int n_copy_headers; + char *harvest_file; + char *milledir; + char *millefile; + int cpu_pin; + int worker; + int worker_id; + char *worker_tmpdir; + int fd_stream; + int fd_mille; + char *queue_sem; + char *shm_name; + + struct taketwo_options **taketwo_opts_ptr; + struct felix_options **felix_opts_ptr; + struct xgandalf_options **xgandalf_opts_ptr; + struct pinkindexer_options **pinkindexer_opts_ptr; + struct fromfile_options **fromfile_opts_ptr; + struct asdf_options **asdf_opts_ptr; +}; + +extern struct indexamajig_arguments *parse_indexamajig_args(int argc, char *argv[]); +extern void cleanup_indexamajig_args(struct indexamajig_arguments *args); + +#endif /* IM_ARGPARSE_H */ diff --git a/src/im-asapo.c b/src/im-asapo.c index 9b57883c..2f483a2d 100644 --- a/src/im-asapo.c +++ b/src/im-asapo.c @@ -74,7 +74,7 @@ static void show_asapo_error(const char *msg, const AsapoErrorHandle err) } -static int create_producer(struct im_asapo *a, struct im_asapo_params *params) +static int create_producer(struct im_asapo *a, const struct im_asapo_params *params) { char *source; AsapoSourceCredentialsHandle cred; @@ -122,7 +122,7 @@ static int create_producer(struct im_asapo *a, struct im_asapo_params *params) } -struct im_asapo *im_asapo_connect(struct im_asapo_params *params) +struct im_asapo *im_asapo_connect(const struct im_asapo_params *params) { struct im_asapo *a; AsapoSourceCredentialsHandle cred; diff --git a/src/im-asapo.h b/src/im-asapo.h index 70f67fd0..87dbce42 100644 --- a/src/im-asapo.h +++ b/src/im-asapo.h @@ -52,7 +52,7 @@ struct im_asapo; #if defined(HAVE_ASAPO) -extern struct im_asapo *im_asapo_connect(struct im_asapo_params *params); +extern struct im_asapo *im_asapo_connect(const struct im_asapo_params *params); extern void im_asapo_shutdown(struct im_asapo *a); @@ -66,7 +66,7 @@ extern void im_asapo_send(struct im_asapo *a, struct image *image, int hit); #else /* defined(HAVE_ASAPO) */ -static UNUSED struct im_asapo *im_asapo_connect(struct im_asapo_params *params) +static UNUSED struct im_asapo *im_asapo_connect(const struct im_asapo_params *params) { ERROR("This installation of CrystFEL was compiled without ASAP::O support.\n"); return NULL; diff --git a/src/im-sandbox.c b/src/im-sandbox.c index fac01999..4d69ab35 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -3,13 +3,13 @@ * * Sandbox for indexing * - * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2024 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2020 Thomas White <taw@physics.org> + * 2010-2024 Thomas White <taw@physics.org> * 2014 Valerio Mariani * 2011 Richard Kirian * 2012 Lorenzo Galli @@ -65,6 +65,7 @@ #endif #include "im-sandbox.h" +#include "im-argparse.h" #include "process_image.h" #include "im-zmq.h" #include "profile.h" @@ -72,10 +73,20 @@ #include "predict-refine.h" +typedef struct +{ + int n_read; + int *fds; + void **buffers; + size_t *buffer_len; + size_t *buffer_pos; +} PipeList; + + struct sandbox { int n_processed_last_stats; - double t_last_stats; + time_t t_last_stats; /* Processing timeout in seconds. After this long without responding * to a ping, the worker will be killed. After 3 times this long @@ -84,6 +95,9 @@ struct sandbox int timeout; struct index_args *iargs; + int argc; + char **argv; + const char *probed_methods; /* Worker processes */ int n_proc; @@ -91,18 +105,20 @@ struct sandbox int *running; time_t *last_response; int last_ping[MAX_NUM_WORKERS]; + int warned_long_running[MAX_NUM_WORKERS]; int profile; /* Whether to do wall-clock time profiling */ int cpu_pin; /* Streams to read from (NB not the same indices as the above) */ - int n_read; - FILE **fhs; - int *fds; + PipeList *st_from_workers; + PipeList *mille_from_workers; int serial; struct sb_shm *shared; + char *shm_name; sem_t *queue_sem; + char *sem_name; const char *tmpdir; @@ -114,6 +130,7 @@ struct sandbox /* Final output */ Stream *stream; + FILE *mille_fh; }; struct get_pattern_ctx @@ -131,11 +148,11 @@ struct get_pattern_ctx #ifdef HAVE_CLOCK_GETTIME -static double get_monotonic_seconds() +time_t get_monotonic_seconds() { struct timespec tp; clock_gettime(CLOCK_MONOTONIC, &tp); - return tp.tv_sec + tp.tv_nsec * 1e-9; + return tp.tv_sec; } #else @@ -143,11 +160,11 @@ static double get_monotonic_seconds() /* Fallback version of the above. The time according to gettimeofday() is not * monotonic, so measuring intervals based on it will screw up if there's a * timezone change (e.g. daylight savings) while the program is running. */ -static double get_monotonic_seconds() +time_t get_monotonic_seconds() { struct timeval tp; gettimeofday(&tp, NULL); - return tp.tv_sec + tp.tv_usec * 1e-6; + return tp.tv_sec; } #endif @@ -164,6 +181,8 @@ static void check_hung_workers(struct sandbox *sb) { int i; time_t tnow = get_monotonic_seconds(); + + pthread_mutex_lock(&sb->shared->debug_lock); for ( i=0; i<sb->n_proc; i++ ) { if ( !sb->running[i] ) continue; @@ -180,7 +199,7 @@ static void check_hung_workers(struct sandbox *sb) } if ( tnow - sb->shared->time_last_start[i] > sb->timeout*3 ) { - if ( !sb->shared->warned_long_running[i] ) { + if ( !sb->warned_long_running[i] ) { STATUS("Worker %i has been working on one " "frame for more than %i seconds (just " "for info).\n", i, sb->timeout); @@ -188,11 +207,12 @@ static void check_hung_workers(struct sandbox *sb) sb->shared->last_ev[i]); STATUS("Task ID is: %s\n", sb->shared->last_task[i]); - sb->shared->warned_long_running[i] = 1; + sb->warned_long_running[i] = 1; } } } + pthread_mutex_unlock(&sb->shared->debug_lock); } @@ -280,7 +300,11 @@ static int get_pattern(struct get_pattern_ctx *gpctx, filename = read_prefixed_filename(gpctx, &evstr); /* Nothing left in file -> we're done */ - if ( filename == NULL ) return 0; + if ( filename == NULL ) { + free(gpctx->filename); + free(gpctx->events); + return 0; + } /* Does the line from the input file contain an event ID? * If so, just send it straight back. */ @@ -316,377 +340,181 @@ static int get_pattern(struct get_pattern_ctx *gpctx, } -static void shuffle_events(struct sb_shm *sb_shared) +/* Like strstr(), but with a length parameter instead of nul-termination */ +static const char *str_in_str(const char *haystack, size_t len, const char *needle) { - int i; + size_t u; + const size_t endl = strlen(needle); - for ( i=1; i<sb_shared->n_events; i++ ) { - memcpy(sb_shared->queue[i-1], sb_shared->queue[i], MAX_EV_LEN); - } - sb_shared->n_events--; -} + if ( len < endl ) return NULL; - -void set_last_task(char *lt, const char *task) -{ - if ( lt == NULL ) return; - assert(strlen(task) < MAX_TASK_LEN-1); - strcpy(lt, task); -} - - -static int run_work(const struct index_args *iargs, Stream *st, - int cookie, const char *tmpdir, struct sandbox *sb) -{ - int allDone = 0; - struct im_zmq *zmqstuff = NULL; - struct im_asapo *asapostuff = NULL; - Mille *mille; - ImageDataArrays *ida; - - if ( sb->profile ) { - profile_init(); - } - - /* Connect via ZMQ */ - if ( sb->zmq_params != NULL ) { - zmqstuff = im_zmq_connect(sb->zmq_params); - if ( zmqstuff == NULL ) { - ERROR("ZMQ setup failed.\n"); - return 1; - } - } - - if ( sb->asapo_params != NULL ) { - asapostuff = im_asapo_connect(sb->asapo_params); - if ( asapostuff == NULL ) { - ERROR("ASAP::O setup failed.\n"); - sb->shared->should_shutdown = 1; - return 1; - } - } - - mille = NULL; - if ( iargs->mille ) { - char tmp[1024]; - snprintf(tmp, 1024, "%s/mille-data-%i.bin", iargs->milledir, cookie); - mille = crystfel_mille_new(tmp); - } - - ida = image_data_arrays_new(); - - while ( !allDone ) { - - struct pattern_args pargs; - int ser; - char *line; - size_t len; - int i; - char *event_str = NULL; - char *ser_str = NULL; + for ( u=0; u<len-endl; u++ ) { + size_t v; int ok = 1; - - /* Wait until an event is ready */ - sb->shared->pings[cookie]++; - set_last_task(sb->shared->last_task[cookie], "wait_event"); - profile_start("wait-queue-semaphore"); - if ( sem_wait(sb->queue_sem) != 0 ) { - ERROR("Failed to wait on queue semaphore: %s\n", - strerror(errno)); - } - profile_end("wait-queue-semaphore"); - - /* Get the event from the queue */ - set_last_task(sb->shared->last_task[cookie], "read_queue"); - pthread_mutex_lock(&sb->shared->queue_lock); - if ( ((sb->shared->n_events==0) && (sb->shared->no_more)) - || (sb->shared->should_shutdown) ) - { - /* Queue is empty and no more are coming, - * or another process has initiated a shutdown. - * Either way, it's time to get out of here. */ - pthread_mutex_unlock(&sb->shared->queue_lock); - allDone = 1; - continue; - } - if ( sb->shared->n_events == 0 ) { - ERROR("Got the semaphore, but no events in queue!\n"); - ERROR("no_more = %i\n", sb->shared->no_more); - pthread_mutex_unlock(&sb->shared->queue_lock); - allDone = 1; - continue; - } - - line = strdup(sb->shared->queue[0]); - - len = strlen(line); - assert(len > 1); - for ( i=len-1; i>0; i-- ) { - if ( line[i] == ' ' ) { - line[i] = '\0'; - ser_str = &line[i+1]; - break; - } - } - len = strlen(line); - assert(len > 1); - for ( i=len-1; i>0; i-- ) { - if ( line[i] == ' ' ) { - line[i] = '\0'; - event_str = &line[i+1]; - break; - } - } - if ( (ser_str != NULL) && (event_str != NULL) ) { - if ( sscanf(ser_str, "%i", &ser) != 1 ) { - STATUS("Invalid serial number '%s'\n", - ser_str); + for ( v=0; v<endl; v++ ) { + if ( haystack[u+v] != needle[v] ) { ok = 0; + break; } } - if ( !ok ) { - STATUS("Invalid event string '%s'\n", - sb->shared->queue[0]); - ok = 0; - } - memcpy(sb->shared->last_ev[cookie], sb->shared->queue[0], - MAX_EV_LEN); - shuffle_events(sb->shared); - pthread_mutex_unlock(&sb->shared->queue_lock); - - if ( !ok ) continue; - - pargs.filename = strdup(line); - pargs.event = safe_strdup(event_str); - - free(line); - ok = 0; - - /* Default values */ - pargs.zmq_data = NULL; - pargs.zmq_data_size = 0; - pargs.asapo_data = NULL; - pargs.asapo_data_size = 0; - pargs.asapo_meta = NULL; - - if ( sb->zmq_params != NULL ) { - - profile_start("zmq-fetch"); - set_last_task(sb->shared->last_task[cookie], "ZMQ fetch"); - pargs.zmq_data = im_zmq_fetch(zmqstuff, - &pargs.zmq_data_size); - profile_end("zmq-fetch"); - - if ( (pargs.zmq_data != NULL) - && (pargs.zmq_data_size > 15) ) ok = 1; - - /* The filename/event, which will be 'fake' values in - * this case, still came via the event queue. More - * importantly, the event queue gave us a unique - * serial number for this image. */ - - } else if ( sb->asapo_params != NULL ) { - - char *filename; - char *event; - int finished = 0; - int asapo_message_id; - - profile_start("asapo-fetch"); - set_last_task(sb->shared->last_task[cookie], "ASAPO fetch"); - pargs.asapo_data = im_asapo_fetch(asapostuff, - &pargs.asapo_data_size, - &pargs.asapo_meta, - &filename, - &event, - &finished, - &asapo_message_id); - profile_end("asapo-fetch"); - if ( pargs.asapo_data != NULL ) { - ok = 1; - - /* ASAP::O provides a meaningful filename, which - * replaces the placeholder. */ - free(pargs.filename); - free(pargs.event); - pargs.filename = filename; - pargs.event = event; - sb->shared->end_of_stream[cookie] = 0; - - /* We will also use ASAP::O's serial number - * instead of our own. */ - ser = asapo_message_id; - } else { - if ( finished ) { - sb->shared->end_of_stream[cookie] = 1; - } - } - - } else { - ok = 1; - } - if ( ok ) { - sb->shared->time_last_start[cookie] = get_monotonic_seconds(); - profile_start("process-image"); - process_image(iargs, &pargs, st, cookie, tmpdir, ser, - sb->shared, sb->shared->last_task[cookie], - asapostuff, mille, ida); - profile_end("process-image"); - - if ( sb->asapo_params != NULL ) { - im_asapo_finalise(asapostuff, ser); - } - - } - - /* NB pargs.zmq_data, pargs.asapo_data and pargs.asapo_meta - * will be copied into the image structure, so - * that it can be queried for "header" values etc. They will - * eventually be freed by image_free() under process_image(). */ - - if ( sb->profile ) { - profile_print_and_reset(cookie); + return &haystack[u]; } - - free(pargs.filename); - free(pargs.event); } + return NULL; +} - image_data_arrays_free(ida); - crystfel_mille_free(mille); - /* These are both no-ops if argument is NULL */ - im_zmq_shutdown(zmqstuff); - im_asapo_shutdown(asapostuff); +static size_t pump_chunk(void *buf, size_t len, struct sandbox *sb) +{ + const char *txt = (char *)buf; + const char *endpos; + size_t chunk_len; - data_template_free(iargs->dtempl); - cleanup_indexing(iargs->ipriv); - cell_free(iargs->cell); - return 0; -} + endpos = str_in_str(txt, len, STREAM_CHUNK_END_MARKER"\n"); + if ( endpos == NULL ) return 0; + chunk_len = (endpos-txt)+strlen(STREAM_CHUNK_END_MARKER"\n"); + fwrite(buf, 1, chunk_len, stream_get_fh(sb->stream)); + fflush(stream_get_fh(sb->stream)); -static ssize_t lwrite(int fd, const char *a) -{ - size_t l = strlen(a); - return write(fd, a, l); + return chunk_len; } -static int pump_chunk(FILE *fh, int ofd) +static size_t pump_mille(void *buf, size_t len, struct sandbox *sb) { - int chunk_started = 0; - - do { + int n; + int ni; + size_t reclen; - char line[1024]; - char *rval; + if ( len < sizeof(int) ) return 0; + ni = *(int *)buf; + n = (ni-2)/2; - rval = fgets(line, 1024, fh); - if ( rval == NULL ) { + reclen = (2+n)*sizeof(int) + (1+n)*sizeof(float); + if ( len < reclen ) return 0; - if ( feof(fh) ) { - /* Whoops, connection lost */ - if ( chunk_started ) { - ERROR("EOF during chunk!\n"); - lwrite(ofd, "Unfinished chunk!\n"); - lwrite(ofd, STREAM_CHUNK_END_MARKER"\n"); - } /* else normal end of output */ - return 1; - } + fwrite(buf, 1, reclen, sb->mille_fh); + fflush(sb->mille_fh); - ERROR("fgets() failed: %s\n", strerror(errno)); - if ( errno != EINTR ) return 1; + return reclen; +} - } - if ( strcmp(line, "FLUSH\n") == 0 ) break; - lwrite(ofd, line); +static PipeList *pipe_list_new() +{ + PipeList *pd = malloc(sizeof(PipeList)); + if ( pd == NULL ) return NULL; + pd->fds = NULL; + pd->buffers = NULL; + pd->buffer_len = NULL; + pd->buffer_pos = NULL; + pd->n_read = 0; + return pd; +} - if ( strcmp(line, STREAM_CHUNK_START_MARKER"\n") == 0 ) { - chunk_started = 1; - } - if ( strcmp(line, STREAM_CHUNK_END_MARKER"\n") == 0 ) break; - } while ( 1 ); - return 0; +static void pipe_list_destroy(PipeList *pd) +{ + int i; + for ( i=0; i<pd->n_read; i++ ) { + free(pd->buffers[i]); + } + free(pd->buffers); + free(pd->buffer_len); + free(pd->buffer_pos); + free(pd->fds); + free(pd); } /* Add an fd to the list of pipes to be read from */ -static void add_pipe(struct sandbox *sb, int fd) +static void add_pipe(PipeList *pd, int fd) { int *fds_new; - FILE **fhs_new; + void **buffers_new; + size_t *buflens_new; + size_t *bufposs_new; int slot; - fds_new = realloc(sb->fds, (sb->n_read+1)*sizeof(int)); - if ( fds_new == NULL ) { + fds_new = realloc(pd->fds, (pd->n_read+1)*sizeof(int)); + buffers_new = realloc(pd->buffers, (pd->n_read+1)*sizeof(void *)); + buflens_new = realloc(pd->buffer_len, (pd->n_read+1)*sizeof(size_t)); + bufposs_new = realloc(pd->buffer_pos, (pd->n_read+1)*sizeof(size_t)); + if ( (fds_new == NULL) || (buffers_new == NULL) + || (buflens_new == NULL) || (bufposs_new == NULL) ) + { ERROR("Failed to allocate memory for new pipe.\n"); return; } - fhs_new = realloc(sb->fhs, (sb->n_read+1)*sizeof(FILE *)); - if ( fhs_new == NULL ) { - ERROR("Failed to allocate memory for new FH.\n"); - free(fds_new); - return; - } + pd->fds = fds_new; + pd->buffers = buffers_new; + pd->buffer_len = buflens_new; + pd->buffer_pos = bufposs_new; - sb->fds = fds_new; - sb->fhs = fhs_new; - slot = sb->n_read; + slot = pd->n_read; + pd->fds[slot] = fd; + pd->buffer_len[slot] = 64*1024; /* Initial buffer size */ + pd->buffers[slot] = malloc(pd->buffer_len[slot]); + if ( pd->buffers[slot] == NULL ) return; + pd->buffer_pos[slot] = 0; - sb->fds[slot] = fd; - - sb->fhs[slot] = fdopen(fd, "r"); - if ( sb->fhs[slot] == NULL ) { - ERROR("Couldn't fdopen() stream!\n"); - return; - } - - sb->n_read++; + pd->n_read++; } -static void remove_pipe(struct sandbox *sb, int d) +static void remove_pipe(PipeList *pd, int d) { int i; - fclose(sb->fhs[d]); + close(pd->fds[d]); + free(pd->buffers[d]); - for ( i=d; i<sb->n_read; i++ ) { - if ( i < sb->n_read-1 ) { - sb->fds[i] = sb->fds[i+1]; - sb->fhs[i] = sb->fhs[i+1]; + for ( i=d; i<pd->n_read; i++ ) { + if ( i < pd->n_read-1 ) { + pd->fds[i] = pd->fds[i+1]; + pd->buffers[i] = pd->buffers[i+1]; + pd->buffer_len[i] = pd->buffer_len[i+1]; + pd->buffer_pos[i] = pd->buffer_pos[i+1]; } /* else don't bother */ } - sb->n_read--; + pd->n_read--; /* We don't bother shrinking the arrays */ } -static void try_read(struct sandbox *sb) +static int find_marked(PipeList *pd) +{ + int i; + for ( i=0; i<pd->n_read; i++ ) { + if ( pd->buffer_len[i] == 0 ) return i; + } + return -1; +} + + +static void check_pipes(PipeList *pd, size_t(*pump)(void *, size_t len, struct sandbox *), + struct sandbox *sb) { int r, i; struct timeval tv; fd_set fds; int fdmax; - const int ofd = stream_get_fd(sb->stream); tv.tv_sec = 0; tv.tv_usec = 500000; FD_ZERO(&fds); fdmax = 0; - for ( i=0; i<sb->n_read; i++ ) { + for ( i=0; i<pd->n_read; i++ ) { int fd; - fd = sb->fds[i]; + fd = pd->fds[i]; FD_SET(fd, &fds); if ( fd > fdmax ) fdmax = fd; @@ -702,33 +530,60 @@ static void try_read(struct sandbox *sb) return; } - for ( i=0; i<sb->n_read; i++ ) { + for ( i=0; i<pd->n_read; i++ ) { + + size_t r; - if ( !FD_ISSET(sb->fds[i], &fds) ) { + if ( !FD_ISSET(pd->fds[i], &fds) ) { continue; } + if ( pd->buffer_len[i] == pd->buffer_pos[i] ) { + const size_t buffer_increment = 64*1024; + void *buf_new = realloc(pd->buffers[i], + pd->buffer_len[i]+buffer_increment); + if ( buf_new == NULL ) { + ERROR("Failed to grow buffer\n"); + continue; + } + pd->buffers[i] = buf_new; + pd->buffer_len[i] += buffer_increment; + } + /* If the chunk cannot be read, assume the connection * is broken and that the process will die soon. */ - if ( pump_chunk(sb->fhs[i], ofd) ) { - remove_pipe(sb, i); + r = read(pd->fds[i], pd->buffers[i]+pd->buffer_pos[i], + pd->buffer_len[i]-pd->buffer_pos[i]); + + if ( r == 0 ) { + /* Mark for deletion */ + pd->buffer_len[i] = 0; + } else { + size_t h; + pd->buffer_pos[i] += r; + h = pump(pd->buffers[i], pd->buffer_pos[i], sb); + assert(h <= pd->buffer_pos[i]); + assert(h >= 0); + if ( h > 0 ) { + memmove(pd->buffers[i], pd->buffers[i]+h, + pd->buffer_pos[i]-h); + pd->buffer_pos[i] -= h; + } } + } + int deleteme = find_marked(pd); + while ( deleteme != -1 ) { + remove_pipe(pd, deleteme); + deleteme = find_marked(pd); } } -static void pin_to_cpu(int slot) +static void try_read(struct sandbox *sb) { - #ifdef HAVE_SCHED_SETAFFINITY - cpu_set_t c; - - CPU_ZERO(&c); - CPU_SET(slot, &c); - if ( sched_setaffinity(0, sizeof(cpu_set_t), &c) ) { - fprintf(stderr, "Failed to set CPU affinity for %i\n", slot); - } - #endif + check_pipes(sb->st_from_workers, pump_chunk, sb); + check_pipes(sb->mille_from_workers, pump_mille, sb); } @@ -736,123 +591,135 @@ static void start_worker_process(struct sandbox *sb, int slot) { pid_t p; int stream_pipe[2]; + int mille_pipe[2]; + char **nargv; + int nargc; + int i; + char *tmpdir_copy; + char *methods_copy = NULL; + char *worker_id; + char *fd_stream; + char *fd_mille; + char buf[1024]; + const char *indexamajig = NULL; + size_t len; if ( pipe(stream_pipe) == - 1 ) { ERROR("pipe() failed!\n"); return; } + if ( pipe(mille_pipe) == - 1 ) { + ERROR("pipe() failed!\n"); + return; + } + pthread_mutex_lock(&sb->shared->queue_lock); - sb->shared->pings[slot] = 0; sb->shared->end_of_stream[slot] = 0; - sb->last_ping[slot] = 0; - sb->shared->time_last_start[slot] = get_monotonic_seconds(); - sb->shared->warned_long_running[slot] = 0; pthread_mutex_unlock(&sb->shared->queue_lock); - p = fork(); - if ( p == -1 ) { - ERROR("fork() failed!\n"); - return; - } - - if ( p == 0 ) { + pthread_mutex_lock(&sb->shared->debug_lock); + sb->shared->pings[slot] = 0; + sb->last_ping[slot] = 0; + sb->shared->time_last_start[slot] = get_monotonic_seconds(); + pthread_mutex_unlock(&sb->shared->debug_lock); - Stream *st; - struct sigaction sa; - int r; - char *tmp; - struct stat s; - size_t ll; - int i; + sb->warned_long_running[slot] = 0; - if ( sb->cpu_pin ) pin_to_cpu(slot); - - /* First, disconnect the signal handlers */ - sa.sa_flags = 0; - sigemptyset(&sa.sa_mask); - sa.sa_handler = SIG_DFL; - r = sigaction(SIGCHLD, &sa, NULL); - if ( r == -1 ) { - ERROR("Failed to set signal handler!\n"); - exit(1); - } - r = sigaction(SIGINT, &sa, NULL); - if ( r == -1 ) { - ERROR("Failed to set signal handler!\n"); - exit(1); - } - r = sigaction(SIGQUIT, &sa, NULL); - if ( r == -1 ) { - ERROR("Failed to set signal handler!\n"); - exit(1); - } - - sa.sa_handler = SIG_IGN; - r = sigaction(SIGUSR1, &sa, NULL); - if ( r == -1 ) { - ERROR("Failed to set signal handler!\n"); - exit(1); - } - - ll = 64 + strlen(sb->tmpdir); - tmp = malloc(ll); - if ( tmp == NULL ) { - ERROR("Failed to allocate temporary dir\n"); - exit(1); - } + /* Set up nargv including "new" args */ + nargc = 0; + nargv = malloc((sb->argc+16)*sizeof(char *)); + if ( nargv == NULL ) return; + for ( i=0; i<sb->argc; i++ ) { + nargv[nargc++] = sb->argv[i]; + } - snprintf(tmp, ll, "%s/worker.%i", sb->tmpdir, slot); + nargv[nargc++] = "--shm-name"; + nargv[nargc++] = sb->shm_name; + + nargv[nargc++] = "--queue-sem"; + nargv[nargc++] = sb->sem_name; + + nargv[nargc++] = "--worker-tmpdir"; + tmpdir_copy = strdup(sb->tmpdir); + if ( tmpdir_copy == NULL ) return; + nargv[nargc++] = tmpdir_copy; + + nargv[nargc++] = "--worker-id"; + worker_id = malloc(64); + if ( worker_id == NULL ) return; + snprintf(worker_id, 64, "%i", slot); + nargv[nargc++] = worker_id; + + nargv[nargc++] = "--fd-stream"; + fd_stream = malloc(64); + snprintf(fd_stream, 64, "%i", stream_pipe[1]); + nargv[nargc++] = fd_stream; + + nargv[nargc++] = "--fd-mille"; + fd_mille = malloc(64); + snprintf(fd_mille, 64, "%i", mille_pipe[1]); + nargv[nargc++] = fd_mille; + + if ( sb->probed_methods != NULL ) { + methods_copy = strdup(sb->probed_methods); + nargv[nargc++] = "--indexing"; + nargv[nargc++] = methods_copy; + } + nargv[nargc++] = NULL; - if ( stat(tmp, &s) == -1 ) { - if ( errno != ENOENT ) { - ERROR("Failed to stat temporary folder.\n"); - exit(1); - } + len = readlink("/proc/self/exe", buf, 1024); - r = mkdir(tmp, S_IRWXU); - if ( r ) { - ERROR("Failed to create temporary folder: %s\n", - strerror(errno)); - exit(1); - } + if ( (len == -1) || (len >= 1023) ) { + ERROR("readlink() failed: %s\n", strerror(errno)); + } else { + buf[len] = '\0'; + if ( strstr(buf, "indexamajig") == NULL ) { + ERROR("Didn't recognise /proc/self/exe (%s)\n", buf); + } else { + indexamajig = buf; } - - /* Free resources which will not be needed by worker */ - free(sb->pids); - for ( i=0; i<sb->n_read; i++ ) { - fclose(sb->fhs[i]); + } + if ( indexamajig == NULL ) { + if ( strstr(sb->argv[0], "indexamajig") == NULL ) { + ERROR("Didn't recognise argv[0] (%s)\n", sb->argv[0]); + } else { + indexamajig = sb->argv[0]; } - free(sb->fhs); - free(sb->fds); - free(sb->running); - /* Not freed because it's not worth passing them down just for - * this purpose: event list file handle, - * main output stream handle - * original temp dir name (without indexamajig.XX) - * prefix - */ - - st = stream_open_fd_for_write(stream_pipe[1], sb->iargs->dtempl); - r = run_work(sb->iargs, st, slot, tmp, sb); - stream_close(st); - - free(tmp); - - munmap(sb->shared, sizeof(struct sb_shm)); - - free(sb); + } + if ( indexamajig == NULL ) { + ERROR("Falling back on shell search path.\n"); + indexamajig = "indexamajig"; + } - exit(r); + p = fork(); + if ( p == -1 ) { + ERROR("fork() failed!\n"); + return; + } + if ( p == 0 ) { + execvp(indexamajig, nargv); + ERROR("Failed to exec!\n"); + return; } + free(tmpdir_copy); + free(methods_copy); + free(worker_id); + free(fd_stream); + free(fd_mille); + free(nargv); + /* Parent process gets the 'write' end of the filename pipe * and the 'read' end of the result pipe. */ sb->pids[slot] = p; sb->running[slot] = 1; + pthread_mutex_lock(&sb->shared->debug_lock); stamp_response(sb, slot); - add_pipe(sb, stream_pipe[0]); + pthread_mutex_unlock(&sb->shared->debug_lock); + add_pipe(sb->st_from_workers, stream_pipe[0]); + add_pipe(sb->mille_from_workers, mille_pipe[0]); close(stream_pipe[1]); } @@ -889,6 +756,15 @@ static void handle_zombie(struct sandbox *sb, int respawn) sb->running[i] = 0; if ( WIFEXITED(status) ) { + if ( WEXITSTATUS(status) != 0 ) { + STATUS("Worker %i returned error code %i\n", + i, WEXITSTATUS(status)); + STATUS("Shutting down.\n"); + /* Error status from worker */ + pthread_mutex_lock(&sb->shared->totals_lock); + sb->shared->should_shutdown = 1; + pthread_mutex_unlock(&sb->shared->totals_lock); + } continue; } @@ -915,37 +791,75 @@ static void handle_zombie(struct sandbox *sb, int respawn) static int setup_shm(struct sandbox *sb) { pthread_mutexattr_t attr; + char tmp[128]; + int shm_fd; - sb->shared = mmap(NULL, sizeof(struct sb_shm), PROT_READ | PROT_WRITE, - MAP_SHARED | MAP_ANON, -1, 0); + snprintf(tmp, 127, "/indexamajig.shm.%i", getpid()); + shm_fd = shm_open(tmp, O_CREAT | O_EXCL | O_RDWR, 0600); + if ( shm_fd == -1 ) { + ERROR("SHM setup failed: %s\n", strerror(errno)); + return 1; + } + sb->shm_name = strdup(tmp); + if ( ftruncate(shm_fd, sizeof(struct sb_shm)) == -1 ) { + ERROR("SHM setup failed: %s\n", strerror(errno)); + free(sb->shm_name); + return 1; + } + + sb->shared = mmap(NULL, sizeof(struct sb_shm), PROT_READ | PROT_WRITE, + MAP_SHARED, shm_fd, 0); if ( sb->shared == MAP_FAILED ) { ERROR("SHM setup failed: %s\n", strerror(errno)); + free(sb->shm_name); return 1; } if ( pthread_mutexattr_init(&attr) ) { ERROR("Failed to initialise mutex attr.\n"); + free(sb->shm_name); return 1; } if ( pthread_mutexattr_setpshared(&attr, PTHREAD_PROCESS_SHARED) ) { ERROR("Failed to set process shared attribute.\n"); + pthread_mutexattr_destroy(&attr); + free(sb->shm_name); return 1; } if ( pthread_mutex_init(&sb->shared->term_lock, &attr) ) { ERROR("Terminal lock setup failed.\n"); + pthread_mutexattr_destroy(&attr); + free(sb->shm_name); return 1; } if ( pthread_mutex_init(&sb->shared->queue_lock, &attr) ) { ERROR("Queue lock setup failed.\n"); + pthread_mutexattr_destroy(&attr); + pthread_mutex_destroy(&sb->shared->term_lock); + free(sb->shm_name); + return 1; + } + + if ( pthread_mutex_init(&sb->shared->debug_lock, &attr) ) { + ERROR("Queue lock setup failed.\n"); + pthread_mutexattr_destroy(&attr); + pthread_mutex_destroy(&sb->shared->term_lock); + pthread_mutex_destroy(&sb->shared->queue_lock); + free(sb->shm_name); return 1; } if ( pthread_mutex_init(&sb->shared->totals_lock, &attr) ) { ERROR("Totals lock setup failed.\n"); + pthread_mutexattr_destroy(&attr); + pthread_mutex_destroy(&sb->shared->term_lock); + pthread_mutex_destroy(&sb->shared->queue_lock); + pthread_mutex_destroy(&sb->shared->debug_lock); + free(sb->shm_name); return 1; } @@ -1013,8 +927,7 @@ static void sigusr1_handler(int sig, siginfo_t *si, void *uc_v) } -static void check_signals(struct sandbox *sb, const char *semname_q, - int respawn) +static void check_signals(struct sandbox *sb, int respawn) { if ( at_zombies ) { at_zombies = 0; @@ -1022,7 +935,8 @@ static void check_signals(struct sandbox *sb, const char *semname_q, } if ( at_interrupt ) { - sem_unlink(semname_q); + sem_unlink(sb->sem_name); + shm_unlink(sb->shm_name); exit(0); } @@ -1040,8 +954,8 @@ static void try_status(struct sandbox *sb, int final) { int r; int n_proc_this; - double tNow; - double time_this; + time_t tNow; + time_t time_this; const char *finalstr; char persec[64]; @@ -1049,7 +963,9 @@ static void try_status(struct sandbox *sb, int final) time_this = tNow - sb->t_last_stats; if ( !final && (time_this < 5) ) return; + pthread_mutex_lock(&sb->shared->totals_lock); n_proc_this = sb->shared->n_processed - sb->n_processed_last_stats; + pthread_mutex_unlock(&sb->shared->totals_lock); r = pthread_mutex_trylock(&sb->shared->term_lock); if ( r ) return; /* No lock -> don't bother */ @@ -1199,7 +1115,8 @@ int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, struct im_zmq_params *zmq_params, struct im_asapo_params *asapo_params, int timeout, int profile, int cpu_pin, - int no_data_timeout) + int no_data_timeout, int argc, char *argv[], + const char *probed_methods, FILE *mille_fh) { int i; struct sandbox *sb; @@ -1208,7 +1125,7 @@ int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int r; int allDone = 0; struct get_pattern_ctx gpctx; - double t_last_data; + time_t t_last_data; if ( n_proc > MAX_NUM_WORKERS ) { ERROR("Number of workers (%i) is too large. Using %i\n", @@ -1248,7 +1165,10 @@ int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, sb->tmpdir = tmpdir; sb->profile = profile; sb->timeout = timeout; - sb->cpu_pin = cpu_pin; + sb->argc = argc; + sb->argv = argv; + sb->probed_methods = probed_methods; + sb->mille_fh = mille_fh; if ( zmq_params->addr != NULL ) { sb->zmq_params = zmq_params; @@ -1268,8 +1188,8 @@ int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, return 0; } - sb->fds = NULL; - sb->fhs = NULL; + sb->st_from_workers = pipe_list_new(); + sb->mille_from_workers = pipe_list_new(); sb->stream = stream; gpctx.fh = fh; @@ -1300,6 +1220,7 @@ int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, ERROR("Failed to create semaphore: %s\n", strerror(errno)); return 0; } + sb->sem_name = strdup(semname_q); sb->pids = calloc(n_proc, sizeof(pid_t)); sb->running = calloc(n_proc, sizeof(int)); @@ -1364,7 +1285,7 @@ int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, try_read(sb); /* Check for interrupt or zombies */ - check_signals(sb, semname_q, 1); + check_signals(sb, 1); /* Check for hung workers */ check_hung_workers(sb); @@ -1381,6 +1302,7 @@ int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, /* Begin exit criterion checking */ pthread_mutex_lock(&sb->shared->queue_lock); + pthread_mutex_lock(&sb->shared->totals_lock); /* Case 1: Queue empty and no more coming? */ if ( sb->shared->no_more && (sb->shared->n_events == 0) ) allDone = 1; @@ -1401,6 +1323,7 @@ int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, t_last_data = get_monotonic_seconds(); } + pthread_mutex_unlock(&sb->shared->totals_lock); pthread_mutex_unlock(&sb->shared->queue_lock); /* End exit criterion checking */ @@ -1422,7 +1345,7 @@ int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, for ( i=0; i<n_proc; i++ ) { while ( any_running(sb) ) { try_read(sb); - check_signals(sb, semname_q, 0); + check_signals(sb, 0); check_hung_workers(sb); try_status(sb, 0); } @@ -1433,11 +1356,8 @@ int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, sem_unlink(semname_q); sem_close(sb->queue_sem); - for ( i=0; i<sb->n_read; i++ ) { - fclose(sb->fhs[i]); - } - free(sb->fhs); - free(sb->fds); + pipe_list_destroy(sb->st_from_workers); + pipe_list_destroy(sb->mille_from_workers); free(sb->running); free(sb->last_response); free(sb->pids); @@ -1448,7 +1368,10 @@ int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, delete_temporary_folder(sb->tmpdir, n_proc); + shm_unlink(sb->shm_name); munmap(sb->shared, sizeof(struct sb_shm)); + free(sb->shm_name); + free(sb->sem_name); free(sb); return r; diff --git a/src/im-sandbox.h b/src/im-sandbox.h index ba3a6d3d..87788668 100644 --- a/src/im-sandbox.h +++ b/src/im-sandbox.h @@ -51,8 +51,7 @@ struct sb_shm; /* Maximum length of an event ID including serial number */ #define MAX_EV_LEN (1024) -/* Maximum length of a task ID, e.g. indexing:xgandalf. - * NB If changing this, also update the value in index.c */ +/* Maximum length of a task ID, e.g. indexing:xgandalf. */ #define MAX_TASK_LEN (32) /* Maximum number of workers */ @@ -66,12 +65,13 @@ struct sb_shm int n_events; char queue[QUEUE_SIZE][MAX_EV_LEN]; int no_more; + int end_of_stream[MAX_NUM_WORKERS]; + + pthread_mutex_t debug_lock; char last_ev[MAX_NUM_WORKERS][MAX_EV_LEN]; char last_task[MAX_NUM_WORKERS][MAX_TASK_LEN]; int pings[MAX_NUM_WORKERS]; - int end_of_stream[MAX_NUM_WORKERS]; time_t time_last_start[MAX_NUM_WORKERS]; - int warned_long_running[MAX_NUM_WORKERS]; pthread_mutex_t totals_lock; int n_processed; @@ -83,7 +83,7 @@ struct sb_shm extern char *create_tempdir(const char *temp_location); -extern void set_last_task(char *lt, const char *task); +extern time_t get_monotonic_seconds(void); extern int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, Stream *stream, @@ -91,6 +91,7 @@ extern int create_sandbox(struct index_args *iargs, int n_proc, char *prefix, struct im_zmq_params *zmq_params, struct im_asapo_params *asapo_params, int timeout, int profile, int cpu_pin, - int no_data_timeout); + int no_data_timeout, int argc, char *argv[], + const char *probed_methods, FILE *mille_fh); #endif /* IM_SANDBOX_H */ diff --git a/src/im-zmq.c b/src/im-zmq.c index ece80662..b7a2ed40 100644 --- a/src/im-zmq.c +++ b/src/im-zmq.c @@ -57,7 +57,7 @@ struct im_zmq }; -struct im_zmq *im_zmq_connect(struct im_zmq_params *params) +struct im_zmq *im_zmq_connect(const struct im_zmq_params *params) { struct im_zmq *z; diff --git a/src/im-zmq.h b/src/im-zmq.h index cfb98ff4..fb99f4a6 100644 --- a/src/im-zmq.h +++ b/src/im-zmq.h @@ -46,13 +46,13 @@ struct im_zmq_params #if defined(HAVE_ZMQ) -extern struct im_zmq *im_zmq_connect(struct im_zmq_params *params); +extern struct im_zmq *im_zmq_connect(const struct im_zmq_params *params); extern void im_zmq_shutdown(struct im_zmq *z); extern void *im_zmq_fetch(struct im_zmq *z, size_t *pdata_size); #else /* defined(HAVE_ZMQ) */ -static UNUSED struct im_zmq *im_zmq_connect(struct im_zmq_params *params) { return NULL; } +static UNUSED struct im_zmq *im_zmq_connect(const struct im_zmq_params *params) { return NULL; } static UNUSED void im_zmq_shutdown(struct im_zmq *z) { } static UNUSED void *im_zmq_fetch(struct im_zmq *z, size_t *psize) { *psize = 0; return NULL; } diff --git a/src/indexamajig.c b/src/indexamajig.c index 29e60d88..c0cb62e7 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -38,16 +38,25 @@ #include <config.h> #endif -#include <stdarg.h> +#ifdef HAVE_SCHED_SETAFFINITY +#define _GNU_SOURCE +#include <sys/sysinfo.h> +#include <sched.h> +#endif + #include <stdlib.h> #include <stdio.h> #include <string.h> #include <unistd.h> -#include <argp.h> #include <gsl/gsl_errno.h> -#include <sys/types.h> #include <sys/stat.h> #include <fcntl.h> +#include <assert.h> +#include <sys/mman.h> + +#ifdef HAVE_FFTW +#include <fftw3.h> +#endif #include <utils.h> #include <index.h> @@ -65,48 +74,12 @@ #include <peakfinder8.h> #include "im-sandbox.h" +#include "im-argparse.h" #include "im-zmq.h" #include "im-asapo.h" #include "version.h" #include "json-utils.h" - - -struct indexamajig_arguments -{ - struct index_args iargs; /* These are the options that will be - * given to process_image */ - char *filename; - char *geom_filename; - char *outfile; - char *prefix; - int check_prefix; - int n_proc; - char *cellfile; - char *indm_str; - int basename; - struct im_zmq_params zmq_params; - struct im_asapo_params asapo_params; - int serial_start; - char *temp_location; - int if_refine; - int if_checkcell; - int if_peaks; - int if_multi; - int if_retry; - int profile; /* Whether to do wall-clock time profiling */ - int no_data_timeout; - char **copy_headers; - int n_copy_headers; - char *harvest_file; - int cpu_pin; - - struct taketwo_options **taketwo_opts_ptr; - struct felix_options **felix_opts_ptr; - struct xgandalf_options **xgandalf_opts_ptr; - struct pinkindexer_options **pinkindexer_opts_ptr; - struct fromfile_options **fromfile_opts_ptr; - struct asdf_options **asdf_opts_ptr; -}; +#include "profile.h" static double nan_if_neg(double n) @@ -275,620 +248,421 @@ static void write_harvest_file(struct index_args *args, } -static void show_version(FILE *fh, struct argp_state *state) +static void shuffle_events(struct sb_shm *sb_shared) { - printf("CrystFEL: %s\n", crystfel_version_string()); - printf("%s\n", crystfel_licence_string()); + int i; + + for ( i=1; i<sb_shared->n_events; i++ ) { + memcpy(sb_shared->queue[i-1], sb_shared->queue[i], MAX_EV_LEN); + } + sb_shared->n_events--; } -static void add_copy_header(struct indexamajig_arguments *args, - const char *header) +static void pin_to_cpu(int slot) { - char **new_copy_headers = realloc(args->copy_headers, - (args->n_copy_headers+1)*sizeof(char *)); - if ( new_copy_headers == NULL ) { - ERROR("Failed to add copy header '%s'\n", header); - return; - } + #ifdef HAVE_SCHED_SETAFFINITY + cpu_set_t c; - args->copy_headers = new_copy_headers; - args->copy_headers[args->n_copy_headers++] = strdup(header); + CPU_ZERO(&c); + CPU_SET(slot, &c); + if ( sched_setaffinity(0, sizeof(cpu_set_t), &c) ) { + fprintf(stderr, "Failed to set CPU affinity for %i\n", slot); + } + #endif } -static DataSourceType parse_data_format(const char *str) +struct sb_shm *shared; +int _worker; + +static void set_last_task_sandbox(const char *task) { - if ( strcmp(str, "hdf5") == 0 ) return DATA_SOURCE_TYPE_HDF5; - if ( strcmp(str, "msgpack") == 0 ) return DATA_SOURCE_TYPE_MSGPACK; - if ( strcmp(str, "seedee") == 0 ) return DATA_SOURCE_TYPE_SEEDEE; - /* CBF and CBFGZ should be added here once image-cbf.c supports - * in-memory access */ - return DATA_SOURCE_TYPE_UNKNOWN; + assert(strlen(task) < MAX_TASK_LEN-1); + pthread_mutex_lock(&shared->debug_lock); + strcpy(shared->last_task[_worker], task); + pthread_mutex_unlock(&shared->debug_lock); } - -static error_t parse_arg(int key, char *arg, struct argp_state *state) +static void notify_alive_sandbox() { - float tmp; - int r; - struct indexamajig_arguments *args = state->input; - - switch ( key ) { - - case ARGP_KEY_INIT : - state->child_inputs[0] = args->taketwo_opts_ptr; - state->child_inputs[1] = args->felix_opts_ptr; - state->child_inputs[2] = args->xgandalf_opts_ptr; - state->child_inputs[3] = args->pinkindexer_opts_ptr; - state->child_inputs[4] = args->fromfile_opts_ptr; - state->child_inputs[5] = args->asdf_opts_ptr; - break; - - case 'h' : - argp_state_help(state, stdout, ARGP_HELP_STD_HELP); - break; /* argp_state_help doesn't return */ - - case 'v' : - show_version(stdout, state); - exit(0); - - case 'i' : - args->filename = strdup(arg); - break; - - case 'o' : - args->outfile = strdup(arg); - break; - - case 'x' : - free(args->prefix); - args->prefix = strdup(arg); - break; - - case 'j' : - args->n_proc = atoi(arg); - break; - - case 'g' : - args->geom_filename = arg; - break; - - case 201 : - args->basename = 1; - break; - - case 202 : - args->check_prefix = 0; - break; - - case 203 : - if ( sscanf(arg, "%f", &tmp) != 1 ) { - ERROR("Invalid value for --highres\n"); - return EINVAL; - } - args->iargs.highres = 1.0 / (tmp/1e10); /* A -> m^-1 */ - break; - - case 204 : - args->profile = 1; - break; - - case 205 : - args->temp_location = strdup(arg); - break; + pthread_mutex_lock(&shared->debug_lock); + shared->pings[_worker]++; + pthread_mutex_unlock(&shared->debug_lock); +} - case 206 : - if (sscanf(arg, "%d", &args->iargs.wait_for_file) != 1) - { - ERROR("Invalid value for --wait-for-file\n"); - return EINVAL; - } - break; - case 207 : - args->zmq_params.addr = strdup(arg); - break; +static int run_work(struct indexamajig_arguments *args) +{ + int allDone = 0; + struct im_zmq *zmqstuff = NULL; + struct im_asapo *asapostuff = NULL; + Mille *mille; + ImageDataArrays *ida; + size_t ll; + char *tmp; + struct stat s; + Stream *st; + int shm_fd; + sem_t *queue_sem; + IndexingFlags flags = 0; + struct pf8_private_data *pf8_data = NULL; - case 208 : - args->iargs.no_image_data = 1; - break; + if ( args->cpu_pin ) pin_to_cpu(args->worker_id); + _worker = args->worker_id; - case 209 : - ERROR("--spectrum-filename is no longer used.\n"); + ll = 64 + strlen(args->worker_tmpdir); + tmp = malloc(ll); + if ( tmp == NULL ) { + ERROR("Failed to allocate temporary dir\n"); return 1; + } - case 210 : - args->iargs.no_mask_data = 1; - break; - - case 211 : - if ( args->zmq_params.n_subscriptions == 256 ) { - ERROR("Too many ZMQ subscriptions.\n"); - return 1; - } - args->zmq_params.subscriptions[args->zmq_params.n_subscriptions++] = strdup(arg); - break; - - case 212 : - args->zmq_params.request = strdup(arg); - break; - - case 213 : - args->asapo_params.endpoint = strdup(arg); - break; - - case 214 : - args->asapo_params.token = strdup(arg); - break; - - case 215 : - args->asapo_params.beamtime = strdup(arg); - break; - - case 217 : - args->asapo_params.group_id = strdup(arg); - break; - - case 218 : - args->asapo_params.source = strdup(arg); - break; - - case 219 : - args->iargs.data_format = parse_data_format(arg); - if ( args->iargs.data_format == DATA_SOURCE_TYPE_UNKNOWN ) { - ERROR("Unrecognised data format '%s'\n", arg); - return EINVAL; - } - break; - - case 220 : - args->asapo_params.stream = strdup(arg); - break; - - case 221 : - args->asapo_params.wait_for_stream = 1; - break; + snprintf(tmp, ll, "%s/worker.%i", args->worker_tmpdir, args->worker_id); - case 222 : - args->asapo_params.write_output_stream = 1; - break; + if ( stat(tmp, &s) == -1 ) { - case 223 : - args->cpu_pin = 1; - break; + int r; - case 224 : - if (sscanf(arg, "%d", &args->no_data_timeout) != 1) - { - ERROR("Invalid value for --no-data-timeout\n"); - return EINVAL; + if ( errno != ENOENT ) { + ERROR("Failed to stat temporary folder.\n"); + exit(1); } - break; - case 225 : - if (sscanf(arg, "%d", &args->asapo_params.consumer_timeout_ms) != 1) - { - ERROR("Invalid value for --asapo-consumer-timeout\n"); - return EINVAL; + r = mkdir(tmp, S_IRWXU); + if ( r ) { + ERROR("Failed to create temporary folder: %s\n", + strerror(errno)); + exit(1); } - break; - - case 226 : - args->asapo_params.use_ack = 1; - break; - - /* ---------- Peak search ---------- */ + } - case 't' : - args->iargs.peak_search.threshold = strtof(arg, NULL); - break; + st = stream_open_fd_for_write(args->fd_stream, args->iargs.dtempl); - case 301 : - args->iargs.peak_search.method = parse_peaksearch(arg); - if ( args->iargs.peak_search.method == PEAK_ERROR ) { - ERROR("Unrecognised peak detection method '%s'\n", arg); - return EINVAL; - } - break; - - case 302 : - r = sscanf(arg, "%f,%f,%f", &args->iargs.peak_search.pk_inn, - &args->iargs.peak_search.pk_mid, &args->iargs.peak_search.pk_out); - if ( (r != 3) || (args->iargs.peak_search.pk_inn < 0) ) { - ERROR("Invalid parameters for '--peak-radius'\n"); - return EINVAL; - } - break; - - case 303 : - if (sscanf(arg, "%d", &args->iargs.min_peaks) != 1) - { - ERROR("Invalid value for --min-peaks\n"); - return EINVAL; - } - break; - - case 304 : - ERROR("The option --hdf5-peaks is no longer used.\n"); - ERROR("Set the peak path in the geometry file.\n"); - break; + if ( args->profile ) { + profile_init(); + } - case 305 : - if (sscanf(arg, "%d", &args->iargs.peak_search.median_filter) != 1) - { - ERROR("Invalid value for --median-filter\n"); - return EINVAL; + /* Load unit cell (if given) */ + if ( args->cellfile != NULL ) { + args->iargs.cell = load_cell_from_file(args->cellfile); + if ( args->iargs.cell == NULL ) { + ERROR("Couldn't read unit cell (from %s)\n", args->cellfile); + return 1; } - break; + } else { + args->iargs.cell = NULL; + } - case 306 : - args->iargs.peak_search.noisefilter = 1; - break; + if ( args->if_checkcell ) { + flags |= INDEXING_CHECK_CELL; + } + if ( args->if_refine ) { + flags |= INDEXING_REFINE; + } + if ( args->if_peaks ) { + flags |= INDEXING_CHECK_PEAKS; + } + if ( args->if_multi ) { + flags |= INDEXING_MULTI; + } + if ( args->if_retry ) { + flags |= INDEXING_RETRY; + } - case 307 : - if (sscanf(arg, "%f", &args->iargs.peak_search.min_sq_gradient) != 1) - { - ERROR("Invalid value for --min-squared-gradient\n"); - return EINVAL; - } - break; + args->iargs.ipriv = setup_indexing(args->indm_str, + args->iargs.cell, + args->iargs.tols, + flags, + args->iargs.wavelength_estimate, + args->iargs.clen_estimate, + args->iargs.n_threads, + *args->taketwo_opts_ptr, + *args->xgandalf_opts_ptr, + *args->pinkindexer_opts_ptr, + *args->felix_opts_ptr, + *args->fromfile_opts_ptr, + *args->asdf_opts_ptr); - case 308 : - if (sscanf(arg, "%f", &args->iargs.peak_search.min_snr) != 1) - { - ERROR("Invalid value for --min-snr\n"); - return EINVAL; - } - break; + if ( args->iargs.ipriv == NULL ) { + ERROR("Failed to set up indexing system\n"); + return 1; + } - case 309 : - if (sscanf(arg, "%d", &args->iargs.peak_search.min_pix_count) != 1) - { - ERROR("Invalid value for --min-pix-count\n"); - return EINVAL; - } - break; + if ( _worker == 0 ) { + print_indexing_info(args->iargs.ipriv); + } - case 310 : - if (sscanf(arg, "%d", &args->iargs.peak_search.max_pix_count) != 1) - { - ERROR("Invalid value for --max-pix-count\n"); - return EINVAL; + if ( args->iargs.peak_search.method == PEAK_PEAKFINDER8 ) { + struct detgeom *dg; + dg = data_template_get_2d_detgeom_if_possible(args->iargs.dtempl); + if ( dg == NULL ) { + ERROR("WARNING: Detector geometry is not static. " + "Peak search will be slower than optimal.\n"); } - break; + pf8_data = prepare_peakfinder8(dg, args->iargs.peak_search.peakfinder8_fast); + args->iargs.pf_private = pf8_data; + detgeom_free(dg); + } - case 311 : - if (sscanf(arg, "%d", &args->iargs.peak_search.local_bg_radius) != 1) - { - ERROR("Invalid value for --local-bg-radius\n"); - return EINVAL; - } - break; + /* Set up SHM */ + shm_fd = shm_open(args->shm_name, O_RDWR, 0); + if ( shm_fd == -1 ) { + ERROR("SHM setup failed: %s\n", strerror(errno)); + return 1; + } + shared = mmap(NULL, sizeof(struct sb_shm), PROT_READ | PROT_WRITE, + MAP_SHARED, shm_fd, 0); + if ( shared == MAP_FAILED ) { + ERROR("SHM setup failed: %s\n", strerror(errno)); + return 1; + } - case 312 : - if (sscanf(arg, "%d", &args->iargs.peak_search.min_res) != 1) - { - ERROR("Invalid value for --min-res\n"); - return EINVAL; - } - break; + queue_sem = sem_open(args->queue_sem, 0); + if ( queue_sem == SEM_FAILED ) { + ERROR("Failed to open semaphore: %s\n", strerror(errno)); + return 1; + } - case 313 : - if (sscanf(arg, "%d", &args->iargs.peak_search.max_res) != 1) - { - ERROR("Invalid value for --max-res\n"); - return EINVAL; + /* Connect via ZMQ */ + if ( args->zmq_params.addr != NULL ) { + zmqstuff = im_zmq_connect(&args->zmq_params); + if ( zmqstuff == NULL ) { + ERROR("ZMQ setup failed.\n"); + return 1; } - break; + } - case 314 : - if (sscanf(arg, "%f", &args->iargs.peak_search.min_snr_biggest_pix) != 1) - { - ERROR("Invalid value for --max-snr-biggest-pix\n"); - return EINVAL; + if ( args->asapo_params.endpoint != NULL ) { + asapostuff = im_asapo_connect(&args->asapo_params); + if ( asapostuff == NULL ) { + ERROR("ASAP::O setup failed.\n"); + pthread_mutex_lock(&shared->totals_lock); + shared->should_shutdown = 1; + pthread_mutex_unlock(&shared->totals_lock); + return 1; } - break; + } - case 315 : - if (sscanf(arg, "%f", &args->iargs.peak_search.min_snr_peak_pix) != 1) - { - ERROR("Invalid value for --max-snr-peak-pix\n"); - return EINVAL; - } - break; + set_debug_funcs(set_last_task_sandbox, notify_alive_sandbox); - case 316 : - if (sscanf(arg, "%f", &args->iargs.peak_search.min_sig) != 1) - { - ERROR("Invalid value for --max-ssig\n"); - return EINVAL; - } - break; + mille = crystfel_mille_new_fd(args->fd_mille); - case 317 : - if (sscanf(arg, "%f", &args->iargs.peak_search.min_peak_over_neighbour) != 1) - { - ERROR("Invalid value for --max-peak-over-neighbour\n"); - return EINVAL; - } - break; - - case 318 : - args->iargs.peak_search.use_saturated = 0; - break; - - case 319 : - args->iargs.peak_search.revalidate = 0; - break; - - case 320 : - args->iargs.peak_search.half_pixel_shift = 0; - break; - - case 321 : - args->iargs.peak_search.check_hdf5_snr = 1; - break; - - case 322: - args->iargs.peak_search.peakfinder8_fast = 1; - break; - - /* ---------- Indexing ---------- */ - - case 400 : - case 'z' : - args->indm_str = strdup(arg); - break; - - case 'p' : - args->cellfile = strdup(arg); - break; - - case 401 : - /* Values in 'tols' are in frac (not %) and rad - * Conversion happens a few lines below */ - r = sscanf(arg, "%f,%f,%f,%f,%f,%f", - &args->iargs.tols[0], &args->iargs.tols[1], &args->iargs.tols[2], - &args->iargs.tols[3], &args->iargs.tols[4], &args->iargs.tols[5]); - if ( r != 6 ) { - /* Try old format */ - r = sscanf(arg, "%f,%f,%f,%f", - &args->iargs.tols[0], &args->iargs.tols[1], - &args->iargs.tols[2], &args->iargs.tols[3]); - if ( r != 4 ) { - ERROR("Invalid parameters for '--tolerance'\n"); - return EINVAL; - } - args->iargs.tols[4] = args->iargs.tols[3]; - args->iargs.tols[5] = args->iargs.tols[3]; - } + ida = image_data_arrays_new(); - /* Percent to fraction */ - args->iargs.tols[0] /= 100.0; - args->iargs.tols[1] /= 100.0; - args->iargs.tols[2] /= 100.0; - args->iargs.tols[3] = deg2rad(args->iargs.tols[3]); - args->iargs.tols[4] = deg2rad(args->iargs.tols[4]); - args->iargs.tols[5] = deg2rad(args->iargs.tols[5]); - break; - - case 402 : - args->if_checkcell = 0; - break; - - case 403 : - args->if_checkcell = 1; /* This is the default */ - break; - - case 404 : - args->if_multi = 1; - break; - - case 405 : - args->if_multi = 0; /* This is the default */ - break; - - case 406 : - args->if_retry = 0; - break; - - case 407 : - args->if_retry = 1; /* This is the default */ - break; - - case 408 : - args->if_refine = 0; - break; - - case 409 : - args->if_refine = 1; /* This is the default */ - break; - - case 410 : - args->if_peaks = 0; - break; - - case 411 : - args->if_peaks = 1; /* This is the default */ - break; - - case 412 : - ERROR("The option --no-cell-combinations is no longer used.\n"); - /* .. but we can still carry on. Results will probably be - * better than the user expected. */ - break; - - case 413 : - if (sscanf(arg, "%f", &args->iargs.wavelength_estimate) != 1) - { - ERROR("Invalid value for --wavelength-estimate\n"); - return EINVAL; - } - break; + while ( !allDone ) { - case 414 : - if (sscanf(arg, "%d", &args->iargs.n_threads) != 1) - { - ERROR("Invalid value for --max-indexer-threads\n"); - return EINVAL; - } - break; + struct pattern_args pargs; + int ser; + char *line; + size_t len; + int i; + char *event_str = NULL; + char *ser_str = NULL; + int ok = 1; + int should_shutdown; - case 415 : - if (sscanf(arg, "%f", &args->iargs.clen_estimate) != 1) - { - ERROR("Invalid value for --camera-length-estimate\n"); - return EINVAL; + /* Wait until an event is ready */ + notify_alive(); + set_last_task("wait_event"); + profile_start("wait-queue-semaphore"); + if ( sem_wait(queue_sem) != 0 ) { + ERROR("Failed to wait on queue semaphore: %s\n", + strerror(errno)); } - break; - - case 416 : - args->iargs.mille = 1; - break; + profile_end("wait-queue-semaphore"); - case 417 : - free(args->iargs.milledir); - args->iargs.milledir = strdup(arg); - break; - - case 418 : - if (sscanf(arg, "%d", &args->iargs.max_mille_level) != 1) + /* Get the event from the queue */ + set_last_task("read_queue"); + pthread_mutex_lock(&shared->totals_lock); + should_shutdown = shared->should_shutdown; + pthread_mutex_unlock(&shared->totals_lock); + pthread_mutex_lock(&shared->queue_lock); + if ( ((shared->n_events==0) && (shared->no_more)) + || should_shutdown ) { - ERROR("Invalid value for --max-mille-level\n"); - return EINVAL; - } - if ( args->iargs.max_mille_level < 0 ) { - ERROR("Invalid value for --max-mille-level\n"); - return EINVAL; - } - break; - - /* ---------- Integration ---------- */ - - case 501 : - args->iargs.int_meth = integration_method(arg, &r); - if ( r ) { - ERROR("Invalid integration method '%s'\n", arg); - return EINVAL; - } - break; - - case 502 : - if ( sscanf(arg, "%f", &args->iargs.fix_profile_r) != 1 ) { - ERROR("Invalid value for --fix-profile-radius\n"); - return EINVAL; - } - break; - - case 503 : - ERROR("The option --fix-bandwidth is no longer used.\n"); - ERROR("Set the bandwidth in the geometry file instead.\n"); - break; - - case 504 : - if ( sscanf(arg, "%f", &args->iargs.fix_divergence) != 1 ) { - ERROR("Invalid value for --fix-divergence\n"); - return EINVAL; - } - break; - - case 505 : - r = sscanf(arg, "%f,%f,%f", &args->iargs.ir_inn, - &args->iargs.ir_mid, &args->iargs.ir_out); - if ( r != 3 ) { - ERROR("Invalid parameters for '--int-radius'\n"); - return EINVAL; - } - break; - - case 506 : - if ( strcmp(arg, "random") == 0 ) { - args->iargs.int_diag = INTDIAG_RANDOM; + /* Queue is empty and no more are coming, + * or another process has initiated a shutdown. + * Either way, it's time to get out of here. */ + pthread_mutex_unlock(&shared->queue_lock); + allDone = 1; + continue; + } + if ( shared->n_events == 0 ) { + ERROR("Got the semaphore, but no events in queue!\n"); + ERROR("no_more = %i\n", shared->no_more); + pthread_mutex_unlock(&shared->queue_lock); + allDone = 1; + continue; + } + + line = strdup(shared->queue[0]); + + len = strlen(line); + assert(len > 1); + for ( i=len-1; i>0; i-- ) { + if ( line[i] == ' ' ) { + line[i] = '\0'; + ser_str = &line[i+1]; + break; + } } - - if ( strcmp(arg, "all") == 0 ) { - args->iargs.int_diag = INTDIAG_ALL; + len = strlen(line); + assert(len > 1); + for ( i=len-1; i>0; i-- ) { + if ( line[i] == ' ' ) { + line[i] = '\0'; + event_str = &line[i+1]; + break; + } } - - if ( strcmp(arg, "negative") == 0 ) { - args->iargs.int_diag = INTDIAG_NEGATIVE; + if ( (ser_str != NULL) && (event_str != NULL) ) { + if ( sscanf(ser_str, "%i", &ser) != 1 ) { + STATUS("Invalid serial number '%s'\n", + ser_str); + ok = 0; + } } + if ( !ok ) { + STATUS("Invalid event string '%s'\n", + shared->queue[0]); + ok = 0; + } + + pthread_mutex_lock(&shared->debug_lock); + memcpy(shared->last_ev[args->worker_id], shared->queue[0], + MAX_EV_LEN); + pthread_mutex_unlock(&shared->debug_lock); + + shuffle_events(shared); + pthread_mutex_unlock(&shared->queue_lock); + + if ( !ok ) continue; + + pargs.filename = strdup(line); + pargs.event = safe_strdup(event_str); + + free(line); + ok = 0; + + /* Default values */ + pargs.zmq_data = NULL; + pargs.zmq_data_size = 0; + pargs.asapo_data = NULL; + pargs.asapo_data_size = 0; + pargs.asapo_meta = NULL; + + if ( args->zmq_params.addr != NULL ) { + + profile_start("zmq-fetch"); + set_last_task("ZMQ fetch"); + pargs.zmq_data = im_zmq_fetch(zmqstuff, + &pargs.zmq_data_size); + profile_end("zmq-fetch"); + + if ( (pargs.zmq_data != NULL) + && (pargs.zmq_data_size > 15) ) ok = 1; + + /* The filename/event, which will be 'fake' values in + * this case, still came via the event queue. More + * importantly, the event queue gave us a unique + * serial number for this image. */ + + } else if ( args->asapo_params.endpoint != NULL ) { + + char *filename; + char *event; + int finished = 0; + int asapo_message_id; + + profile_start("asapo-fetch"); + set_last_task("ASAPO fetch"); + pargs.asapo_data = im_asapo_fetch(asapostuff, + &pargs.asapo_data_size, + &pargs.asapo_meta, + &filename, + &event, + &finished, + &asapo_message_id); + profile_end("asapo-fetch"); + if ( pargs.asapo_data != NULL ) { + ok = 1; + + /* ASAP::O provides a meaningful filename, which + * replaces the placeholder. */ + free(pargs.filename); + free(pargs.event); + pargs.filename = filename; + pargs.event = event; + pthread_mutex_lock(&shared->queue_lock); + shared->end_of_stream[args->worker_id] = 0; + pthread_mutex_unlock(&shared->queue_lock); + + /* We will also use ASAP::O's serial number + * instead of our own. */ + ser = asapo_message_id; + } else { + if ( finished ) { + pthread_mutex_lock(&shared->queue_lock); + shared->end_of_stream[args->worker_id] = 1; + pthread_mutex_unlock(&shared->queue_lock); + } + } - if ( strcmp(arg, "implausible") == 0 ) { - args->iargs.int_diag = INTDIAG_IMPLAUSIBLE; + } else { + ok = 1; } - if ( strcmp(arg, "strong") == 0 ) { - args->iargs.int_diag = INTDIAG_STRONG; - } + if ( ok ) { + pthread_mutex_lock(&shared->debug_lock); + shared->time_last_start[args->worker_id] = get_monotonic_seconds(); + pthread_mutex_unlock(&shared->debug_lock); + profile_start("process-image"); + process_image(&args->iargs, &pargs, st, args->worker_id, + args->worker_tmpdir, ser, + shared, asapostuff, mille, ida); + profile_end("process-image"); - r = sscanf(arg, "%i,%i,%i", &args->iargs.int_diag_h, - &args->iargs.int_diag_k, &args->iargs.int_diag_l); - if ( r == 3 ) { - args->iargs.int_diag = INTDIAG_INDICES; - } + if ( asapostuff != NULL ) { + im_asapo_finalise(asapostuff, ser); + } - if ( (args->iargs.int_diag == INTDIAG_NONE) - && (strcmp(arg, "none") != 0) ) - { - ERROR("Invalid value for --int-diag.\n"); - return EINVAL; } - break; + /* NB pargs.zmq_data, pargs.asapo_data and pargs.asapo_meta + * will be copied into the image structure, so + * that it can be queried for "header" values etc. They will + * eventually be freed by image_free() under process_image(). */ - case 507 : - if ( sscanf(arg, "%f", &args->iargs.push_res) != 1 ) { - ERROR("Invalid value for --push-res\n"); - return EINVAL; + if ( args->profile ) { + profile_print_and_reset(args->worker_id); } - args->iargs.push_res *= 1e9; /* nm^-1 -> m^-1 */ - break; - - case 508 : - args->iargs.overpredict = 1; - break; - - case 509 : - args->iargs.cell_params_only = 1; - break; - - /* ---------- Output ---------- */ - - case 601 : - args->iargs.stream_nonhits = 0; - break; - - case 602 : - add_copy_header(args, arg); - break; - case 603 : - args->iargs.stream_flags = CLEAR_BIT(args->iargs.stream_flags, - STREAM_PEAKS); - break; + free(pargs.filename); + free(pargs.event); + } - case 604 : - args->iargs.stream_flags = CLEAR_BIT(args->iargs.stream_flags, - STREAM_REFLECTIONS); - break; + stream_close(st); + free(tmp); + munmap(shared, sizeof(struct sb_shm)); + sem_close(queue_sem); - case 605 : - if ( sscanf(arg, "%d", &args->serial_start) != 1 ) { - ERROR("Invalid value for --serial-start\n"); - return EINVAL; - } - break; + image_data_arrays_free(ida); + crystfel_mille_free(mille); - case 606 : - args->harvest_file = strdup(arg); - break; + /* These are both no-ops if argument is NULL */ + im_zmq_shutdown(zmqstuff); + im_asapo_shutdown(asapostuff); - default : - return ARGP_ERR_UNKNOWN; + data_template_free(args->iargs.dtempl); + if ( pf8_data != NULL ) free_pf8_private_data(pf8_data); + cleanup_indexing(args->iargs.ipriv); + cell_free(args->iargs.cell); + cleanup_indexamajig_args(args); - } +#ifdef HAVE_FFTW + fftw_cleanup(); +#endif return 0; } @@ -898,440 +672,176 @@ int main(int argc, char *argv[]) { FILE *fh = NULL; Stream *st; - struct indexamajig_arguments args; + struct indexamajig_arguments *args; char *tmpdir; /* e.g. /tmp/indexamajig.12345 */ char *rn; /* e.g. /home/taw/indexing */ int r; int timeout = 240; - struct taketwo_options *taketwo_opts = NULL; - struct felix_options *felix_opts = NULL; - struct xgandalf_options *xgandalf_opts = NULL; - struct pinkindexer_options *pinkindexer_opts = NULL; - struct fromfile_options *fromfile_opts = NULL; - struct asdf_options *asdf_opts = NULL; double wl_from_dt; double clen_from_dt; int err = 0; + char *probed_methods = NULL; + size_t mille_fn_len; + char *mille_filename; + FILE *mille_fh; + + args = parse_indexamajig_args(argc, argv); + if ( args == NULL ) return 1; + + /* Load data template (new API) */ + args->iargs.dtempl = data_template_new_from_file(args->geom_filename); + if ( args->iargs.dtempl == NULL ) { + ERROR("Failed to read detector geometry from '%s'\n", + args->geom_filename); + return 1; + } + + /* Add any headers we need to copy */ + for ( r=0; r<args->n_copy_headers; r++ ) { + data_template_add_copy_header(args->iargs.dtempl, + args->copy_headers[r]); + } + - /* Defaults for "top level" arguments */ - args.filename = NULL; - args.geom_filename = NULL; - args.outfile = NULL; - args.temp_location = strdup("."); - args.prefix = strdup(""); - args.check_prefix = 1; - args.n_proc = 1; - args.cellfile = NULL; - args.indm_str = NULL; - args.basename = 0; - args.zmq_params.addr = NULL; - args.zmq_params.request = NULL; - args.zmq_params.n_subscriptions = 0; - args.asapo_params.endpoint = NULL; - args.asapo_params.token = NULL; - args.asapo_params.beamtime = NULL; - args.asapo_params.group_id = NULL; - args.asapo_params.source = NULL; - args.asapo_params.stream = NULL; - args.asapo_params.wait_for_stream = 0; - args.asapo_params.write_output_stream = 0; - args.asapo_params.consumer_timeout_ms = 500; - args.asapo_params.use_ack = 0; - args.cpu_pin = 0; - args.serial_start = 1; - args.if_peaks = 1; - args.if_multi = 0; - args.if_retry = 1; - args.if_refine = 1; - args.if_checkcell = 1; - args.profile = 0; - args.no_data_timeout = 60; - args.copy_headers = NULL; - args.n_copy_headers = 0; - args.harvest_file = NULL; - args.taketwo_opts_ptr = &taketwo_opts; - args.felix_opts_ptr = &felix_opts; - args.xgandalf_opts_ptr = &xgandalf_opts; - args.pinkindexer_opts_ptr = &pinkindexer_opts; - args.fromfile_opts_ptr = &fromfile_opts; - args.asdf_opts_ptr = &asdf_opts; - - /* Defaults for process_image arguments */ - args.iargs.cell = NULL; - args.iargs.peak_search.noisefilter = 0; - args.iargs.peak_search.median_filter = 0; - args.iargs.tols[0] = 0.05; /* frac (not %) */ - args.iargs.tols[1] = 0.05; /* frac (not %) */ - args.iargs.tols[2] = 0.05; /* frac (not %) */ - args.iargs.tols[3] = deg2rad(1.5); /* radians */ - args.iargs.tols[4] = deg2rad(1.5); /* radians */ - args.iargs.tols[5] = deg2rad(1.5); /* radians */ - args.iargs.peak_search.threshold = 800.0; - args.iargs.peak_search.min_sq_gradient = 100000.0; - args.iargs.peak_search.min_snr = 5.0; - args.iargs.peak_search.min_pix_count = 2; - args.iargs.peak_search.max_pix_count = 200; - args.iargs.peak_search.min_res = 0; - args.iargs.peak_search.max_res = 1200; - args.iargs.peak_search.local_bg_radius = 3; - args.iargs.peak_search.min_snr_biggest_pix = 7.0; /* peak finder 9 */ - args.iargs.peak_search.min_snr_peak_pix = 6.0; - args.iargs.peak_search.min_sig = 11.0; - args.iargs.peak_search.min_peak_over_neighbour = -INFINITY; - args.iargs.peak_search.check_hdf5_snr = 0; - args.iargs.peak_search.peakfinder8_fast = 0; - args.iargs.pf_private = NULL; - args.iargs.dtempl = NULL; - args.iargs.peak_search.method = PEAK_ZAEF; - args.iargs.peak_search.half_pixel_shift = 1; - args.iargs.peak_search.pk_inn = -1.0; - args.iargs.peak_search.pk_mid = -1.0; - args.iargs.peak_search.pk_out = -1.0; - args.iargs.ir_inn = -1.0; - args.iargs.ir_mid = -1.0; - args.iargs.ir_out = -1.0; - args.iargs.peak_search.use_saturated = 1; - args.iargs.peak_search.revalidate = 1; - args.iargs.stream_flags = STREAM_PEAKS | STREAM_REFLECTIONS; - args.iargs.stream_nonhits = 1; - args.iargs.int_diag = INTDIAG_NONE; - args.iargs.min_peaks = 0; - args.iargs.overpredict = 0; - args.iargs.cell_params_only = 0; - args.iargs.wait_for_file = 0; - args.iargs.ipriv = NULL; /* No default */ - args.iargs.int_meth = integration_method("rings-nocen-nosat-nograd", NULL); - args.iargs.push_res = +INFINITY; - args.iargs.highres = +INFINITY; - args.iargs.fix_profile_r = -1.0; - args.iargs.fix_divergence = -1.0; - args.iargs.no_image_data = 0; - args.iargs.no_mask_data = 0; - args.iargs.wavelength_estimate = NAN; - args.iargs.clen_estimate = NAN; - args.iargs.n_threads = 1; - args.iargs.data_format = DATA_SOURCE_TYPE_UNKNOWN; - args.iargs.mille = 0; - args.iargs.milledir = strdup("."); - args.iargs.max_mille_level = 99; - - argp_program_version_hook = show_version; - - static char doc[] = "Index and integrate snapshot diffraction images.\v" - "For more information including a tutorial, visit " - "https://www.desy.de/~twhite/crystfel"; - - static struct argp_option options[] = { - - {NULL, 0, 0, OPTION_DOC, "Basic options:", 2}, - - {NULL, 'h', NULL, OPTION_HIDDEN, NULL}, - {NULL, 'v', NULL, OPTION_HIDDEN, NULL}, - - {"input", 'i', "infile", 0, "List of input image filenames"}, - {"output", 'o', "filename.stream", 0, "Output stream filename"}, - {"geometry",'g', "experiment.geom", 0, "Detector geometry filename"}, - {"prefix", 'x', "/path/to/images/", OPTION_NO_USAGE, "Prefix filenames from input " - "file"}, - {NULL, 'j', "nproc", 0, "Run this many analyses in parallel, default 1"}, - {"basename", 201, NULL, OPTION_NO_USAGE, "Remove director parts from the " - "filenames"}, - {"no-check-prefix", 202, NULL, OPTION_NO_USAGE, "Don't attempt to correct the " - "--prefix"}, - {"highres", 203, "res", OPTION_NO_USAGE, "Absolute resolution cutoff in Angstroms"}, - {"profile", 204, NULL, OPTION_NO_USAGE, "Show timing data for performance " - "monitoring"}, - {"temp-dir", 205, "path", OPTION_NO_USAGE, "Location for temporary folder"}, - {"wait-for-file", 206, "seconds", OPTION_NO_USAGE, "Wait for each file before " - "processing"}, - {"zmq-input", 207, "addr", OPTION_NO_USAGE, "Receive data over ZeroMQ from " - "this location"}, - {"no-image-data", 208, NULL, OPTION_NO_USAGE, "Do not load image data"}, - {"spectrum-file", 209, "fn", OPTION_NO_USAGE | OPTION_HIDDEN, - "File containing radiation spectrum"}, - {"no-mask-data", 210, NULL, OPTION_NO_USAGE, "Do not load mask data"}, - {"zmq-subscribe", 211, "tag", OPTION_NO_USAGE, "Subscribe to ZMQ message" - "type"}, - {"zmq-request", 212, "str", OPTION_NO_USAGE, "Request messages using" - "this string."}, - {"asapo-endpoint", 213, "str", OPTION_NO_USAGE, "ASAP::O endpoint"}, - {"asapo-token", 214, "str", OPTION_NO_USAGE, "ASAP::O token"}, - {"asapo-beamtime", 215, "str", OPTION_NO_USAGE, "ASAP::O beamtime ID"}, - {"asapo-group", 217, "str", OPTION_NO_USAGE, "ASAP::O group ID"}, - {"asapo-source", 218, "str", OPTION_NO_USAGE, "ASAP::O data source"}, - {"data-format", 219, "str", OPTION_NO_USAGE, "Streamed data format"}, - {"asapo-stream", 220, "str", OPTION_NO_USAGE, "ASAP::O stream name"}, - {"asapo-wait-for-stream", 221, NULL, OPTION_NO_USAGE, - "Wait for ASAP::O stream to appear"}, - {"asapo-output-stream", 222, NULL, OPTION_NO_USAGE, - "Create an ASAP::O hits-only stream"}, - {"cpu-pin", 223, NULL, OPTION_NO_USAGE, "Pin worker processes to CPUs"}, - {"no-data-timeout", 224, "s", OPTION_NO_USAGE, - "Shut down after this many seconds without ASAP::O data"}, - {"asapo-consumer-timeout", 225, "ms", OPTION_NO_USAGE, - "ASAP::O get_next timeout for one frame (milliseconds)"}, - {"asapo-acks", 226, NULL, OPTION_NO_USAGE, "Use ASAP::O acknowledgements"}, - - {NULL, 0, 0, OPTION_DOC, "Peak search options:", 3}, - {"peaks", 301, "method", 0, "Peak search method. Default: zaef"}, - {"peak-radius", 302, "r1,r2,r3", OPTION_NO_USAGE, "Radii for peak search"}, - {"min-peaks", 303, "n", OPTION_NO_USAGE, "Minimum number of peaks for indexing"}, - {"hdf5-peaks", 304, "p", OPTION_HIDDEN, "Location of peak table in HDF5 file"}, - {"median-filter", 305, "n", OPTION_NO_USAGE, "Apply median filter to image data"}, - {"filter-noise", 306, NULL, OPTION_NO_USAGE, "Apply noise filter to image data"}, - {"threshold", 't', "adu", OPTION_NO_USAGE, "Threshold for peak detection " - "(zaef only, default 800)"}, - {"min-squared-gradient", 307, "n", OPTION_NO_USAGE, "Minimum squared gradient " - "(zaef only, default 100000)"}, - {"min-gradient", 307, "n", OPTION_ALIAS | OPTION_HIDDEN, NULL}, - {"min-snr", 308, "n", OPTION_NO_USAGE, "Minimum signal/noise ratio for peaks " - "(zaef,peakfinder8,peakfinder9 only, default 5)"}, - {"min-pix-count", 309, "n", OPTION_NO_USAGE, "Minimum number of pixels per peak " - "(peakfinder8 only, default 2)"}, - {"max-pix-count", 310, "n", OPTION_NO_USAGE, "Maximum number of pixels per peak " - "(peakfinder8 only, default 2)"}, - {"local-bg-radius", 311, "n", OPTION_NO_USAGE, "Radius (pixels) for local " - "background estimation (peakfinder8/9 only, default 3)"}, - {"min-res", 312, "n", OPTION_NO_USAGE, "Minimum resoultion (pixels) for peak " - "search (peakfinder8 only, default 0)"}, - {"max-res", 313, "n", OPTION_NO_USAGE, "Maximum resoultion (pixels) for peak " - "search (peakfinder8 only, default 1200)"}, - {"min-snr-biggest-pix", 314, "n", OPTION_NO_USAGE, "Minimum SNR of the biggest " - "pixel in the peak (peakfinder9 only)"}, - {"min-snr-peak-pix", 315, "n", OPTION_NO_USAGE, "Minimum SNR of peak pixel " - "(peakfinder9 only)"}, - {"min-sig", 316, "n", OPTION_NO_USAGE, "Minimum standard deviation of the " - "background (peakfinder9 only)"}, - {"min-peak-over-neighbour", 317, "n", OPTION_NO_USAGE, "Minimum difference between " - "highest pixel and neighbours (peakfinder9 only, just for speed)"}, - {"no-use-saturated", 318, NULL, OPTION_NO_USAGE, "Reject saturated peaks"}, - {"no-revalidate", 319, NULL, OPTION_NO_USAGE, "Don't re-integrate and check HDF5 " - "or MsgPack peaks"}, - {"no-half-pixel-shift", 320, NULL, OPTION_NO_USAGE, "Don't offset HDF5 peak " - "locations by 0.5 pixels"}, - {"check-hdf5-snr", 321, NULL, OPTION_NO_USAGE, "Check SNR for peaks from HDF5, " - "CXI or MsgPack (see --min-snr)"}, - {"peakfinder8-fast", 322, NULL, OPTION_NO_USAGE, "peakfinder8 fast execution"}, - - {NULL, 0, 0, OPTION_DOC, "Indexing options:", 4}, - {"indexing", 400, "method", 0, "List of indexing methods"}, - {NULL, 'z', "method", OPTION_HIDDEN | OPTION_ALIAS, NULL}, - {"pdb", 'p', "parameters.cell", 0, "PDB or CrystFEL Unit Cell File"}, - {"tolerance", 401, "a,b,c,al,be,ga", OPTION_NO_USAGE, "Tolerances for cell " - "comparison in percent and degrees, default 5,5,5,1.5,1.5,1.5"}, - {"no-check-cell", 402, NULL, OPTION_NO_USAGE, "Don't check cell parameters " - "against target cell"}, - {"check-cell", 403, NULL, OPTION_HIDDEN, NULL}, - {"multi", 404, NULL, OPTION_NO_USAGE, "Repeat indexing to index multiple hits"}, - {"no-multi", 405, NULL, OPTION_HIDDEN, NULL}, - {"no-retry", 406, NULL, OPTION_NO_USAGE, "Don't repeat indexing to increase " - "indexing rate"}, - {"retry", 407, NULL, OPTION_HIDDEN, NULL}, - {"no-refine", 408, NULL, OPTION_NO_USAGE, "Skip prediction refinement"}, - {"refine", 409, NULL, OPTION_HIDDEN, NULL}, - {"no-check-peaks", 410, NULL, OPTION_NO_USAGE, "Don't check that most peaks can be " - "accounted for by the indexing solution"}, - {"check-peaks", 411, NULL, OPTION_HIDDEN, NULL}, - {"no-cell-combinations", 412, NULL, OPTION_HIDDEN, NULL}, - {"wavelength-estimate", 413, "metres", 0, - "Estimate of the incident radiation wavelength, in metres."}, - {"max-indexer-threads", 414, "n", 0, - "Maximum number of threads allowed for indexing engines."}, - {"camera-length-estimate", 415, "metres", 0, - "Estimate of the camera length, in metres."}, - {"mille", 416, NULL, 0, - "Generate data for detector geometry refinement using Millepede"}, - {"mille-dir", 417, "dirname", 0, "Save Millepede data in folder"}, - {"max-mille-level", 418, "n", 0, "Maximum geometry refinement level"}, - - {NULL, 0, 0, OPTION_DOC, "Integration options:", 5}, - {"integration", 501, "method", OPTION_NO_USAGE, "Integration method"}, - {"fix-profile-radius", 502, "r", OPTION_NO_USAGE, "Fix profile radius for spot " - "prediction, instead of automatically determining"}, - {"fix-bandwidth", 503, "bw", OPTION_NO_USAGE, "Set the bandwidth for spot " - "prediction"}, - {"fix-divergence", 504, "deg", OPTION_NO_USAGE, "Set the divergence (full angle) " - "for spot prediction"}, - {"int-radius", 505, "r1,r2,r3", 0, "Set the integration radii (inner,mid,outer)"}, - {"int-diag", 506, "condition", 0, "Show debugging information about reflections"}, - {"push-res", 507, "dist", 0, "Integrate higher than apparent resolution cutoff (m^-1)"}, - {"overpredict", 508, NULL, 0, "Over-predict reflections"}, - {"cell-parameters-only", 509, NULL, 0, "Don't predict reflections at all"}, - - {NULL, 0, 0, OPTION_DOC, "Output options:", 6}, - {"no-non-hits-in-stream", 601, NULL, OPTION_NO_USAGE, "Don't include non-hits in " - "stream (see --min-peaks)"}, - {"copy-hdf5-field", 602, "f", OPTION_HIDDEN, NULL}, - {"copy-header", 602, "f", OPTION_NO_USAGE, "Put the value of this image header " - "field into the stream"}, - {"no-peaks-in-stream", 603, NULL, OPTION_NO_USAGE, "Don't put peak search results " - "in stream"}, - {"no-refls-in-stream", 604, NULL, OPTION_NO_USAGE, "Don't put integration results " - "in stream"}, - {"serial-start", 605, "n", OPTION_NO_USAGE, "Start the serial numbers in the stream " - "here"}, - {"harvest-file", 606, "filename", OPTION_NO_USAGE, "Write the actual parameters " - "used in JSON format"}, - - {NULL, 0, 0, OPTION_DOC, "More information:", 99}, - - {0} - - }; - - static struct argp_child argp_children[] = { - {&taketwo_argp, 0, NULL, -2}, - {&felix_argp, 0, NULL, -2}, - {&xgandalf_argp, 0, NULL, -2}, - {&pinkIndexer_argp, 0, NULL, -2}, - {&fromfile_argp, 0, NULL, -2}, - {&asdf_argp, 0, NULL, -2}, - {0} - }; - - #if defined(__APPLE__) - int q; - for ( q=0; q<argc; q++ ) { - if ( strcmp(argv[q], "--help") == 0 ) { - fprintf(stderr, "\n" - "WARNING: 'indexamajig --help' crashes on some Mac OS versions.\n" - "This is a known problem, and is due to a bug in external library " - "code.\n\n"); - fflush(stderr); - break; + wl_from_dt = data_template_get_wavelength_if_possible(args->iargs.dtempl); + if ( !isnan(wl_from_dt) ) { + if ( !isnan(args->iargs.wavelength_estimate) && !args->worker) { + ERROR("WARNING: Ignoring your value for " + "--wavelength-estimate because the geometry file " + "already contains a static value.\n"); } + args->iargs.wavelength_estimate = wl_from_dt; } - #endif - static struct argp argp = { options, parse_arg, NULL, doc, - argp_children, NULL, NULL }; - if ( argp_parse(&argp, argc, argv, 0, NULL, &args) ) return 1; + clen_from_dt = data_template_get_clen_if_possible(args->iargs.dtempl); + if ( !isnan(clen_from_dt) ) { + if ( !isnan(args->iargs.clen_estimate) && !args->worker) { + ERROR("WARNING: Ignoring your value for " + "--camera-length-estimate because the geometry file " + "already contains a static value.\n"); + } + args->iargs.clen_estimate = clen_from_dt; + } + + /* If no integration radii were given, apply the defaults */ + if ( args->iargs.ir_inn < 0 ) { + args->iargs.ir_inn = 4.0; + args->iargs.ir_mid = 5.0; + args->iargs.ir_out = 7.0; + } + + /* If no peak radii were given, copy the integration radii */ + if ( args->iargs.peak_search.pk_inn < 0.0 ) { + args->iargs.peak_search.pk_inn = args->iargs.ir_inn; + args->iargs.peak_search.pk_mid = args->iargs.ir_mid; + args->iargs.peak_search.pk_out = args->iargs.ir_out; + } + + if ( args->worker ) { + /* I am a worker process */ + return run_work(args); + } /* else I am the main process */ /* Check for minimal information */ - if ( (args.filename == NULL) - && (args.zmq_params.addr == NULL) - && (args.asapo_params.endpoint == NULL) ) { + if ( (args->filename == NULL) + && (args->zmq_params.addr == NULL) + && (args->asapo_params.endpoint == NULL) ) { ERROR("You need to provide the input filename (use -i)\n"); return 1; } - if ( args.geom_filename == NULL ) { + if ( args->geom_filename == NULL ) { ERROR("You need to specify the geometry filename (use -g)\n"); return 1; } - if ( args.outfile == NULL ) { + if ( args->outfile == NULL ) { ERROR("You need to specify the output filename (use -o)\n"); return 1; } - if ( (args.filename != NULL) && (args.zmq_params.addr != NULL) ) { + if ( (args->filename != NULL) && (args->zmq_params.addr != NULL) ) { ERROR("The options --input and --zmq-input are mutually " "exclusive.\n"); return 1; } - if ( (args.filename != NULL) && (args.asapo_params.endpoint != NULL) ) { + if ( (args->filename != NULL) && (args->asapo_params.endpoint != NULL) ) { ERROR("The options --input and --asapo-endpoint are mutually " "exclusive.\n"); return 1; } - if ( (args.asapo_params.endpoint != NULL) && (args.zmq_params.addr != NULL) ) { + if ( (args->asapo_params.endpoint != NULL) && (args->zmq_params.addr != NULL) ) { ERROR("The options --asapo-endpoint and --zmq-input are mutually " "exclusive.\n"); return 1; } - if ( (args.zmq_params.request != NULL) && (args.zmq_params.n_subscriptions > 0) ) { + if ( (args->zmq_params.request != NULL) && (args->zmq_params.n_subscriptions > 0) ) { ERROR("The options --zmq-request and --zmq-subscribe are " "mutually exclusive.\n"); return 1; } - if ( (args.filename != NULL) - && (strcmp(args.filename, "-") != 0) - && is_hdf5_file(args.filename, &err) ) + if ( (args->filename != NULL) + && (strcmp(args->filename, "-") != 0) + && is_hdf5_file(args->filename, &err) ) { ERROR("Your input file appears to be an HDF5 file.\n"); ERROR("The input file should be a list of data files, not the " "data file itself.\n"); ERROR("If you have only one input file, try the following:\n"); - ERROR(" echo %s > files.lst\n", args.filename); + ERROR(" echo %s > files.lst\n", args->filename); ERROR(" indexamajig -i files.lst -o %s -g %s ...\n", - args.outfile, args.geom_filename); + args->outfile, args->geom_filename); return 1; } else if ( err ) { - ERROR("Couldn't open '%s'\n", args.filename); + ERROR("Couldn't open '%s'\n", args->filename); return 1; } /* Open input */ - if ( args.filename != NULL ) { - if ( strcmp(args.filename, "-") == 0 ) { + if ( args->filename != NULL ) { + if ( strcmp(args->filename, "-") == 0 ) { fh = stdin; } else { - fh = fopen(args.filename, "r"); + fh = fopen(args->filename, "r"); } if ( fh == NULL ) { - ERROR("Failed to open input file '%s'\n", args.filename); + ERROR("Failed to open input file '%s'\n", args->filename); return 1; } } /* Check prefix (if given) */ - if ( args.check_prefix ) { - args.prefix = check_prefix(args.prefix); + if ( args->check_prefix ) { + args->prefix = check_prefix(args->prefix); } /* Check number of processes */ - if ( args.n_proc == 0 ) { + if ( args->n_proc == 0 ) { ERROR("Invalid number of processes.\n"); return 1; } - /* Load data template (new API) */ - args.iargs.dtempl = data_template_new_from_file(args.geom_filename); - if ( args.iargs.dtempl == NULL ) { - ERROR("Failed to read detector geometry from '%s'\n", - args.geom_filename); - return 1; - } - - /* Add any headers we need to copy */ - for ( r=0; r<args.n_copy_headers; r++ ) { - data_template_add_copy_header(args.iargs.dtempl, - args.copy_headers[r]); - } - - /* If no integration radii were given, apply the defaults */ - if ( args.iargs.ir_inn < 0 ) { - STATUS("WARNING: You did not specify --int-radius.\n"); - STATUS("WARNING: I will use the default values, which are" - " probably not appropriate for your patterns.\n"); - args.iargs.ir_inn = 4.0; - args.iargs.ir_mid = 5.0; - args.iargs.ir_out = 7.0; - } - - /* If no peak radii were given, copy the integration radii */ - if ( args.iargs.peak_search.pk_inn < 0.0 ) { - args.iargs.peak_search.pk_inn = args.iargs.ir_inn; - args.iargs.peak_search.pk_mid = args.iargs.ir_mid; - args.iargs.peak_search.pk_out = args.iargs.ir_out; - } - /* Load unit cell (if given) */ - if ( args.cellfile != NULL ) { - args.iargs.cell = load_cell_from_file(args.cellfile); - if ( args.iargs.cell == NULL ) { - ERROR("Couldn't read unit cell (from %s)\n", args.cellfile); + if ( args->cellfile != NULL ) { + args->iargs.cell = load_cell_from_file(args->cellfile); + if ( args->iargs.cell == NULL ) { + ERROR("Couldn't read unit cell (from %s)\n", args->cellfile); return 1; } - free(args.cellfile); } else { - args.iargs.cell = NULL; + args->iargs.cell = NULL; + } + + if ( args->iargs.cell != NULL ) { + STATUS("This is what I understood your unit cell to be:\n"); + cell_print(args->iargs.cell); + } else { + STATUS("No reference unit cell provided.\n"); } - tmpdir = create_tempdir(args.temp_location); + tmpdir = create_tempdir(args->temp_location); if ( tmpdir == NULL ) return 1; /* Change into temporary folder, temporarily, to control the crap @@ -1345,7 +855,7 @@ int main(int argc, char *argv[]) } /* Auto-detect indexing methods if 'requested' */ - if ( args.indm_str == NULL ) { + if ( args->indm_str == NULL ) { STATUS("No indexing methods specified. I will try to "); STATUS("automatically detect the available methods.\n"); @@ -1353,104 +863,30 @@ int main(int argc, char *argv[]) STATUS("which methods to use with --indexing=<methods>.\n"); STATUS("Use --indexing=none to disable indexing and integration.\n"); - args.indm_str = detect_indexing_methods(args.iargs.cell); - - } - - wl_from_dt = data_template_get_wavelength_if_possible(args.iargs.dtempl); - if ( !isnan(wl_from_dt) ) { - if ( !isnan(args.iargs.wavelength_estimate) ) { - ERROR("WARNING: Ignoring your value for " - "--wavelength-estimate because the geometry file " - "already contains a static value.\n"); - } - args.iargs.wavelength_estimate = wl_from_dt; - } - - clen_from_dt = data_template_get_clen_if_possible(args.iargs.dtempl); - if ( !isnan(clen_from_dt) ) { - if ( !isnan(args.iargs.clen_estimate) ) { - ERROR("WARNING: Ignoring your value for " - "--camera-length-estimate because the geometry file " - "already contains a static value.\n"); + probed_methods = detect_indexing_methods(args->iargs.cell); + if ( probed_methods == NULL ) { + ERROR("No indexing method specified, and no usable indexing "); + ERROR("methods auto-detected.\n"); + ERROR("Install some indexing programs (mosflm,dirax etc), or "); + ERROR("try again with --indexing=none.\n"); + return 1; + } else { + STATUS("Auto-determined indexing methods: %s\n", probed_methods); } - args.iargs.clen_estimate = clen_from_dt; - } - /* Prepare the indexing system */ - if ( args.indm_str == NULL ) { - - ERROR("No indexing method specified, and no usable indexing "); - ERROR("methods auto-detected.\n"); - ERROR("Install some indexing programs (mosflm,dirax etc), or "); - ERROR("try again with --indexing=none.\n"); - return 1; - - } else if ( strcmp(args.indm_str, "none") == 0 ) { + } else if ( strcmp(args->indm_str, "none") == 0 ) { STATUS("Indexing/integration disabled.\n"); - if ( args.iargs.cell != NULL ) { + if ( args->iargs.cell != NULL ) { STATUS("Ignoring your unit cell.\n"); } - args.iargs.ipriv = NULL; + args->iargs.ipriv = NULL; } else { - int i, n; - const IndexingMethod *methods; - IndexingFlags flags = 0; - - if ( args.iargs.cell != NULL ) { - STATUS("This is what I understood your unit cell to be:\n"); - cell_print(args.iargs.cell); - } else { - STATUS("No reference unit cell provided.\n"); - } - - if ( args.if_checkcell ) { - flags |= INDEXING_CHECK_CELL; - } - if ( args.if_refine ) { - flags |= INDEXING_REFINE; - } - if ( args.if_peaks ) { - flags |= INDEXING_CHECK_PEAKS; - } - if ( args.if_multi ) { - flags |= INDEXING_MULTI; - } - if ( args.if_retry ) { - flags |= INDEXING_RETRY; - } - - args.iargs.ipriv = setup_indexing(args.indm_str, - args.iargs.cell, - args.iargs.tols, - flags, - args.iargs.wavelength_estimate, - args.iargs.clen_estimate, - args.iargs.n_threads, - taketwo_opts, - xgandalf_opts, - pinkindexer_opts, - felix_opts, - fromfile_opts, - asdf_opts); - - free(args.filename); - - if ( args.iargs.ipriv == NULL ) { - ERROR("Failed to set up indexing system\n"); - return 1; - } - - methods = indexing_methods(args.iargs.ipriv, &n); - for ( i=0; i<n; i++ ) { - if ( (methods[i] & INDEXING_METHOD_MASK) == INDEXING_PINKINDEXER ) { - /* Extend timeout if using pinkIndexer */ - timeout = 3000; - break; - } + if ( strstr(args->indm_str, "pinkindexer") != NULL ) { + /* Extend timeout if using pinkIndexer */ + timeout = 3000; } } @@ -1466,61 +902,60 @@ int main(int argc, char *argv[]) free(rn); /* Open output stream */ - st = stream_open_for_write(args.outfile, args.iargs.dtempl); + st = stream_open_for_write(args->outfile, args->iargs.dtempl); if ( st == NULL ) { - ERROR("Failed to open stream '%s'\n", args.outfile); + ERROR("Failed to open stream '%s'\n", args->outfile); return 1; } /* Write audit info */ stream_write_commandline_args(st, argc, argv); - stream_write_geometry_file(st, args.geom_filename); - stream_write_target_cell(st, args.iargs.cell); - stream_write_indexing_methods(st, args.indm_str); - - if ( (args.harvest_file != NULL) && (args.serial_start <= 1) ) { - write_harvest_file(&args.iargs, args.harvest_file, - args.if_multi, args.if_refine, args.if_retry, - args.if_peaks, args.if_checkcell); - } + stream_write_geometry_file(st, args->geom_filename); + stream_write_target_cell(st, args->iargs.cell); + stream_write_indexing_methods(st, args->indm_str); - free(args.outfile); - free(args.indm_str); + if ( (args->harvest_file != NULL) && (args->serial_start <= 1) ) { + write_harvest_file(&args->iargs, args->harvest_file, + args->if_multi, args->if_refine, args->if_retry, + args->if_peaks, args->if_checkcell); + } - r = mkdir(args.iargs.milledir, S_IRWXU); + r = mkdir(args->milledir, S_IRWXU); if ( r ) { if ( errno != EEXIST ) { ERROR("Failed to create folder for Millepede data: %s\n", strerror(errno)); - exit(1); + return 1; } } - - gsl_set_error_handler_off(); - - struct pf8_private_data *pf8_data = NULL; - if ( args.iargs.peak_search.method == PEAK_PEAKFINDER8 ) { - struct detgeom *detgeom = data_template_get_2d_detgeom_if_possible(args.iargs.dtempl); - pf8_data = prepare_peakfinder8(detgeom, args.iargs.peak_search.peakfinder8_fast); - args.iargs.pf_private = pf8_data; - detgeom_free(detgeom); + mille_fn_len = strlen(args->milledir) +strlen(args->millefile)+2; + mille_filename = malloc(mille_fn_len); + if ( mille_filename == NULL ) { + ERROR("Failed to generate Millepede filename\n"); + return 1; } + snprintf(mille_filename, mille_fn_len, "%s/%s", args->milledir, args->millefile); + mille_fh = fopen(mille_filename, "wb"); + if ( mille_fh == NULL ) { + ERROR("Failed to open Millepede file: %s\n", strerror(errno)); + return 1; + } + free(mille_filename); - r = create_sandbox(&args.iargs, args.n_proc, args.prefix, args.basename, - fh, st, tmpdir, args.serial_start, - &args.zmq_params, &args.asapo_params, - timeout, args.profile, args.cpu_pin, - args.no_data_timeout); + r = create_sandbox(&args->iargs, args->n_proc, args->prefix, args->basename, + fh, st, tmpdir, args->serial_start, + &args->zmq_params, &args->asapo_params, + timeout, args->profile, args->cpu_pin, + args->no_data_timeout, argc, argv, + probed_methods, mille_fh); - if ( pf8_data != NULL ) free_pf8_private_data(pf8_data); - cell_free(args.iargs.cell); - free(args.iargs.milledir); - free(args.prefix); - free(args.temp_location); + fclose(mille_fh); free(tmpdir); - data_template_free(args.iargs.dtempl); + free(probed_methods); + data_template_free(args->iargs.dtempl); + cell_free(args->iargs.cell); stream_close(st); - cleanup_indexing(args.iargs.ipriv); + cleanup_indexamajig_args(args); return r; } diff --git a/src/merge.c b/src/merge.c index 4a2dae66..14935f9a 100644 --- a/src/merge.c +++ b/src/merge.c @@ -491,7 +491,10 @@ double log_residual(RefList *list, Crystal *cr, /* Has to match run_merge_job to be useful */ -void write_unmerged(const char *fn, struct crystal_refls *crystals, int n_crystals) +void write_unmerged(const char *fn, + struct crystal_refls *crystals, + struct image **images, + int n_crystals) { FILE *fh; int i; @@ -510,6 +513,7 @@ void write_unmerged(const char *fn, struct crystal_refls *crystals, int n_crysta UnitCell *cell; fprintf(fh, "Crystal %i\n", i); + fprintf(fh, "Filename: %s %s\n", images[i]->filename, images[i]->ev); if ( crystal_get_user_flag(crystals[i].cr) ) { fprintf(fh, "Flagged: yes\n"); } else { diff --git a/src/merge.h b/src/merge.h index e304db78..4453a2e1 100644 --- a/src/merge.h +++ b/src/merge.h @@ -35,6 +35,7 @@ #endif +#include "image.h" #include "crystal.h" #include "reflist.h" #include "geometry.h" @@ -59,6 +60,9 @@ extern double residual(RefList *list, Crystal *cr, const RefList *full, int free extern double log_residual(RefList *list, Crystal *cr, const RefList *full, int free, int *pn_used, const char *filename); -extern void write_unmerged(const char *fn, struct crystal_refls *crystals, int n_crystals); +extern void write_unmerged(const char *fn, + struct crystal_refls *crystals, + struct image **images, + int n_crystals); #endif /* MERGE */ diff --git a/src/partialator.c b/src/partialator.c index 7d3f19f6..93e78abb 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -1809,7 +1809,7 @@ int main(int argc, char *argv[]) } if ( unmerged_filename != NULL ) { - write_unmerged(unmerged_filename, crystals, n_crystals); + write_unmerged(unmerged_filename, crystals, images, n_crystals); } /* Write final figures of merit (no rejection any more) */ diff --git a/src/process_image.c b/src/process_image.c index 9586d471..9bf601b9 100644 --- a/src/process_image.c +++ b/src/process_image.c @@ -106,7 +106,6 @@ static struct image *file_wait_open_read(const char *filename, const char *event, DataTemplate *dtempl, struct sb_shm *sb_shared, - char *last_task, signed int wait_for_file, int cookie, int no_image_data, @@ -119,13 +118,13 @@ static struct image *file_wait_open_read(const char *filename, int r; struct image *image; - set_last_task(last_task, "wait for file"); + set_last_task("wait for file"); do { struct stat statbuf; - sb_shared->pings[cookie]++; + notify_alive(); r = stat(filename, &statbuf); if ( r ) { @@ -152,8 +151,8 @@ static struct image *file_wait_open_read(const char *filename, do { - set_last_task(last_task, "read file"); - sb_shared->pings[cookie]++; + set_last_task("read file"); + notify_alive(); profile_start("image-read"); image = image_read(dtempl, filename, event, @@ -181,7 +180,7 @@ static struct image *file_wait_open_read(const char *filename, void process_image(const struct index_args *iargs, struct pattern_args *pargs, Stream *st, int cookie, const char *tmpdir, int serial, struct sb_shm *sb_shared, - char *last_task, struct im_asapo *asapostuff, + struct im_asapo *asapostuff, Mille *mille, ImageDataArrays *ida) { struct image *image; @@ -194,7 +193,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, if ( pargs->zmq_data != NULL ) { - set_last_task(last_task, "unpacking ZMQ data"); + set_last_task("unpacking ZMQ data"); profile_start("read-zmq-data"); image = image_read_data_block(iargs->dtempl, pargs->zmq_data, @@ -217,7 +216,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, } else if ( pargs->asapo_data != NULL ) { - set_last_task(last_task, "unpacking ASAP::O data"); + set_last_task("unpacking ASAP::O data"); profile_start("read-asapo-data"); image = image_read_data_block(iargs->dtempl, pargs->asapo_data, @@ -241,7 +240,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, profile_start("file-wait-open-read"); image = file_wait_open_read(pargs->filename, pargs->event, iargs->dtempl, - sb_shared, last_task, + sb_shared, iargs->wait_for_file, cookie, iargs->no_image_data, @@ -261,9 +260,9 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, image->serial = serial; /* Take snapshot of image before applying horrible noise filters */ - set_last_task(last_task, "image filter"); + set_last_task("image filter"); profile_start("image-filter"); - sb_shared->pings[cookie]++; + notify_alive(); if ( (iargs->peak_search.median_filter > 0) || iargs->peak_search.noisefilter ) { profile_start("data-backup"); @@ -286,11 +285,11 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, } profile_end("image-filter"); - set_last_task(last_task, "resolution range"); - sb_shared->pings[cookie]++; + set_last_task("resolution range"); + notify_alive(); mark_resolution_range_as_bad(image, iargs->highres, +INFINITY); - sb_shared->pings[cookie]++; + notify_alive(); profile_start("peak-search"); switch ( iargs->peak_search.method ) { @@ -298,7 +297,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, case PEAK_HDF5: case PEAK_CXI: - set_last_task(last_task, "peaksearch:hdf5orcxi"); + set_last_task("peaksearch:hdf5orcxi"); peaks = image_read_peaks(iargs->dtempl, pargs->filename, pargs->event, @@ -319,7 +318,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, break; case PEAK_ZAEF: - set_last_task(last_task, "peaksearch:zaef"); + set_last_task("peaksearch:zaef"); image->features = search_peaks(image, iargs->peak_search.threshold, iargs->peak_search.min_sq_gradient, iargs->peak_search.min_snr, @@ -330,7 +329,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, break; case PEAK_PEAKFINDER8: - set_last_task(last_task, "peaksearch:pf8"); + set_last_task("peaksearch:pf8"); image->features = peakfinder8(image, 2048, iargs->peak_search.threshold, iargs->peak_search.min_snr, @@ -345,7 +344,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, break; case PEAK_PEAKFINDER9: - set_last_task(last_task, "peaksearch:pf9"); + set_last_task("peaksearch:pf9"); image->features = search_peaks_peakfinder9(image, iargs->peak_search.min_snr_biggest_pix, iargs->peak_search.min_snr_peak_pix, @@ -370,6 +369,8 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, iargs->peak_search.check_hdf5_snr); image_feature_list_free(peaks); image->features = npeaks; + } else { + image->features = peaks; } break; @@ -428,10 +429,9 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, image->hit = 1; /* Index the pattern */ - set_last_task(last_task, "indexing"); + set_last_task("indexing"); profile_start("index"); - index_pattern_4(image, iargs->ipriv, &sb_shared->pings[cookie], - last_task, mille, iargs->max_mille_level); + index_pattern_5(image, iargs->ipriv, mille, iargs->max_mille_level); profile_end("index"); r = chdir(rn); @@ -442,7 +442,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, free(rn); /* Set beam/crystal parameters */ - set_last_task(last_task, "prediction params"); + set_last_task("prediction params"); if ( iargs->fix_profile_r >= 0.0 ) { for ( i=0; i<image->n_crystals; i++ ) { crystal_set_profile_radius(image->crystals[i].cr, @@ -466,9 +466,9 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, /* Integrate! */ if ( !iargs->cell_params_only ) { - set_last_task(last_task, "integration"); + set_last_task("integration"); profile_start("integration"); - sb_shared->pings[cookie]++; + notify_alive(); integrate_all_5(image, iargs->int_meth, PMODEL_XSPHERE, iargs->push_res, iargs->ir_inn, iargs->ir_mid, iargs->ir_out, @@ -479,9 +479,9 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, } streamwrite: - set_last_task(last_task, "stream write"); + set_last_task("stream write"); profile_start("write-stream"); - sb_shared->pings[cookie]++; + notify_alive(); ret = stream_write_chunk(st, image, iargs->stream_flags); if ( ret != 0 ) { ERROR("Error writing stream file.\n"); @@ -501,8 +501,8 @@ streamwrite: out: /* Count crystals which are still good */ - set_last_task(last_task, "process_image finalisation"); - sb_shared->pings[cookie]++; + set_last_task("process_image finalisation"); + notify_alive(); pthread_mutex_lock(&sb_shared->totals_lock); any_crystals = 0; for ( i=0; i<image->n_crystals; i++ ) { @@ -521,5 +521,5 @@ out: image_free(image); profile_end("free-image"); - set_last_task(last_task, "sandbox"); + set_last_task("sandbox"); } diff --git a/src/process_image.h b/src/process_image.h index 10c9ec31..e61d79ed 100644 --- a/src/process_image.h +++ b/src/process_image.h @@ -42,7 +42,6 @@ struct index_args; #endif #include "integration.h" -#include "im-sandbox.h" #include "peaks.h" #include "image.h" #include "im-asapo.h" @@ -75,7 +74,6 @@ struct index_args float clen_estimate; int n_threads; int mille; - char *milledir; int max_mille_level; /* Integration */ @@ -99,6 +97,8 @@ struct index_args }; +#include "im-sandbox.h" + /* Information about the indexing process for one pattern */ struct pattern_args { @@ -118,7 +118,7 @@ struct pattern_args extern void process_image(const struct index_args *iargs, struct pattern_args *pargs, Stream *st, int cookie, const char *tmpdir, int serial, - struct sb_shm *sb_shared, char *last_task, + struct sb_shm *sb_shared, struct im_asapo *asapostuff, Mille *mille, ImageDataArrays *ida); diff --git a/tests/file-wait b/tests/file-wait index 22e86f52..dc4efa0f 100755 --- a/tests/file-wait +++ b/tests/file-wait @@ -27,7 +27,7 @@ fi make_link & echo "file_wait_input.h5" > file_wait_input.lst -$INDEXAMAJIG -i file_wait_input.lst -o file_wait_input.stream -g $GEOM --wait-for-file=10 +$INDEXAMAJIG -i file_wait_input.lst -o file_wait_input.stream -g $GEOM --indexing=none --wait-for-file=10 OUTVAL=$? wait |