diff --git a/image_analysis/structure_refinement/RefinementMinimiser.h b/image_analysis/structure_refinement/RefinementMinimiser.h index 42bb3afe1..82ac9e118 100644 --- a/image_analysis/structure_refinement/RefinementMinimiser.h +++ b/image_analysis/structure_refinement/RefinementMinimiser.h @@ -27,6 +27,7 @@ public: void Reset() { s_.clear(); y_.clear(); + sy_.clear(); } // The estimated diagonal of the Hessian; every element > 0. void SetDiagonal(std::vector h) { h_ = std::move(h); } @@ -94,9 +95,11 @@ public: if (sy > 1e-10 * std::sqrt(Dot(s, s) * Dot(y, y))) { s_.push_back(std::move(s)); y_.push_back(std::move(y)); + sy_.push_back(sy); if (s_.size() > memory_) { s_.pop_front(); y_.pop_front(); + sy_.pop_front(); } } const double decrease = f - f_new; @@ -124,7 +127,7 @@ private: const size_t n = g.size(), m = s_.size(); std::vector q = g, a(m), rho(m); for (size_t i = m; i-- > 0;) { - rho[i] = 1.0 / Dot(s_[i], y_[i]); + rho[i] = 1.0 / sy_[i]; a[i] = rho[i] * Dot(s_[i], q); for (size_t k = 0; k < n; k++) q[k] -= a[i] * y_[i][k]; @@ -134,7 +137,7 @@ private: double yhy = 0; for (size_t k = 0; k < n; k++) yhy += y_[m - 1][k] * y_[m - 1][k] / h_[k]; - gamma = Dot(s_[m - 1], y_[m - 1]) / yhy; + gamma = sy_[m - 1] / yhy; } for (size_t k = 0; k < n; k++) q[k] *= gamma / h_[k]; @@ -149,5 +152,6 @@ private: size_t memory_; std::deque> s_, y_; + std::deque sy_; // s.y of each pair, as taken when it was kept std::vector h_; }; diff --git a/image_analysis/structure_refinement/RefinementXray.cpp b/image_analysis/structure_refinement/RefinementXray.cpp index 1ce87e07d..3cf0506bb 100644 --- a/image_analysis/structure_refinement/RefinementXray.cpp +++ b/image_analysis/structure_refinement/RefinementXray.cpp @@ -282,7 +282,14 @@ gemmi::Grid BulkSolventMask::MaskGrid(const gemmi::UnitCell &cell, const grid.unit_cell = cell; grid.spacegroup = &sg; const double spacing = std::max(std::min(d_min / 3.0, MASK_SPACING_MAX), std::cbrt(cell.volume / MASK_POINTS_MAX)); - grid.set_size_from_spacing(spacing, gemmi::GridSizeRounding::Up); + // set_size_from_spacing()'s size, without allocating the values + const std::array n = gemmi::good_grid_size({{cell.a / spacing, cell.b / spacing, cell.c / spacing}}, + gemmi::GridSizeRounding::Up, &sg); + grid.nu = n[0]; + grid.nv = n[1]; + grid.nw = n[2]; + grid.calculate_spacing(); + grid.axis_order = gemmi::AxisOrder::XYZ; return grid; } @@ -291,6 +298,7 @@ void BulkSolventMask::Fmask(const gemmi::Model &model, const std::vector mask = grid_; + mask.data.resize(mask.point_count()); PutMaskOnGrid(mask, model, orbit_leaders_, nthreads_); const gemmi::FPhiGrid fm = MapToFPhi(mask); f.clear();