11910868632

Committed 19 Nov 2024 10:14AM UTC coverage: 73.049% (-0.003%) from 73.052%

Build # 11910868632

Build Type

push

github

Committed by

web-flow

Commit Message

Make fewer copies of arrays (#2114)

Run Details

18 of 18 new or added lines in 2 files covered. (100.0%)

4 existing lines in 3 files now uncovered.

55987 of 76643 relevant lines covered (73.05%)

8661450.08 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

92.31

/ql/math/optimization/simplex.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2001, 2002, 2003 Sadruddin Rejeb
 Copyright (C) 2006 Ferdinando Ametrano
 Copyright (C) 2007 Marco Bianchetti
 Copyright (C) 2007 François du Vignaud

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

/* The implementation of the algorithm was highly inspired by
 * "Numerical Recipes in C", 2nd edition, Press, Teukolsky, Vetterling,
 * Flannery, chapter 10.
 * Modified may 2007: end criteria set on x instead on fx,
 * inspired by bad behaviour found with test function fx=x*x+x+1,
 * xStart = -100, lambda = 1.0, ftol = 1.e-16
 * (it reports x=0 as the minimum!)
 * and by GSL implementation, v. 1.9 (http://www.gnu.org/software/gsl/)
 */

#include <ql/math/optimization/simplex.hpp>
#include <ql/math/optimization/constraint.hpp>

namespace QuantLib {

    namespace {
    // Computes the size of the simplex
        Real computeSimplexSize (const std::vector<Array>& vertices) {
            Array center(vertices.front().size(),0);
            for (const auto& vertice : vertices)
                center += vertice;
            center *=1/Real(vertices.size());
            Real result = 0;
            for (const auto& vertice : vertices) {
                Array temp = vertice - center;
                result += Norm2(temp);
            }
            return result/Real(vertices.size());
        }
    }

    Real Simplex::extrapolate(Problem& P,
                              Size iHighest,
                              Real &factor) const {

        Array pTry;
        do {
            Size dimensions = values_.size() - 1;
            Real factor1 = (1.0 - factor)/dimensions;
            Real factor2 = factor1 - factor;
            pTry = sum_*factor1 - vertices_[iHighest]*factor2;
            factor *= 0.5;
        } while (!P.constraint().test(pTry) && std::fabs(factor) > QL_EPSILON);
        if (std::fabs(factor) <= QL_EPSILON) {
            return values_[iHighest];
        }
        factor *= 2.0;
        Real vTry = P.value(pTry);
        if (vTry < values_[iHighest]) {
            values_[iHighest] = vTry;
            sum_ += pTry - vertices_[iHighest];
            vertices_[iHighest] = pTry;
        }
        return vTry;

    }


    EndCriteria::Type Simplex::minimize(Problem& P,
                                        const EndCriteria& endCriteria) {
        // set up of the problem
        //Real ftol = endCriteria.functionEpsilon();    // end criteria on f(x) (see Numerical Recipes in C++, p.410)
        Real xtol = endCriteria.rootEpsilon();          // end criteria on x (see GSL v. 1.9, http://www.gnu.org/software/gsl/)
        Size maxStationaryStateIterations_
            = endCriteria.maxStationaryStateIterations();
        EndCriteria::Type ecType = EndCriteria::None;
        P.reset();

        Array x_ = P.currentValue();
        if (!P.constraint().test(x_))
            QL_FAIL("Initial guess " << x_ << " is not in the feasible region.");

        Integer iterationNumber_=0;

        // Initialize vertices of the simplex
        Size n = x_.size();
        vertices_ = std::vector<Array>(n+1, x_);
        for (Size i=0; i<n; ++i) {
            Array direction(n, 0.0);
            direction[i] = 1.0;
            P.constraint().update(vertices_[i+1], direction, lambda_);
        }
        // Initialize function values at the vertices of the simplex
        values_ = Array(n+1, 0.0);
        for (Size i=0; i<=n; ++i)
            values_[i] = P.value(vertices_[i]);
        // Loop looking for minimum
        do {
            sum_ = Array(n, 0.0);
            Size i;
            for (i=0; i<=n; i++)
                sum_ += vertices_[i];
            // Determine the best (iLowest), worst (iHighest)
            // and 2nd worst (iNextHighest) vertices
            Size iLowest = 0;
            Size iHighest, iNextHighest;
            if (values_[0]<values_[1]) {
                iHighest = 1;
                iNextHighest = 0;
            } else {
                iHighest = 0;
                iNextHighest = 1;
            }
            for (i=1;i<=n; i++) {
                if (values_[i]>values_[iHighest]) {
                    iNextHighest = iHighest;
                    iHighest = i;
                } else {
                    if ((values_[i]>values_[iNextHighest]) && i!=iHighest)
                        iNextHighest = i;
                }
                if (values_[i]<values_[iLowest])
                    iLowest = i;
            }
            // Now compute accuracy, update iteration number and check end criteria
            //// Numerical Recipes exit strategy on fx (see NR in C++, p.410)
            //Real low = values_[iLowest];
            //Real high = values_[iHighest];
            //Real rtol = 2.0*std::fabs(high - low)/
            //    (std::fabs(high) + std::fabs(low) + QL_EPSILON);
            //++iterationNumber_;
            //if (rtol < ftol ||
            //    endCriteria.checkMaxIterations(iterationNumber_, ecType)) {
            // GSL exit strategy on x (see GSL v. 1.9, http://www.gnu.org/software/gsl
            Real simplexSize = computeSimplexSize(vertices_);
            ++iterationNumber_;
            if (simplexSize < xtol ||
                endCriteria.checkMaxIterations(iterationNumber_, ecType)) {
                endCriteria.checkStationaryPoint(0.0, 0.0,
                    maxStationaryStateIterations_, ecType);
                endCriteria.checkMaxIterations(iterationNumber_, ecType);
                x_ = vertices_[iLowest];
                Real low = values_[iLowest];
                P.setFunctionValue(low);
                P.setCurrentValue(x_);
                return ecType;
            }
            // If end criteria is not met, continue
            Real factor = -1.0;
            Real vTry = extrapolate(P, iHighest, factor);
            if ((vTry <= values_[iLowest]) && (factor == -1.0)) {
                factor = 2.0;
                extrapolate(P, iHighest, factor);
            } else if (std::fabs(factor) > QL_EPSILON) {
                if (vTry >= values_[iNextHighest]) {
                    Real vSave = values_[iHighest];
                    factor = 0.5;
                    vTry = extrapolate(P, iHighest, factor);
                    if (vTry >= vSave && std::fabs(factor) > QL_EPSILON) {
                        for (Size i=0; i<=n; i++) {
                            if (i!=iLowest) {
                                vertices_[i] =
                                    0.5*(vertices_[i] + vertices_[iLowest]);
                                values_[i] = P.value(vertices_[i]);
                            }
                        }
                    }
                }
            }
            // If can't extrapolate given the constraints, exit
            if (std::fabs(factor) <= QL_EPSILON) {
                x_ = vertices_[iLowest];
                Real low = values_[iLowest];
                P.setFunctionValue(low);
                P.setCurrentValue(x_);
                return EndCriteria::StationaryFunctionValue;
            }
        } while (true);
        QL_FAIL("optimization failed: unexpected behaviour");
    }
}

1	/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3	/*
4	Copyright (C) 2001, 2002, 2003 Sadruddin Rejeb
5	Copyright (C) 2006 Ferdinando Ametrano
6	Copyright (C) 2007 Marco Bianchetti
7	Copyright (C) 2007 François du Vignaud
8
9	This file is part of QuantLib, a free-software/open-source library
10	for financial quantitative analysts and developers - http://quantlib.org/
11
12	QuantLib is free software: you can redistribute it and/or modify it
13	under the terms of the QuantLib license. You should have received a
14	copy of the license along with this program; if not, please email
15	<quantlib-dev@lists.sf.net>. The license is also available online at
16	<http://quantlib.org/license.shtml>.
17
18	This program is distributed in the hope that it will be useful, but WITHOUT
19	ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20	FOR A PARTICULAR PURPOSE. See the license for more details.
21	*/
22
23	/* The implementation of the algorithm was highly inspired by
24	* "Numerical Recipes in C", 2nd edition, Press, Teukolsky, Vetterling,
25	* Flannery, chapter 10.
26	* Modified may 2007: end criteria set on x instead on fx,
27	* inspired by bad behaviour found with test function fx=x*x+x+1,
28	* xStart = -100, lambda = 1.0, ftol = 1.e-16
29	* (it reports x=0 as the minimum!)
30	* and by GSL implementation, v. 1.9 (http://www.gnu.org/software/gsl/)
31	*/
32
33	#include <ql/math/optimization/simplex.hpp>
34	#include <ql/math/optimization/constraint.hpp>
35
36	namespace QuantLib {
37
38	namespace {
39	// Computes the size of the simplex
40	Real computeSimplexSize (const std::vector<Array>& vertices) {	243,081✔
41	Array center(vertices.front().size(),0);	243,081✔
42	for (const auto& vertice : vertices)	1,009,919✔
43	center += vertice;	766,838✔
44	center *=1/Real(vertices.size());	243,081✔
45	Real result = 0;
46	for (const auto& vertice : vertices) {	1,009,919✔
47	Array temp = vertice - center;	766,838✔
48	result += Norm2(temp);	766,838✔
49	}
50	return result/Real(vertices.size());	486,162✔
51	}
52	}
53
54	Real Simplex::extrapolate(Problem& P,	458,340✔
55	Size iHighest,
56	Real &factor) const {
57
58	Array pTry;
59	do {
60	Size dimensions = values_.size() - 1;	458,434✔
61	Real factor1 = (1.0 - factor)/dimensions;	458,434✔
62	Real factor2 = factor1 - factor;	458,434✔
63	pTry = sum_factor1 - vertices_[iHighest]factor2;	1,375,302✔
64	factor *= 0.5;	458,434✔
65	} while (!P.constraint().test(pTry) && std::fabs(factor) > QL_EPSILON);	458,434✔
66	if (std::fabs(factor) <= QL_EPSILON) {	458,340✔
67	return values_[iHighest];	×
68	}
69	factor *= 2.0;	458,340✔
70	Real vTry = P.value(pTry);	458,340✔
71	if (vTry < values_[iHighest]) {	458,340✔
72	values_[iHighest] = vTry;	97,694✔
73	sum_ += pTry - vertices_[iHighest];	195,388✔
74	vertices_[iHighest] = pTry;	97,694✔
75	}
76	return vTry;
77
78	}
79
80
81	EndCriteria::Type Simplex::minimize(Problem& P,	533✔
82	const EndCriteria& endCriteria) {
83	// set up of the problem
84	//Real ftol = endCriteria.functionEpsilon(); // end criteria on f(x) (see Numerical Recipes in C++, p.410)
85	Real xtol = endCriteria.rootEpsilon(); // end criteria on x (see GSL v. 1.9, http://www.gnu.org/software/gsl/)	533✔
86	Size maxStationaryStateIterations_
87	= endCriteria.maxStationaryStateIterations();	533✔
88	EndCriteria::Type ecType = EndCriteria::None;	533✔
89	P.reset();
90
91	Array x_ = P.currentValue();	533✔
92	if (!P.constraint().test(x_))	533✔
93	QL_FAIL("Initial guess " << x_ << " is not in the feasible region.");	×
94
95	Integer iterationNumber_=0;
96
97	// Initialize vertices of the simplex
98	Size n = x_.size();
99	vertices_ = std::vector<Array>(n+1, x_);	533✔
100	for (Size i=0; i<n; ++i) {	1,364✔
101	Array direction(n, 0.0);	831✔
102	direction[i] = 1.0;	831✔
103	P.constraint().update(vertices_[i+1], direction, lambda_);	831✔
104	}
105	// Initialize function values at the vertices of the simplex
106	values_ = Array(n+1, 0.0);	533✔
107	for (Size i=0; i<=n; ++i)	1,897✔
108	values_[i] = P.value(vertices_[i]);	1,364✔
109	// Loop looking for minimum
110	do {
111	sum_ = Array(n, 0.0);	243,081✔
112	Size i;
113	for (i=0; i<=n; i++)	1,009,919✔
114	sum_ += vertices_[i];	766,838✔
115	// Determine the best (iLowest), worst (iHighest)
116	// and 2nd worst (iNextHighest) vertices
117	Size iLowest = 0;
118	Size iHighest, iNextHighest;
119	if (values_[0]<values_[1]) {	243,081✔
120	iHighest = 1;
121	iNextHighest = 0;
122	} else {
123	iHighest = 0;
124	iNextHighest = 1;
125	}
126	for (i=1;i<=n; i++) {	766,838✔
127	if (values_[i]>values_[iHighest]) {	523,757✔
128	iNextHighest = iHighest;
129	iHighest = i;
130	} else {
131	if ((values_[i]>values_[iNextHighest]) && i!=iHighest)	459,881✔
132	iNextHighest = i;
133	}
134	if (values_[i]<values_[iLowest])	523,757✔
135	iLowest = i;
136	}
137	// Now compute accuracy, update iteration number and check end criteria
138	//// Numerical Recipes exit strategy on fx (see NR in C++, p.410)
139	//Real low = values_[iLowest];
140	//Real high = values_[iHighest];
141	//Real rtol = 2.0*std::fabs(high - low)/
142	// (std::fabs(high) + std::fabs(low) + QL_EPSILON);
143	//++iterationNumber_;
144	//if (rtol < ftol \|\|
145	// endCriteria.checkMaxIterations(iterationNumber_, ecType)) {
146	// GSL exit strategy on x (see GSL v. 1.9, http://www.gnu.org/software/gsl
147	Real simplexSize = computeSimplexSize(vertices_);	243,081✔
148	++iterationNumber_;	243,081✔
149	if (simplexSize < xtol \|\|	485,795✔
150	endCriteria.checkMaxIterations(iterationNumber_, ecType)) {	242,714✔
151	endCriteria.checkStationaryPoint(0.0, 0.0,	533✔
152	maxStationaryStateIterations_, ecType);
153	endCriteria.checkMaxIterations(iterationNumber_, ecType);	533✔
154	x_ = vertices_[iLowest];	533✔
155	Real low = values_[iLowest];	533✔
156	P.setFunctionValue(low);
157	P.setCurrentValue(x_);	533✔
158	return ecType;	533✔
159	}
160	// If end criteria is not met, continue
161	Real factor = -1.0;	242,548✔
162	Real vTry = extrapolate(P, iHighest, factor);	242,548✔
163	if ((vTry <= values_[iLowest]) && (factor == -1.0)) {	242,548✔
164	factor = 2.0;	185,594✔
165	extrapolate(P, iHighest, factor);	185,594✔
166	} else if (std::fabs(factor) > QL_EPSILON) {	56,954✔
167	if (vTry >= values_[iNextHighest]) {	56,954✔
168	Real vSave = values_[iHighest];	30,198✔
169	factor = 0.5;	30,198✔
170	vTry = extrapolate(P, iHighest, factor);	30,198✔
171	if (vTry >= vSave && std::fabs(factor) > QL_EPSILON) {	30,198✔
172	for (Size i=0; i<=n; i++) {	4,138✔
173	if (i!=iLowest) {	3,348✔
174	vertices_[i] =
175	0.5*(vertices_[i] + vertices_[iLowest]);	5,116✔
176	values_[i] = P.value(vertices_[i]);	2,558✔
177	}
178	}
179	}
180	}
181	}
182	// If can't extrapolate given the constraints, exit
183	if (std::fabs(factor) <= QL_EPSILON) {	242,548✔
184	x_ = vertices_[iLowest];	×
185	Real low = values_[iLowest];	×
186	P.setFunctionValue(low);
UNCOV 187	P.setCurrentValue(x_);	×
UNCOV 188	return EndCriteria::StationaryFunctionValue;	×
189	}
190	} while (true);	242,548✔
191	QL_FAIL("optimization failed: unexpected behaviour");
192	}
193	}

lballabio / QuantLib / 11910868632

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous