30255006850

Committed 27 Jul 2026 09:41AM UTC coverage: 74.898% (-0.05%) from 74.944%

Build # 30255006850

Build Type

push

github

Committed by

lballabio

Commit Message

Avoid deprecated Ubuntu runner

Coverage Stats

60031 of 80150 relevant lines covered (74.9%)

8637062.36 hits per line

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

85.14

/ql/math/optimization/spherecylinder.cpp

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

/*
 Copyright (C) 2007 Mark Joshi

 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
 <https://www.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.
*/

#include <ql/math/optimization/spherecylinder.hpp>
#include <ql/errors.hpp>
#include <algorithm>

namespace QuantLib {

    namespace {

        template<class F>
        Real BrentMinimize(Real low,
                           Real mid,
                           Real high,
                           Real tolerance,
                           Size maxIt,
                           const F& objectiveFunction) {
            Real W = 0.5*(3.0-std::sqrt(5.0));
            Real x = W*low+(1-W)*high;
            if (mid > low && mid < high)
                x = mid;

            Real midValue = objectiveFunction(x);

            Size iterations = 0;
            while (high-low > tolerance && iterations < maxIt) {
                if (x - low > high -x) { // left interval is bigger
                    Real tentativeNewMid = W*low+(1-W)*x;
                    Real tentativeNewMidValue =
                        objectiveFunction(tentativeNewMid);

                    if (tentativeNewMidValue < midValue) { // go left
                        high =x;
                        x = tentativeNewMid;
                        midValue = tentativeNewMidValue;
                    } else { // go right
                        low = tentativeNewMid;
                    }
                } else {
                    Real tentativeNewMid = W*x+(1-W)*high;
                    Real tentativeNewMidValue =
                        objectiveFunction(tentativeNewMid);

                    if (tentativeNewMidValue < midValue) { // go right
                        low =x;
                        x = tentativeNewMid;
                        midValue = tentativeNewMidValue;
                    } else { // go left
                        high = tentativeNewMid;
                    }
                }
                ++iterations;
            }
            return x;
        }
    }

    SphereCylinderOptimizer::SphereCylinderOptimizer(Real r,
                                                     Real s,
                                                     Real alpha,
                                                     Real z1,
                                                     Real z2,
                                                     Real z3,
                                                     Real zweight)
    : r_(r), s_(s), alpha_(alpha), z1_(z1), z2_(z2), z3_(z3), zweight_(zweight)
    {

        QL_REQUIRE(r > 0, "sphere must have positive radius");

        s = std::max(s, 0.0);
        QL_REQUIRE(alpha > 0, "cylinder centre must have positive coordinate");

        nonEmpty_ = std::fabs(alpha - s) <= r;

        Real cylinderInside = r * r - (s + alpha) * (s + alpha);

        if (cylinderInside > 0.0) {
            topValue_ = alpha + s;
            bottomValue_ = alpha - s;
        } else {
            bottomValue_ = alpha - s;
            Real tmp = r * r - (s * s + alpha * alpha);

            if (tmp <= 0) { // max to left of maximum
                Real topValue2 = std::sqrt(s * s - tmp * tmp / (4 * alpha * alpha));
                topValue_ = alpha - std::sqrt(s * s - topValue2 * topValue2);
            } else {
                topValue_ = alpha + tmp / (2.0 * alpha);
            }
        }
    }

    bool SphereCylinderOptimizer::isIntersectionNonEmpty() const {
        return nonEmpty_;
    }

    void SphereCylinderOptimizer::findClosest(Size maxIterations,
                                              Real tolerance,
                                              Real& y1,
                                              Real& y2,
                                              Real& y3) const
    {
         Real x1,x2,x3;
         findByProjection(x1,x2,x3);

         y1 = BrentMinimize(
                bottomValue_, x1, topValue_,tolerance, maxIterations,
                [&](Real x){ return objectiveFunction(x); });
         y2 =std::sqrt(s_*s_ - (y1-alpha_)*(y1-alpha_));
         y3= std::sqrt(r_*r_ - y1*y1-y2*y2);
    }

    Real SphereCylinderOptimizer::objectiveFunction(Real x1) const
    {
   //     Real x1 = alpha_ - std::sqrt(s_*s_-x2*x2);

        Real x2sq = s_*s_ - (x1-alpha_)*(x1-alpha_);
         // a negative number will be minuscule and a result of rounding error
        Real x2 = x2sq >= 0.0 ? Real(std::sqrt(x2sq)) : 0.0;
        Real x3= std::sqrt(r_*r_ - x1*x1-x2*x2);

        Real err=0.0;
        err+= (x1-z1_)*(x1-z1_);
        err+= (x2-z2_)*(x2-z2_);
        err+= (x3-z3_)*(x3-z3_)*zweight_;

        return err;
    }

