aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel')
-rw-r--r--libcrystfel/src/cell-utils.c4
-rw-r--r--libcrystfel/src/crystfel-mille.c30
-rw-r--r--libcrystfel/src/detgeom.c1
-rw-r--r--libcrystfel/src/index.c7
-rw-r--r--libcrystfel/src/indexers/pinkindexer.c1
-rw-r--r--libcrystfel/src/integration.c5
-rw-r--r--libcrystfel/src/predict-refine.c31
-rw-r--r--libcrystfel/src/utils.c26
-rw-r--r--libcrystfel/src/utils.h1
9 files changed, 73 insertions, 33 deletions
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/crystfel-mille.c b/libcrystfel/src/crystfel-mille.c
index 72f89425..c2fe2a57 100644
--- a/libcrystfel/src/crystfel-mille.c
+++ b/libcrystfel/src/crystfel-mille.c
@@ -72,6 +72,10 @@ struct mille
int *int_arr;
int max_entries;
int n;
+
+ int *have_local;
+ int n_local;
+
FILE *fh;
};
@@ -115,6 +119,17 @@ static void mille_add_measurement(Mille *m,
m->max_entries = new_max_entries;
}
+ if ( NLC > m->n_local ) {
+ int i;
+ int *new_have_local = cfrealloc(m->have_local, NLC*sizeof(int));
+ if ( new_have_local == NULL ) return;
+ m->have_local = new_have_local;
+ for ( i=m->n_local; i<NLC; i++ ) {
+ m->have_local[i] = 0;
+ }
+ m->n_local = NLC;
+ }
+
/* The measurement */
m->float_arr[m->n] = rMeas;
m->int_arr[m->n] = 0;
@@ -126,6 +141,7 @@ static void mille_add_measurement(Mille *m,
m->float_arr[m->n] = derLc[i];
m->int_arr[m->n] = i+1;
m->n++;
+ m->have_local[i] = 1;
}
}
@@ -292,6 +308,8 @@ Mille *crystfel_mille_new(const char *outFileName)
m->n = 0;
m->float_arr = NULL;
m->int_arr = NULL;
+ m->have_local = NULL;
+ m->n_local = 0;
m->fh = fopen(outFileName, "wb");
if ( m->fh == NULL ) {
@@ -318,11 +336,15 @@ void crystfel_mille_free(Mille *m)
void crystfel_mille_delete_last_record(Mille *m)
{
m->n = 0;
+ cffree(m->have_local);
+ m->have_local = NULL;
+ m->n_local = 0;
}
void crystfel_mille_write_record(Mille *m)
{
+ int i;
float nf = 0.0;
int ni = 0;
int nw = (m->n * 2)+2;
@@ -330,6 +352,14 @@ void crystfel_mille_write_record(Mille *m)
/* Don't write empty records */
if ( m->n == 0 ) return;
+ /* Don't write records with incomplete local gradients */
+ for ( i=0; i<m->n_local; i++ ) {
+ if ( !m->have_local[i] ) {
+ crystfel_mille_delete_last_record(m);
+ return;
+ }
+ }
+
fwrite(&nw, sizeof(int), 1, m->fh);
fwrite(&nf, sizeof(float), 1, m->fh);
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..083d86c5 100644
--- a/libcrystfel/src/predict-refine.c
+++ b/libcrystfel/src/predict-refine.c
@@ -142,15 +142,15 @@ double r_gradient(int param, Reflection *refl, UnitCell *cell, double wavelength
/* Spot position gradients for diffraction physics (anything that changes the
* diffracted ray direction) */
int fs_ss_gradient_physics(int param, Reflection *refl, UnitCell *cell,
- struct detgeom_panel *p, gsl_matrix *Minv,
- double fs, double ss, double mu,
- float *fsg, float *ssg)
+ struct detgeom_panel *p, gsl_matrix *Minv,
+ double fs, double ss, double mu,
+ float *fsg, float *ssg)
{
signed int h, k, l;
gsl_vector *dRdp;
gsl_vector *v;
- get_indices(refl, &h, &k, &l);
+ get_symmetric_indices(refl, &h, &k, &l);
dRdp = gsl_vector_calloc(3);
@@ -302,7 +302,7 @@ int fs_ss_gradient(int param, Reflection *refl, UnitCell *cell,
double mu;
double fs, ss;
- get_indices(refl, &h, &k, &l);
+ get_symmetric_indices(refl, &h, &k, &l);
kpred = get_kpred(refl);
cell_get_reciprocal(cell, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
@@ -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;
@@ -871,7 +871,8 @@ int refine_prediction(struct image *image, Crystal *cr,
total_shifts[0], total_shifts[1], total_shifts[2]);
crystal_add_notes(cr, tmp);
- if ( mille != NULL ) {
+ if ( (mille != NULL) && (n>4) ) {
+ 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..a35ba272 100644
--- a/libcrystfel/src/utils.c
+++ b/libcrystfel/src/utils.c
@@ -135,39 +135,53 @@ gsl_matrix *matrix_invert(gsl_matrix *m)
perm = gsl_permutation_alloc(m->size1);
if ( perm == NULL ) {
ERROR("Couldn't allocate permutation\n");
- gsl_matrix_free(m);
return NULL;
}
inv = gsl_matrix_alloc(m->size1, m->size2);
if ( inv == NULL ) {
ERROR("Couldn't allocate inverse\n");
- gsl_matrix_free(m);
gsl_permutation_free(perm);
return NULL;
}
if ( gsl_linalg_LU_decomp(m, perm, &s) ) {
ERROR("Couldn't decompose matrix\n");
- gsl_matrix_free(m);
gsl_permutation_free(perm);
return NULL;
}
if ( gsl_linalg_LU_invert(m, perm, inv) ) {
- ERROR("Couldn't invert cell matrix:\n");
- gsl_matrix_free(m);
+ ERROR("Couldn't invert matrix:\n");
gsl_permutation_free(perm);
return NULL;
}
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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-26 15:25:23 +0000