Skip to content

Commit

Permalink
Make fewer copies of arrays (#2114)
Browse files Browse the repository at this point in the history
  • Loading branch information
@lballabio
lballabio authored Nov 19, 2024
2 parents a07bf7b + 0fb7e68 commit a83be9e
Show file tree
Hide file tree
Showing 4 changed files with 48 additions and 45 deletions.
28 changes: 13 additions & 15 deletions ql/math/optimization/levenbergmarquardt.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,17 +36,16 @@ namespace QuantLib {
EndCriteria::Type LevenbergMarquardt::minimize(Problem& P,
const EndCriteria& endCriteria) {
P.reset();
Array x_ = P.currentValue();
const Array& initX = P.currentValue();
currentProblem_ = &P;
initCostValues_ = P.costFunction().values(x_);
initCostValues_ = P.costFunction().values(initX);
int m = initCostValues_.size();
int n = x_.size();
if(useCostFunctionsJacobian_) {
int n = initX.size();
if (useCostFunctionsJacobian_) {
initJacobian_ = Matrix(m,n);
P.costFunction().jacobian(initJacobian_, x_);
P.costFunction().jacobian(initJacobian_, initX);
}
std::unique_ptr<Real[]> xx(new Real[n]);
std::copy(x_.begin(), x_.end(), xx.get());
Array xx = initX;
std::unique_ptr<Real[]> fvec(new Real[m]);
std::unique_ptr<Real[]> diag(new Real[n]);
int mode = 1;
Expand Down Expand Up @@ -81,15 +80,15 @@ namespace QuantLib {
// in n variables by the Levenberg-Marquardt algorithm.
MINPACK::LmdifCostFunction lmdifCostFunction =
[this](const auto m, const auto n, const auto x, const auto fvec, const auto iflag) {
this->fcn(m, n, x, fvec, iflag);
this->fcn(m, n, x, fvec);
};
MINPACK::LmdifCostFunction lmdifJacFunction =
useCostFunctionsJacobian_
? [this](const auto m, const auto n, const auto x, const auto fjac, const auto iflag) {
this->jacFcn(m, n, x, fjac, iflag);
this->jacFcn(m, n, x, fjac);
}
: MINPACK::LmdifCostFunction();
MINPACK::lmdif(m, n, xx.get(), fvec.get(),
MINPACK::lmdif(m, n, xx.begin(), fvec.get(),
endCriteria.functionEpsilon(),
xtol_,
gtol_,
Expand Down Expand Up @@ -132,14 +131,13 @@ namespace QuantLib {
QL_FAIL("unknown MINPACK result: " << info);
}
// set problem
std::copy(xx.get(), xx.get()+n, x_.begin());
P.setCurrentValue(x_);
P.setFunctionValue(P.costFunction().value(x_));
P.setCurrentValue(std::move(xx));
P.setFunctionValue(P.costFunction().value(P.currentValue()));

return ecType;
}

void LevenbergMarquardt::fcn(int, int n, Real* x, Real* fvec, int*) {
void LevenbergMarquardt::fcn(int, int n, Real* x, Real* fvec) {
Array xt(n);
std::copy(x, x+n, xt.begin());
// constraint handling needs some improvement in the future:
Expand All @@ -152,7 +150,7 @@ namespace QuantLib {
}
}

void LevenbergMarquardt::jacFcn(int m, int n, Real* x, Real* fjac, int*) {
void LevenbergMarquardt::jacFcn(int m, int n, Real* x, Real* fjac) {
Array xt(n);
std::copy(x, x+n, xt.begin());
// constraint handling needs some improvement in the future:
Expand Down
26 changes: 15 additions & 11 deletions ql/math/optimization/levenbergmarquardt.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,23 +55,27 @@ namespace QuantLib {
EndCriteria::Type minimize(Problem& P,
const EndCriteria& endCriteria) override;

void fcn(int m,
int n,
Real* x,
Real* fvec,
int* iflag);
void jacFcn(int m,
int n,
Real* x,
Real* fjac,
int* iflag);

/*! \deprecated Don't use this method; inspect the result of minimize instead.
Deprecated in version 1.36.
*/
[[deprecated("Don't use this method; inspect the result of minimize instead")]]
virtual Integer getInfo() const { return info_; }

/*! \deprecated Don't use this method; it is for internal use.
Deprecated in version 1.37.
*/
[[deprecated("Don't use this method; it is for internal use")]]
void fcn(int m, int n, Real* x, Real* fvec, int*) { fcn(m, n, x, fvec); }

/*! \deprecated Don't use this method; it is for internal use.
Deprecated in version 1.37.
*/
[[deprecated("Don't use this method; it is for internal use")]]
void jacFcn(int m, int n, Real* x, Real* fjac, int*) { jacFcn(m, n, x, fjac); }
private:
void fcn(int m, int n, Real* x, Real* fvec);
void jacFcn(int m, int n, Real* x, Real* fjac);

Problem* currentProblem_;
Array initCostValues_;
Matrix initJacobian_;
Expand Down
4 changes: 2 additions & 2 deletions ql/math/optimization/problem.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,8 +73,8 @@ namespace QuantLib {
//! Cost function
CostFunction& costFunction() const { return costFunction_; }

void setCurrentValue(const Array& currentValue) {
currentValue_=currentValue;
void setCurrentValue(Array currentValue) {
currentValue_ = std::move(currentValue);
}

//! current value of the local minimum
Expand Down
35 changes: 18 additions & 17 deletions ql/termstructures/globalbootstrap.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -254,15 +254,21 @@ template <class Curve> void GlobalBootstrap<Curve>::calculate() const {
// setup cost function
class TargetFunction : public CostFunction {
public:
TargetFunction(const Size firstHelper,
const Size numberHelpers,
std::function<Array()> additionalErrors,
// NOTE: TargetFunction retains passed pointers, so they must point to objects
// that outlive the call to optimizer->minimize(). We use pointers instead of
// const refs to avoid accidental binding to temporaries.
TargetFunction(Size firstHelper,
Size numberHelpers,
const std::function<Array()>* additionalErrors,
Curve* ts,
std::vector<Real> lowerBounds,
std::vector<Real> upperBounds)
const std::vector<Real>* lowerBounds,
const std::vector<Real>* upperBounds)
: firstHelper_(firstHelper), numberHelpers_(numberHelpers),
additionalErrors_(std::move(additionalErrors)), ts_(ts),
lowerBounds_(std::move(lowerBounds)), upperBounds_(std::move(upperBounds)) {}
additionalErrors_(*additionalErrors), ts_(ts),
lowerBounds_(*lowerBounds), upperBounds_(*upperBounds) {
QL_REQUIRE(additionalErrors != nullptr, "null additionalErrors");
QL_REQUIRE(lowerBounds != nullptr && upperBounds != nullptr, "null bounds");
}

Real transformDirect(const Real x, const Size i) const {
return (std::atan(x) + M_PI_2) / M_PI * (upperBounds_[i] - lowerBounds_[i]) + lowerBounds_[i];
Expand All @@ -272,12 +278,6 @@ template <class Curve> void GlobalBootstrap<Curve>::calculate() const {
return std::tan((y - lowerBounds_[i]) * M_PI / (upperBounds_[i] - lowerBounds_[i]) - M_PI_2);
}

Real value(const Array& x) const override {
Array v = values(x);
std::transform(v.begin(), v.end(), v.begin(), [](Real x) -> Real { return x*x; });
return std::sqrt(std::accumulate(v.begin(), v.end(), Real(0.0)) / static_cast<Real>(v.size()));
}

Array values(const Array& x) const override {
for (Size i = 0; i < x.size(); ++i) {
Traits::updateGuess(ts_->data_, transformDirect(x[i], i), i + 1);
Expand All @@ -300,11 +300,12 @@ template <class Curve> void GlobalBootstrap<Curve>::calculate() const {

private:
Size firstHelper_, numberHelpers_;
std::function<Array()> additionalErrors_;
Curve *ts_;
const std::vector<Real> lowerBounds_, upperBounds_;
const std::function<Array()>& additionalErrors_;
Curve* ts_;
const std::vector<Real> &lowerBounds_, &upperBounds_;
};
TargetFunction cost(firstHelper_, numberHelpers_, additionalErrors_, ts_, lowerBounds, upperBounds);
TargetFunction cost(firstHelper_, numberHelpers_, &additionalErrors_, ts_,
&lowerBounds, &upperBounds);

// setup guess
Array guess(numberBounds);
Expand Down

0 comments on commit a83be9e

Please sign in to comment.