    bool SphereCylinderOptimizer::findByProjection(Real& y1,
                                                   Real& y2,
                                                   Real& y3) const {
        Real z1moved = z1_-alpha_;
        Real distance = std::sqrt( z1moved*z1moved + z2_*z2_);
        Real scale = s_/distance;
        Real y1moved = z1moved*scale;
        y1 = alpha_+ y1moved;
        y2 = scale*z2_;
        Real residual = r_*r_ - y1*y1 -y2*y2;
        if (residual >=0.0) {
            y3 = std::sqrt(residual);
            return true;
        }
        // we are outside the sphere
        if (!isIntersectionNonEmpty()) {
            y3=0.0;
            return false;
        }

       // intersection is non-empty but projection point is outside sphere
       // so take rightmost point
       y3 = 0.0;
       y1 = topValue_;
       y2 = std::sqrt(r_*r_ -y1*y1);

       return true;
    }

    std::vector<Real> sphereCylinderOptimizerClosest(Real r,
                                                     Real s,
                                                     Real alpha,
                                                     Real z1,
                                                     Real z2,
                                                     Real z3,
                                                     Natural maxIterations,
                                                     Real tolerance,
                                                     Real zweight)
    {

        SphereCylinderOptimizer optimizer(r, s, alpha, z1, z2, z3, zweight);
        std::vector<Real> y(3);

        QL_REQUIRE(optimizer.isIntersectionNonEmpty(),
                   "intersection empty so no solution");

        if (maxIterations ==0)
            optimizer.findByProjection(y[0], y[1], y[2]);
        else
            optimizer.findClosest(maxIterations, tolerance, y[0], y[1], y[2]);

        return y;
     }

}

