Clean-up

image_analysis/scale_merge/Merge.cpp

@@ -293,7 +293,7 @@ MergeResult MergeReflections(const std::vector<std::vector<Reflection> > &observ
 void MergeStatistics::Print(Logger &logger) const {
     logger.Info("");
     logger.Info("  {:>8s} {:>8s} {:>8s} {:>8s}", "d_min", "N_obs", "N_uniq", "<I/sig>");
-    logger.Info("  {:->8s} {:->8s} {:->8s} {:->8s}", "", "", "", "", "", "");
+    logger.Info("  {:->8s} {:->8s} {:->8s} {:->8s}", "", "", "", "");
     for (const auto &sh: shells) {
         if (sh.unique_reflections == 0)
             continue;

image_analysis/scale_merge/ScaleOnTheFly.cpp

@@ -58,11 +58,11 @@ namespace {
                 return false;
 
             const T half_wedge = wedge[0] / T(2.0);
-            const T arg_plus = T(delta_phi_deg + half_wedge) * T(c1) / mosaicity[0];
-            const T arg_minus = T(delta_phi_deg - half_wedge) * T(c1) / mosaicity[0];
+            const T arg_plus = (T(delta_phi_deg) + half_wedge) * T(c1) / mosaicity[0];
+            const T arg_minus = (T(delta_phi_deg) - half_wedge) * T(c1) / mosaicity[0];
             const T partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0);
             const T B_term = ceres::exp(B[0] * T(b_resolution_coeff));
-            residual[0] = (G[0] * partiality * B_term * T(lp) * Itrue - T(Iobs)) * T(weight);
+            residual[0] = (G[0] * partiality * B_term * T(lp) * T(Itrue) - T(Iobs)) * T(weight);
             return true;
         }
 
@@ -126,7 +126,10 @@ ScaleOnTheFlyResult ScaleOnTheFly::Scale(std::vector<Reflection> &reflections, s
     };
 
     if (model == PartialityModel::Rotation) {
-        result.mos = mosaicity_deg.value_or(s.GetDefaultMosaicity());
+        if (mosaicity_deg && std::isfinite(*mosaicity_deg) && *mosaicity_deg > 0.0)
+            result.mos = *mosaicity_deg;
+        else
+            result.mos = s.GetDefaultMosaicity();
         result.wedge = rot_wedge_deg.value_or(0.0);
     } else {
         result.mos = NAN;
@@ -189,7 +192,7 @@ ScaleOnTheFlyResult ScaleOnTheFly::Scale(std::vector<Reflection> &reflections, s
     }
 
     if (model == PartialityModel::Rotation) {
-        if (s.GetRefineWedge()) {
+        if (refine_rot_wedge) {
             problem.SetParameterLowerBound(&result.wedge, 0, s.GetMinWedge());
             problem.SetParameterUpperBound(&result.wedge, 0, s.GetMaxWedge());
         } else {
@@ -223,7 +226,14 @@ ScaleOnTheFlyResult ScaleOnTheFly::Scale(std::vector<Reflection> &reflections, s
                 // For fixed partiality there is no need to change anything
                 break;
         }
-        r.scaling_correction = static_cast<float>(r.rlp / (B_term * r.partiality * result.G));
+        const double denom = B_term * r.partiality * result.G;
+        if (std::isfinite(r.rlp) &&
+            std::isfinite(denom) &&
+            denom > 0.0) {
+                r.scaling_correction = static_cast<float>(r.rlp / denom);
+            } else {
+                r.scaling_correction = NAN;
+            }
     }
 
     auto end = std::chrono::steady_clock::now();
@@ -242,7 +252,10 @@ ScalingResult ScaleOnTheFly::Scale(std::vector<std::vector<Reflection> > &reflec
     if (nthreads <= 1) {
         for (int i = 0; i < reflections.size(); i++) {
             std::optional<double> mos_val;
-            if (model == PartialityModel::Rotation && mosaicity.size() > i)
+            if (model == PartialityModel::Rotation
+                && mosaicity.size() > i
+                && std::isfinite(mosaicity[i])
+                && mosaicity[i] > 0.0)
                 mos_val = mosaicity[i];
 
             auto local_result = Scale(reflections[i], mos_val);
@@ -262,7 +275,10 @@ ScalingResult ScaleOnTheFly::Scale(std::vector<std::vector<Reflection> > &reflec
                 size_t i = curr_image.fetch_add(1);
                 while (i < reflections.size()) {
                     std::optional<double> mos_val;
-                    if (model == PartialityModel::Rotation && mosaicity.size() > i)
+                    if (model == PartialityModel::Rotation
+                        && mosaicity.size() > i
+                        && std::isfinite(mosaicity[i])
+                        && mosaicity[i] > 0.0)
                         mos_val = mosaicity[i];
 
                     auto local_result = Scale(reflections[i], mos_val);