aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-08-12 10:41:30 +0200
committerThomas White <taw@physics.org>2013-09-13 14:31:23 +0200
commitd05eecd314661fad45db4fb48df179f861d77d24 (patch)
treee53976a387b240d59008fab2e8ef7fe39c0091d9
parent896ced2c5826b71cf270d75310b8ac0fd495063b (diff)
partialator: Debugging stuff
-rw-r--r--src/partialator.c25
-rw-r--r--src/post-refinement.c22
-rw-r--r--src/post-refinement.h2
-rw-r--r--src/scaling-report.c3
4 files changed, 45 insertions, 7 deletions
diff --git a/src/partialator.c b/src/partialator.c
index 88912226..ad6d9a49 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -309,6 +309,16 @@ static void display_progress(int n_images, int n_crystals)
}
+static const char *str_flags(Crystal *cr)
+{
+ if ( crystal_get_user_flag(cr) ) {
+ return "N";
+ }
+
+ return "-";
+}
+
+
int main(int argc, char *argv[])
{
int c;
@@ -671,6 +681,21 @@ int main(int argc, char *argv[])
/* Output results */
write_reflist(outfile, full);
+ /* Dump parameters */
+ FILE *fh;
+ fh = fopen("partialator.params", "w");
+ if ( fh == NULL ) {
+ ERROR("Couldn't open partialator.params!\n");
+ } else {
+ for ( i=0; i<n_crystals; i++ ) {
+ fprintf(fh, "%4i %5.2f %8.5e %s\n", i,
+ crystal_get_osf(crystals[i]),
+ crystal_get_image(crystals[i])->div,
+ str_flags(crystals[i]));
+ }
+ fclose(fh);
+ }
+
/* Clean up */
for ( i=0; i<n_crystals; i++ ) {
reflist_free(crystal_get_reflections(crystals[i]));
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 94c8a816..88c609c1 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -314,13 +314,12 @@ static void apply_shift(Crystal *cr, int k, double shift)
}
-static int check_eigen(gsl_vector *e_val)
+static int check_eigen(gsl_matrix *M, gsl_vector *e_val, int verbose)
{
int i;
double vmax, vmin;
const int n = e_val->size;
const double max_condition = 1e6;
- const int verbose = 0;
int n_filt = 0;
if ( verbose ) STATUS("Eigenvalues:\n");
@@ -374,6 +373,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt,
gsl_matrix *AS;
gsl_matrix *SAS;
int i;
+ gsl_matrix *SAS_copy;
n = v->size;
if ( v->size != M->size1 ) return NULL;
@@ -393,6 +393,18 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt,
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, S, AS, 0.0, SAS);
gsl_matrix_free(AS);
+ /* Calculate the vector SB, which is the RHS of the equation */
+ SB = gsl_vector_calloc(n);
+ gsl_blas_dgemv(CblasNoTrans, 1.0, S, v, 0.0, SB);
+
+ if ( verbose ) {
+ STATUS("The equation after rescaling:\n");
+ show_matrix_eqn(SAS, SB);
+ }
+
+ SAS_copy = gsl_matrix_alloc(SAS->size1, SAS->size2);
+ gsl_matrix_memcpy(SAS_copy, SAS);
+
/* Do the SVD */
s_val = gsl_vector_calloc(n);
s_vec = gsl_matrix_calloc(n, n);
@@ -408,11 +420,9 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt,
/* "SAS" is now "U" */
/* Filter the eigenvalues */
- *n_filt = check_eigen(s_val);
+ *n_filt = check_eigen(SAS_copy, s_val, verbose);
- /* Calculate the vector SB, which is the RHS of the equation */
- SB = gsl_vector_calloc(n);
- gsl_blas_dgemv(CblasNoTrans, 1.0, S, v, 0.0, SB);
+ gsl_matrix_free(SAS_copy);
/* Solve the equation SAS.SinvX = SB */
SinvX = gsl_vector_calloc(n);
diff --git a/src/post-refinement.h b/src/post-refinement.h
index b964975e..63b6ff66 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -56,8 +56,8 @@ enum {
REF_CSX,
REF_CSY,
REF_CSZ,
- NUM_PARAMS,
REF_DIV,
+ NUM_PARAMS,
REF_R,
};
diff --git a/src/scaling-report.c b/src/scaling-report.c
index acb007ae..17268489 100644
--- a/src/scaling-report.c
+++ b/src/scaling-report.c
@@ -889,6 +889,9 @@ void sr_iteration(SRContext *sr, int iteration, struct srdata *d)
cairo_restore(sr->cr);
}
+
+ STATUS("%i filtered total, %5.2f filtered per crystal\n",
+ d->n_filtered, (double)d->n_filtered / d->n);
}