1	/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3	/*
4	Copyright (C) 2007 Mark Joshi
5
6	This file is part of QuantLib, a free-software/open-source library
7	for financial quantitative analysts and developers - http://quantlib.org/
8
9	QuantLib is free software: you can redistribute it and/or modify it
10	under the terms of the QuantLib license. You should have received a
11	copy of the license along with this program; if not, please email
12	<quantlib-dev@lists.sf.net>. The license is also available online at
13	<https://www.quantlib.org/license.shtml>.
14
15	This program is distributed in the hope that it will be useful, but WITHOUT
16	ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17	FOR A PARTICULAR PURPOSE. See the license for more details.
18	*/
19
20	#include <ql/math/optimization/spherecylinder.hpp>
21	#include <ql/errors.hpp>
22	#include <algorithm>
23
24	namespace QuantLib {
25
26	namespace {
27
28	template<class F>
29	Real BrentMinimize(Real low,	114✔
30	Real mid,
31	Real high,
32	Real tolerance,
33	Size maxIt,
34	const F& objectiveFunction) {
35	Real W = 0.5*(3.0-std::sqrt(5.0));
36	Real x = Wlow+(1-W)high;	114✔
37	if (mid > low && mid < high)	114✔
38	x = mid;
39
40	Real midValue = objectiveFunction(x);
41
42	Size iterations = 0;
43	while (high-low > tolerance && iterations < maxIt) {	4,456✔
44	if (x - low > high -x) { // left interval is bigger	4,342✔
45	Real tentativeNewMid = Wlow+(1-W)x;	1,934✔
46	Real tentativeNewMidValue =
47	objectiveFunction(tentativeNewMid);
48
49	if (tentativeNewMidValue < midValue) { // go left	1,934✔
50	high =x;
51	x = tentativeNewMid;
52	midValue = tentativeNewMidValue;
53	} else { // go right
54	low = tentativeNewMid;
55	}
56	} else {
57	Real tentativeNewMid = Wx+(1-W)high;	2,408✔
58	Real tentativeNewMidValue =
59	objectiveFunction(tentativeNewMid);
60
61	if (tentativeNewMidValue < midValue) { // go right	2,408✔
62	low =x;
63	x = tentativeNewMid;
64	midValue = tentativeNewMidValue;
65	} else { // go left
66	high = tentativeNewMid;
67	}
68	}
69	++iterations;	4,342✔
70	}
71	return x;	114✔
72	}
73	}
74
75	SphereCylinderOptimizer::SphereCylinderOptimizer(Real r,	114✔
76	Real s,
77	Real alpha,
78	Real z1,
79	Real z2,
80	Real z3,
81	Real zweight)	114✔
82	: r_(r), s_(s), alpha_(alpha), z1_(z1), z2_(z2), z3_(z3), zweight_(zweight)	114✔
83	{
84
85	QL_REQUIRE(r > 0, "sphere must have positive radius");	114✔
86
87	s = std::max(s, 0.0);	114✔
88	QL_REQUIRE(alpha > 0, "cylinder centre must have positive coordinate");	114✔
89
90	nonEmpty_ = std::fabs(alpha - s) <= r;	114✔
91
92	Real cylinderInside = r * r - (s + alpha) * (s + alpha);	114✔
93
94	if (cylinderInside > 0.0) {	114✔
95	topValue_ = alpha + s;	15✔
96	bottomValue_ = alpha - s;	15✔
97	} else {
98	bottomValue_ = alpha - s;	99✔
99	Real tmp = r * r - (s * s + alpha * alpha);	99✔
100
101	if (tmp <= 0) { // max to left of maximum	99✔
102	Real topValue2 = std::sqrt(s * s - tmp * tmp / (4 * alpha * alpha));	85✔
103	topValue_ = alpha - std::sqrt(s * s - topValue2 * topValue2);	85✔
104	} else {
105	topValue_ = alpha + tmp / (2.0 * alpha);	14✔
106	}
107	}
108	}	114✔
109
110	bool SphereCylinderOptimizer::isIntersectionNonEmpty() const {	114✔
111	return nonEmpty_;	114✔
112	}
113
114	void SphereCylinderOptimizer::findClosest(Size maxIterations,	114✔
115	Real tolerance,
116	Real& y1,
117	Real& y2,
118	Real& y3) const
119	{
120	Real x1,x2,x3;
121	findByProjection(x1,x2,x3);	114✔
122
123	y1 = BrentMinimize(	228✔
124	bottomValue_, x1, topValue_,tolerance, maxIterations,	114✔
125	[&](Real x){ return objectiveFunction(x); });	4,570✔
126	y2 =std::sqrt(s_s_ - (y1-alpha_)(y1-alpha_));	114✔
127	y3= std::sqrt(r_r_ - y1y1-y2*y2);	114✔
128	}	114✔
129
130	Real SphereCylinderOptimizer::objectiveFunction(Real x1) const	4,456✔
131	{
132	// Real x1 = alpha_ - std::sqrt(s_s_-x2x2);
133
134	Real x2sq = s_s_ - (x1-alpha_)(x1-alpha_);	4,456✔
135	// a negative number will be minuscule and a result of rounding error
136	Real x2 = x2sq >= 0.0 ? Real(std::sqrt(x2sq)) : 0.0;	4,456✔
137	Real x3= std::sqrt(r_r_ - x1x1-x2*x2);	4,456✔
138
139	Real err=0.0;
140	err+= (x1-z1_)*(x1-z1_);	4,456✔
141	err+= (x2-z2_)*(x2-z2_);	4,456✔
142	err+= (x3-z3_)(x3-z3_)zweight_;	4,456✔
143
144	return err;	4,456✔
145	}
146
147	bool SphereCylinderOptimizer::findByProjection(Real& y1,	116✔
148	Real& y2,
149	Real& y3) const {
150	Real z1moved = z1_-alpha_;	116✔
151	Real distance = std::sqrt( z1movedz1moved + z2_z2_);	116✔
152	Real scale = s_/distance;	116✔
153	Real y1moved = z1moved*scale;	116✔
154	y1 = alpha_+ y1moved;	116✔
155	y2 = scale*z2_;	116✔
156	Real residual = r_r_ - y1y1 -y2*y2;	116✔
157	if (residual >=0.0) {	116✔
158	y3 = std::sqrt(residual);	114✔
159	return true;	114✔
160	}
161	// we are outside the sphere
162	if (!isIntersectionNonEmpty()) {	2✔
163	y3=0.0;	×
164	return false;	×
165	}
166
167	// intersection is non-empty but projection point is outside sphere
168	// so take rightmost point
169	y3 = 0.0;	2✔
170	y1 = topValue_;	2✔
171	y2 = std::sqrt(r_r_ -y1y1);	2✔
172
173	return true;	2✔
174	}
175
176	std::vector<Real> sphereCylinderOptimizerClosest(Real r,	×
177	Real s,
178	Real alpha,
179	Real z1,
180	Real z2,
181	Real z3,
182	Natural maxIterations,
183	Real tolerance,
184	Real zweight)
185	{
186
187	SphereCylinderOptimizer optimizer(r, s, alpha, z1, z2, z3, zweight);	×
188	std::vector<Real> y(3);	×
189
190	QL_REQUIRE(optimizer.isIntersectionNonEmpty(),	×
191	"intersection empty so no solution");
192
193	if (maxIterations ==0)	×
194	optimizer.findByProjection(y[0], y[1], y[2]);	×
195	else
196	optimizer.findClosest(maxIterations, tolerance, y[0], y[1], y[2]);	×
197
198	return y;	×
199	}	×
200
201	}

lballabio / QuantLib / 30255006850

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