diff options
| -rw-r--r-- | .gitlab-ci.yml | 83 | ||||
| -rw-r--r-- | doc/articles/tutorial-images/detector-shift-2.png | bin | 48161 -> 0 bytes | |||
| -rw-r--r-- | doc/articles/tutorial-images/detector-shift.png | bin | 70526 -> 0 bytes | |||
| -rw-r--r-- | doc/articles/tutorial.rst | 99 | ||||
| -rw-r--r-- | doc/man/indexamajig.1.md | 3 | ||||
| -rw-r--r-- | julia/CrystFEL/src/CrystFEL.jl | 2 | ||||
| -rw-r--r-- | julia/CrystFEL/src/detgeom.jl | 140 | ||||
| -rw-r--r-- | julia/CrystFEL/src/image.jl | 7 | ||||
| -rw-r--r-- | julia/CrystFEL/src/indexing.jl | 3 | ||||
| -rw-r--r-- | julia/plot-geom.jl | 162 | ||||
| -rw-r--r-- | libcrystfel/src/cell-utils.c | 4 | ||||
| -rw-r--r-- | libcrystfel/src/detgeom.c | 1 | ||||
| -rw-r--r-- | libcrystfel/src/index.c | 7 | ||||
| -rw-r--r-- | libcrystfel/src/indexers/pinkindexer.c | 1 | ||||
| -rw-r--r-- | libcrystfel/src/integration.c | 5 | ||||
| -rw-r--r-- | libcrystfel/src/predict-refine.c | 19 | ||||
| -rw-r--r-- | libcrystfel/src/utils.c | 20 | ||||
| -rw-r--r-- | libcrystfel/src/utils.h | 1 | ||||
| -rw-r--r-- | src/im-asapo.c | 20 | ||||
| -rw-r--r-- | src/im-asapo.h | 1 | ||||
| -rw-r--r-- | src/indexamajig.c | 6 |
21 files changed, 398 insertions, 186 deletions
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 7f9e8160..92a33525 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -3,22 +3,6 @@ stages: - deploy - build_container -.build-and-test: - image: fedora:32 - script: - - dnf group install -y 'Development Tools' - - dnf install -y hdf5-devel gsl-devel flex bison gcc-c++ - - dnf install -y gtk3-devel cairo-devel pango-devel gdk-pixbuf2-devel fftw-devel libpng-devel diffutils - - dnf install -y meson pandoc - - meson setup build && ninja -C build - - ninja -C build test - artifacts: - when: always - paths: - - build/meson-logs/testlog.junit.xml - reports: - junit: build/meson-logs/testlog.junit.xml - build-and-test: image: gitlab.desy.de:5555/thomas.white/crystfel/fedora32/dev_env script: @@ -31,22 +15,6 @@ build-and-test: reports: junit: build/meson-logs/testlog.junit.xml -.build-and-test-minimal: - image: fedora:37 - script: - - dnf group install -y 'Development Tools' - - dnf install -y gsl-devel flex bison - - dnf install -y gtk3-devel cairo-devel pango-devel gdk-pixbuf2-devel fftw-devel libpng-devel diffutils - - dnf install -y meson - - meson setup build -Dhdf5=disabled && ninja -C build - - ninja -C build test - artifacts: - when: always - paths: - - build/meson-logs/testlog.junit.xml - reports: - junit: build/meson-logs/testlog.junit.xml - build-and-test-minimal: image: gitlab.desy.de:5555/thomas.white/crystfel/fedora37/dev_env script: @@ -120,28 +88,6 @@ build_container_image_tag: - if: $CI_COMMIT_TAG when: manual -.build-native-macos: - tags: - - macOS - variables: - GIT_STRATEGY: clone - script: - - brew update - - brew upgrade -f - - brew install gsl hdf5 flex bison argp-standalone pkg-config doxygen gtk+3 cairo pango gdk-pixbuf fftw meson ninja - - export PATH="$(brew --prefix)/opt/bison/bin:$(brew --prefix)/opt/flex/bin:$PATH" - - export LDFLAGS="-L$(brew --prefix)/opt/bison/lib -L$(brew --prefix)/opt/flex/lib -L$(brew --prefix)/opt/argp-standalone/lib -largp $LDFLAGS" - - export CFLAGS="-I$(brew --prefix)/opt/flex/include -I$(brew --prefix)/opt/argp-standalone/include/ $CFLAGS" - - meson setup build - - ninja -C build - - ninja -C build test - artifacts: - when: always - paths: - - build/meson-logs/testlog.junit.xml - reports: - junit: build/meson-logs/testlog.junit.xml - build-brew-macos: dependencies: [] tags: @@ -155,31 +101,8 @@ build-brew-macos: - brew install -v -s ./crystfel.rb --HEAD - brew test crystfel - brew uninstall -v -f crystfel - -.build-centos7: - image: centos:centos7 - #image:centos:latest - #image: centos:centos7.9.2009 - variables: - SECURE_FILES_DOWNLOAD_PATH: './downloads' - script: - - yum update -y - - yum group install -y "Development Tools" - - yum install -y wget lz4 lz4-devel bzip2-libs bzip2-devel python3 gtk3-devel curl which libcurl-devel pandoc - - curl --silent "https://gitlab.com/gitlab-org/incubation-engineering/mobile-devops/download-secure-files/-/raw/main/installer" | bash - - python3 -m pip install meson ninja - - wget https://gitlab.desy.de/thomas.white/install-crystfel/-/raw/main/run-me-desy-maxwell-centos7 - - chmod 755 ./run-me-desy-maxwell-centos7 - - ./run-me-desy-maxwell-centos7 --yes-really - - cp -Rp /software $CI_PROJECT_DIR - artifacts: - paths: - - software/ - - temp/crystfel/build/meson-logs/testlog.junit.xml - when: on_success - expire_in: 1 day - reports: - junit: temp/crystfel/build/meson-logs/testlog.junit.xml + rules: + - if: $CI_COMMIT_BRANCH == "master" build-centos7: image: gitlab.desy.de:5555/thomas.white/crystfel/centos7/dev_env @@ -203,7 +126,7 @@ build-centos7: reports: junit: temp/crystfel/build/meson-logs/testlog.junit.xml rules: - - if: $CI_PIPELINE_SOURCE != 'merge_request_event' + - if: $CI_COMMIT_BRANCH == "master" build-deploy-local-maxwell: tags: diff --git a/doc/articles/tutorial-images/detector-shift-2.png b/doc/articles/tutorial-images/detector-shift-2.png Binary files differdeleted file mode 100644 index 3c17148f..00000000 --- a/doc/articles/tutorial-images/detector-shift-2.png +++ /dev/null diff --git a/doc/articles/tutorial-images/detector-shift.png b/doc/articles/tutorial-images/detector-shift.png Binary files differdeleted file mode 100644 index 1bd50738..00000000 --- a/doc/articles/tutorial-images/detector-shift.png +++ /dev/null diff --git a/doc/articles/tutorial.rst b/doc/articles/tutorial.rst index 15c23ae0..9e02b2f3 100644 --- a/doc/articles/tutorial.rst +++ b/doc/articles/tutorial.rst @@ -760,69 +760,44 @@ parameteters for just one frame, before launching the big processing job! Advanced: Check and optimise the detector geometry -------------------------------------------------- -This step is marked as *Advanced* because it needs some work outside the -CrystFEL GUI, but it's a very good idea to do this check. - Even if the geometry file is supposedly correct for the experiment, it's best -to check that, for example, the beam position hasn't drifted. Fortunately, -CrystFEL has already done most of the work for you. After indexing each -pattern, CrystFEL runs a short optimisation procedure, adjusting the unit cell -parameters, orientation and beam centre position to get the best possible -agreement between the observed and predicted peak locations at the same time as -making sure that the observed peaks correspond to Bragg positions as closely as -possible. The required beam shift (equivalently considered as a shift of the -detector) is recorded in the stream for each pattern (see the *Advanced* part -of section 8):: - - $ grep "predict_refine/det_shift" index-all-cell-1/crystfel.stream | head -n 5 - predict_refine/det_shift x = 0.020 y = 0.010 mm - predict_refine/det_shift x = -0.028 y = -0.017 mm - predict_refine/det_shift x = 0.028 y = 0.048 mm - predict_refine/det_shift x = 0.022 y = -0.014 mm - predict_refine/det_shift x = -0.016 y = -0.002 mm - $ - -A program provided with CrystFEL called detector-shift will plot these values. -Simply run it on the stream:: - - $ detector-shift index-all-cell-1/crystfel.stream - Mean shifts: dx = 0.0062 mm, dy = -0.0097 mm - -You can also run ``detector-shift`` conveniently from within the CrystFEL GUI: -Simply choose **Check detector shift** from the **Tools** menu. This will run -``detector-shift`` on whichever indexing result is currently selected (see the -start of section 7). - -A window should open, which shows the detector shifts as a scatter plot: - -.. image:: tutorial-images/detector-shift.png - -The cyan point marks the origin (0,0), and the pink point marks the mean of all -the offsets. The pixel size of the `Jungfrau detector -<https://www.psi.ch/en/lxn/jungfrau>`_, which was used for this experiment, is -75 µm, so almost all of the offsets are less than 1 pixel, and the average -offset is very much less than 1 pixel. Therefore, no further refinement is -required. - -Just for reference, here is how the graph might look if the offset were larger: - -.. image:: tutorial-images/detector-shift-2.png - -Notice that the cluster of points is significantly displaced from the origin. -This offset has already been taken into account by CrystFEL when calculating -the position of Bragg peaks, but the results will be better overall if the -geometry is correct from the start of the process. In this case, it would be a -good idea to update the geometry file. The detector-shift program can fix the -geometry file for you:: - - $ detector-shift index-all-cell-1/crystfel.stream jf-swissfel-16M.geom - Mean shifts: dx = -0.14 mm, dy = -0.25 mm - Applying corrections to jf-swissfel-16M.geom, output filename jf-swissfel-16M-predrefine.geom - default res 9097.525473 - -The updated geometry file is called ``jf-swissfel-16M-predrefine.geom``, as the -script tells you. If you process this dataset again, use this new geometry -file. +to check that, for example, the beam position hasn't drifted. CrystFEL +includes a tool for refining the detector geometry. This tool looks for small +alterations to the geometry which improve the average fit between observed and +calculated peak locations, for all patterns, while taking into account that +changing the detector position will also change the orientation reported by the +indexing procedure. The topic of detector geometry refinement is deep and +complex, but only a basic application is needed here. + +Click **Refine detector geometry** in the bar at the left of the CrystFEL GUI. +You will be presented with a dialogue box to choose the filename for the new, +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 +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 +comprising the detector are refined individually, but this usually needs more +data than have been included in this tutorial. + +Once you click **Save**, the updated geometry file should appear after a moment +of calculation. The corrections applied to the detector position are reported +in the terminal:: + + Millepede succeeded. + + Group all: + x-translation +0.019464 mm + y-translation +0.001548 mm + +If you process this dataset again, use this new geometry file and you should +find that the indexing results are slightly better. + +The command-line tool for detector geometry refinement is called +``align_detector``, and you can read more about how it works in its +`manual page <../man/align_detector.1.md>`_. diff --git a/doc/man/indexamajig.1.md b/doc/man/indexamajig.1.md index 1955ef76..42d555e7 100644 --- a/doc/man/indexamajig.1.md +++ b/doc/man/indexamajig.1.md @@ -359,6 +359,9 @@ BASIC OPTIONS : Set the timeout used for "get next" calls from ASAP::O, in ms. The default : is 500 ms. +**--asapo-acks* +: Use ASAP::O acknowledgements for more reliable message delivery. + **--data-format=format** : Specify the data format for data received over ZeroMQ or ASAP::O. Possible : values in this version are msgpack, hdf5 and seedee. diff --git a/julia/CrystFEL/src/CrystFEL.jl b/julia/CrystFEL/src/CrystFEL.jl index 9ab5d15f..bb1f1534 100644 --- a/julia/CrystFEL/src/CrystFEL.jl +++ b/julia/CrystFEL/src/CrystFEL.jl @@ -38,7 +38,7 @@ export rotatecell, compare_reindexed_cell_parameters include("detgeom.jl") using .DetGeoms -export Panel, DetGeom +export Panel, DetGeom, DetGeomPanel, DetGeomGroup, findgroup include("symmetry.jl") using .Symmetry diff --git a/julia/CrystFEL/src/detgeom.jl b/julia/CrystFEL/src/detgeom.jl index ba943d98..d8dc5972 100644 --- a/julia/CrystFEL/src/detgeom.jl +++ b/julia/CrystFEL/src/detgeom.jl @@ -1,7 +1,8 @@ module DetGeoms -export Panel +export DetGeom, DetGeomPanel, DetGeomGroup, InternalDetGeom, copydetgeom +export findgroup -mutable struct Panel +mutable struct InternalPanel name::Cstring cx::Cdouble cy::Cdouble @@ -21,27 +22,144 @@ mutable struct Panel end -mutable struct DetGeom - panels::Ptr{Panel} +mutable struct InternalDetGeomPanelGroup + name::Cstring + n_children::Cint + parent::Ptr{InternalDetGeomPanelGroup} + serial::Cint + cx::Cdouble + cy::Cdouble + cz::Cdouble + children::Ptr{Ptr{InternalDetGeomPanelGroup}} + panel::Ptr{InternalPanel} +end + + +mutable struct InternalDetGeom + panels::Ptr{InternalPanel} n_panels::Cint - top_group::Ptr{Cvoid} + top_group::Ptr{InternalDetGeomPanelGroup} end -function Base.show(io::IO, p::Panel) +mutable struct DetGeomPanel + name + corner + fsvec + ssvec + pixel_pitch + adu_per_photon + max_adu + w + h +end + +struct DetGeomGroup + serial + name + children + center # For rotation +end + + +mutable struct DetGeom + panels::Vector{DetGeomPanel} + topgroup::DetGeomGroup +end + + +function copy_dg_panel(p::InternalPanel) + DetGeomPanel(unsafe_string(p.name), + [p.cx, p.cy, p.cz], + [p.fsx, p.fsy, p.fsz], + [p.ssx, p.ssy, p.ssz], + p.pixel_pitch, p.adu_per_photon, p.max_adu, + p.w, p.h) +end + + +function copydggroup(panels, grp) + grpdata = unsafe_load(grp, 1) + name = unsafe_string(grpdata.name) + children = DetGeomGroup[] + if grpdata.n_children > 1 + for i in 1:grpdata.n_children + cdata = unsafe_load(grpdata.children, i) + push!(children, copydggroup(panels, cdata)) + end + DetGeomGroup(grpdata.serial, name, children, + [grpdata.cx, grpdata.cy, grpdata.cz]) + else + name = unsafe_string(grpdata.name) + pn = findfirst(x->x.name==name, panels) + DetGeomGroup(grpdata.serial, name, [panels[pn]], + [grpdata.cx, grpdata.cy, grpdata.cz]) + end +end + + +function copydetgeom(dg::Ptr{InternalDetGeom}) + dgdata = unsafe_load(dg, 1) + panels = DetGeomPanel[] + for i in 1:dgdata.n_panels + paneldata = unsafe_load(dgdata.panels, i) + push!(panels, copy_dg_panel(paneldata)) + end + topgroup = copydggroup(panels, dgdata.top_group) + DetGeom(panels, topgroup) +end + + +function seq100(ser, head) + if ser == 0 + () + elseif ser < 100 + (head...,ser % 100) + else + seq100(ser ÷ 100, (head..., ser % 100)) + end +end + +findgroup(dg::DetGeom, ser) = findgroup(dg.topgroup, ser, seq100(ser÷100, ())) + +function findgroup(group::DetGeomGroup, ser, path) + if length(path) == 0 + @assert(group.serial == ser) + return group + end + findgroup(group.children[path[1]], ser, path[2:end]) +end + + +function Base.show(io::IO, p::DetGeomPanel) write(io, "Panel(") write(io, "name=\"") - write(io, unsafe_string(p.name)) - write(io, "\", center=(") - show(io, p.cx); write(io, ", "); show(io, p.cy); write(io, ", "); show(io, p.cz) + write(io, p.name) + write(io, "\", corner=") + show(io, p.corner) write(io, "), fs=(") - show(io, p.fsx); write(io, ", "); show(io, p.fsy); write(io, ", "); show(io, p.fsz) + show(io, p.fsvec) write(io, "), ss=(") - show(io, p.ssx); write(io, ", "); show(io, p.ssy); write(io, ", "); show(io, p.ssz) + show(io, p.ssvec) write(io, "), size=(") show(io, p.w); write(io, ", "); show(io, p.h) write(io, "))") end +function showgroup(io, p::DetGeomPanel, prefix) + println(io, prefix, "Panel ", p.name, " (",p.w,"×",p.h,")") +end + +function showgroup(io, grp::DetGeomGroup, prefix) + println(io, prefix, "Group ", grp.name, " (serial number ", grp.serial, ")") + for child in grp.children + showgroup(io, child, prefix*" ") + end +end + +function Base.show(io::IO, ::MIME"text/plain", dg::DetGeom) + println(io, "Detector geometry structure with ", length(dg.panels), " panels") + showgroup(io, dg.topgroup, " ") +end end # of module diff --git a/julia/CrystFEL/src/image.jl b/julia/CrystFEL/src/image.jl index d648dba4..cd24ce01 100644 --- a/julia/CrystFEL/src/image.jl +++ b/julia/CrystFEL/src/image.jl @@ -4,7 +4,7 @@ using Printf import ..CrystFEL: libcrystfel import ..CrystFEL.DataTemplates: DataTemplate, InternalDataTemplate -import ..CrystFEL.DetGeoms: DetGeom +import ..CrystFEL.DetGeoms: InternalDetGeom, copydetgeom import ..CrystFEL.PeakLists: PeakList, InternalPeakList import ..CrystFEL.Crystals: Crystal, InternalCrystal import ..CrystFEL.RefLists: RefList, InternalRefList, UnmergedReflection @@ -29,7 +29,7 @@ mutable struct InternalImage n_crystals::Cint indexed_by::Cint n_indexing_tries::Cint - detgeom::Ptr{DetGeom} + detgeom::Ptr{InternalDetGeom} data_source_type::Cint filename::Cstring ev::Cstring @@ -135,6 +135,9 @@ function Base.getproperty(image::Image, name::Symbol) getfield(idata, :crystals), getfield(idata, :n_crystals)) + elseif name === :detgeom + return copydetgeom(getfield(idata, :detgeom)) + else getfield(idata, name) end diff --git a/julia/CrystFEL/src/indexing.jl b/julia/CrystFEL/src/indexing.jl index 543246ec..81a59b52 100644 --- a/julia/CrystFEL/src/indexing.jl +++ b/julia/CrystFEL/src/indexing.jl @@ -144,7 +144,8 @@ function index(image::Image, idxr::Indexer; mille=nothing) idxr.indexingpriv::Ptr{IndexingPriv}, C_NULL::Ptr{Cvoid}, C_NULL::Ptr{Cvoid}, - imille::Ptr{Cvoid})::Cvoid + imille::Ptr{Cvoid}, + 99::Cint)::Cvoid end diff --git a/julia/plot-geom.jl b/julia/plot-geom.jl new file mode 100644 index 00000000..754111e0 --- /dev/null +++ b/julia/plot-geom.jl @@ -0,0 +1,162 @@ +using CrystFEL +using GLMakie + +function readmillepede(filename) + fh = open(filename, "r") + readline(fh) # discard header line + motions = map(eachline(fh)) do line + sgroupser,sdelta,_ = split(line) + (parse(Int,sgroupser), parse(Float64,sdelta)) + end + return motions +end + +function findeigenvector(fh, num) + while !eof(fh) + l = readline(fh) + bits = split(l) + if length(bits) == 5 + if bits[1] == "Eigenvector" + if parse(Int, bits[2]) == num + return + end + end + end + end + error("Eigenvector not found") +end + +function readeigenvector(filename, num) + fh = open(filename, "r") + findeigenvector(fh, num) + motions = [] + while !eof(fh) + l = readline(fh) + length(l) < 3 && return motions + bits = split(l) + for i in 1:2:length(bits) + push!(motions, (parse(Int,bits[i]), parse(Float64,bits[i+1]))) + end + end + return motions +end + +function translate!(panel::DetGeomPanel, vec) + panel.corner += vec/panel.pixel_pitch +end + +function translate!(group::DetGeomGroup, vec) + for child in group.children + translate!(child, vec) + end +end + +function rotate2d!(p, c, ang) + nx = c[1] + (p[1]-c[1])*cos(ang) - (p[2]-c[2])*sin(ang); + ny = c[2] + (p[1]-c[1])*sin(ang) + (p[2]-c[2])*cos(ang); + p[1] = nx + p[2] = ny +end + +function rotatex!(panel::DetGeomPanel, center, angle) + rotate2d!(view(panel.corner, 2:3), center[2:3]/panel.pixel_pitch, angle); + rotate2d!(view(panel.fsvec, 2:3), [0,0], angle); + rotate2d!(view(panel.ssvec, 2:3), [0,0], angle); +end + +function rotatey!(panel::DetGeomPanel, center, angle) + rotate2d!(view(panel.corner, 3:-2:1), center[3:-2:1]/panel.pixel_pitch, angle); + rotate2d!(view(panel.fsvec, 3:-2:1), [0,0], angle); + rotate2d!(view(panel.ssvec, 3:-2:1), [0,0], angle); +end + +function rotatez!(panel::DetGeomPanel, center, angle) + rotate2d!(view(panel.corner, 1:2), center[1:2]/panel.pixel_pitch, angle); + rotate2d!(view(panel.fsvec, 1:2), [0,0], angle); + rotate2d!(view(panel.ssvec, 1:2), [0,0], angle); +end + +function rotatex!(group::DetGeomGroup, center, angle) + for child in group.children + rotatex!(child, center, angle) + end +end + +function rotatey!(group::DetGeomGroup, center, angle) + for child in group.children + rotatey!(child, center, angle) + end +end + +function rotatez!(group::DetGeomGroup, center, angle) + for child in group.children + rotatez!(child, center, angle) + end +end + +function plotgeom!(dtempl, motions) + + f = Figure(size=(1024,1024)) + ax = Axis3(f[1,1], aspect=:equal, limits=(nothing,nothing,nothing,nothing,0,0.3), + xlabel="x (meters)", ylabel="y (meters)", zlabel="z (meters)") + sl = Slider(f[2,1], range=0:0.01:1) + btoggle = Toggle(f, halign=:left, tellwidth=false) + f[3,1] = grid!([1,1]=>Label(f, "Show beams", halign=:right, tellwidth=false), [1,2]=>btoggle) + + cfimage = CrystFEL.Image(dtempl) + triangles = lift(sl.value) do trscale + + dgmoved = cfimage.detgeom # Makes a fresh copy of underlying C structure + + # Apply all `motions` in the same way as align_detector + for m in motions + param = m[1] % 100 + group = findgroup(dgmoved, m[1]-param) + if param == 1 + translate!(group, trscale*[-m[2],0,0]) + elseif param == 2 + translate!(group, trscale*[0,-m[2],0]) + elseif param == 3 + translate!(group, trscale*[0,0,-m[2]]) + elseif param == 4 + rotatex!(group, group.center, trscale*-m[2]) + elseif param == 5 + rotatey!(group, group.center, trscale*-m[2]) + elseif param == 6 + rotatez!(group, group.center, trscale*-m[2]) + end + end + + # Generate triangle mesh from panel descriptions + points = Point3f[] + for panel in dgmoved.panels + + corner = panel.pixel_pitch*panel.corner + fs = panel.pixel_pitch*panel.w*panel.fsvec + ss = panel.pixel_pitch*panel.h*panel.ssvec + + push!(points, corner) + push!(points, corner+fs) + push!(points, corner+fs+ss) + push!(points, corner+fs+ss) + push!(points, corner+ss) + push!(points, corner) + + end + + return points + + end + + beams = map(cfimage.detgeom.panels) do p + cen = (p.corner+p.fsvec*p.w/2+p.ssvec*p.h/2)*1.5*p.pixel_pitch + ((0,0,0), Tuple(cen)) + end + + mesh!(triangles, color="#7f7f7f") + beamv = lines!(collect(Iterators.flatten(beams[1:4:end])), color="#008080") + connect!(beamv.visible, btoggle.active) + + return f + +end diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c index 80cdd955..a25accbc 100644 --- a/libcrystfel/src/cell-utils.c +++ b/libcrystfel/src/cell-utils.c @@ -2245,9 +2245,9 @@ static IntegerMatrix *check_permutations(UnitCell *cell_reduced, UnitCell *refer best_diff = diff; n_best = 0; - } + best_m[n_best++] = intmat_copy(m); - if ( diff < 1.001*best_diff ) { + } else if ( diff < 1.001*best_diff ) { /* If the new solution is the same as the * previous one, add it to the list */ diff --git a/libcrystfel/src/detgeom.c b/libcrystfel/src/detgeom.c index 39b00663..2b38664c 100644 --- a/libcrystfel/src/detgeom.c +++ b/libcrystfel/src/detgeom.c @@ -269,6 +269,7 @@ gsl_matrix **make_panel_minvs(struct detgeom *dg) gsl_matrix_set(M, 2, 2, p->pixel_pitch*p->ssz); Minvs[i] = matrix_invert(M); + gsl_matrix_free(M); if ( Minvs[i] == NULL ) { ERROR("Failed to calculate inverse panel matrix for %s\n", p->name); diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 12cd9190..5b0f742f 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -50,6 +50,7 @@ #include "cell-utils.h" #include "predict-refine.h" #include "profile.h" +#include "crystfel-mille.h" #include "indexers/dirax.h" #include "indexers/asdf.h" #include "indexers/mosflm.h" @@ -774,6 +775,12 @@ static int try_indexer(struct image *image, IndexingMethod indm, } profile_end("cell-compare-to-others"); + if ( (mille != NULL) && !crystal_get_user_flag(cr) ) { + profile_start("mille-write"); + crystfel_mille_write_record(mille); + profile_end("mille-write"); + } + } n_bad = remove_flagged_crystals(image); diff --git a/libcrystfel/src/indexers/pinkindexer.c b/libcrystfel/src/indexers/pinkindexer.c index 180246c2..14e18548 100644 --- a/libcrystfel/src/indexers/pinkindexer.c +++ b/libcrystfel/src/indexers/pinkindexer.c @@ -33,7 +33,6 @@ #include <stdlib.h> -#include <sys/errno.h> #include <argp.h> #include "utils.h" diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c index d2c06db7..6d37116b 100644 --- a/libcrystfel/src/integration.c +++ b/libcrystfel/src/integration.c @@ -1467,11 +1467,6 @@ int integrate_rings_once(Reflection *refl, * appropriate size */ if ( intensity < -aduph ) { sig2_poisson = -aduph*intensity; - } else if ( intensity < 0.0 ) { - /* If the intensity is negative (by less than one - * photon), assume the reflection is very weak and - * therefore has a Poisson error of one photon. */ - sig2_poisson = aduph; } sigma = sqrt(sig2_poisson + bx->m*sig2_bg); diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index 6da0d3d1..cc740db3 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -811,7 +811,7 @@ int refine_prediction(struct image *image, Crystal *cr, reflist = reflist_new(); n = pair_peaks(image, cr, reflist, rps); - if ( n < 10 ) { + if ( n < 3 ) { cffree(rps); reflist_free(reflist); return 1; @@ -872,6 +872,7 @@ int refine_prediction(struct image *image, Crystal *cr, crystal_add_notes(cr, tmp); if ( mille != NULL ) { + crystfel_mille_delete_last_record(mille); profile_start("mille-calc"); write_mille(mille, n, crystal_get_cell(cr), rps, image, max_mille_level, Minvs); @@ -887,18 +888,10 @@ int refine_prediction(struct image *image, Crystal *cr, n = pair_peaks(image, cr, NULL, rps); free_rps_noreflist(rps, n); - if ( n < 10 ) { - if ( mille != NULL ) { - crystfel_mille_delete_last_record(mille); - } - return 1; - } - if ( mille != NULL ) { - profile_start("mille-write"); - crystfel_mille_write_record(mille); - profile_end("mille-write"); + if ( n < 3 ) { + return 1; + } else { + return 0; } - - return 0; } diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c index 7bfac487..0ff09b58 100644 --- a/libcrystfel/src/utils.c +++ b/libcrystfel/src/utils.c @@ -162,12 +162,30 @@ gsl_matrix *matrix_invert(gsl_matrix *m) } gsl_permutation_free(perm); - gsl_matrix_free(m); return inv; } +gsl_vector *solve_inv(gsl_vector *v, gsl_matrix *M) +{ + gsl_matrix *Minv; + gsl_vector *s; + int n; + + n = v->size; + if ( v->size != M->size1 ) return NULL; + if ( v->size != M->size2 ) return NULL; + + Minv = matrix_invert(M); + if ( Minv == NULL ) return NULL; + + s = gsl_vector_calloc(n); + gsl_blas_dgemv(CblasNoTrans, 1.0, M, v, 0.0, s); + return s; +} + + static int check_eigen(gsl_vector *e_val, int verbose) { int i; diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index 30a7c4ef..08e059a7 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -75,6 +75,7 @@ extern "C" { extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v); extern void show_matrix(gsl_matrix *M); extern void show_vector(gsl_vector *M); +extern gsl_vector *solve_inv(gsl_vector *v, gsl_matrix *M); extern gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt, int verbose); extern gsl_matrix *matrix_mult2(gsl_matrix *A, gsl_matrix *B); diff --git a/src/im-asapo.c b/src/im-asapo.c index 6df14e93..9b57883c 100644 --- a/src/im-asapo.c +++ b/src/im-asapo.c @@ -55,6 +55,7 @@ struct im_asapo AsapoProducerHandle producer; AsapoStringHandle group_id; int wait_for_stream; + int use_ack; }; @@ -173,7 +174,10 @@ struct im_asapo *im_asapo_connect(struct im_asapo_params *params) a->group_id = asapo_string_from_c_str(params->group_id); a->wait_for_stream = params->wait_for_stream; - asapo_consumer_set_resend_nacs(a->consumer, 1, 10000, 3); + a->use_ack = params->use_ack; + if ( a->use_ack ) { + asapo_consumer_set_resend_nacs(a->consumer, 1, 60000, 3); + } asapo_free_handle(&cred); @@ -275,13 +279,15 @@ void *im_asapo_fetch(struct im_asapo *a, size_t *pdata_size, void im_asapo_finalise(struct im_asapo *a, uint64_t message_id) { - AsapoErrorHandle err = asapo_new_handle(); - asapo_consumer_acknowledge(a->consumer, a->group_id, message_id, - a->stream, &err); - if ( asapo_is_error(err) ) { - show_asapo_error("Couldn't acknowledge ASAP::O message", err); + if ( a->use_ack ) { + AsapoErrorHandle err = asapo_new_handle(); + asapo_consumer_acknowledge(a->consumer, a->group_id, message_id, + a->stream, &err); + if ( asapo_is_error(err) ) { + show_asapo_error("Couldn't acknowledge ASAP::O message", err); + } + asapo_free_handle(&err); } - asapo_free_handle(&err); } diff --git a/src/im-asapo.h b/src/im-asapo.h index d06f9527..70f67fd0 100644 --- a/src/im-asapo.h +++ b/src/im-asapo.h @@ -45,6 +45,7 @@ struct im_asapo_params int wait_for_stream; int write_output_stream; int consumer_timeout_ms; + int use_ack; }; struct im_asapo; diff --git a/src/indexamajig.c b/src/indexamajig.c index 9416d234..00b15b55 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -474,6 +474,10 @@ static error_t parse_arg(int key, char *arg, struct argp_state *state) } break; + case 226 : + args->asapo_params.use_ack = 1; + break; + /* ---------- Peak search ---------- */ case 't' : @@ -932,6 +936,7 @@ int main(int argc, char *argv[]) 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; @@ -1065,6 +1070,7 @@ int main(int argc, char *argv[]) "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"}, |
