QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
spherecylinder.cpp
Go to the documentation of this file.
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 <http://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
21#include <ql/errors.hpp>
22#include <algorithm>
23
24namespace QuantLib {
25
26 namespace {
27
28 template<class F>
29 Real BrentMinimize(Real low,
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 = W*low+(1-W)*high;
37 if (mid > low && mid < high)
38 x = mid;
39
40 Real midValue = objectiveFunction(x);
41
42 Size iterations = 0;
43 while (high-low > tolerance && iterations < maxIt) {
44 if (x - low > high -x) { // left interval is bigger
45 Real tentativeNewMid = W*low+(1-W)*x;
46 Real tentativeNewMidValue =
47 objectiveFunction(tentativeNewMid);
48
49 if (tentativeNewMidValue < midValue) { // go left
50 high =x;
51 x = tentativeNewMid;
52 midValue = tentativeNewMidValue;
53 } else { // go right
54 low = tentativeNewMid;
55 }
56 } else {
57 Real tentativeNewMid = W*x+(1-W)*high;
58 Real tentativeNewMidValue =
59 objectiveFunction(tentativeNewMid);
60
61 if (tentativeNewMidValue < midValue) { // go right
62 low =x;
63 x = tentativeNewMid;
64 midValue = tentativeNewMidValue;
65 } else { // go left
66 high = tentativeNewMid;
67 }
68 }
69 ++iterations;
70 }
71 return x;
72 }
73 }
74
76 Real s,
77 Real alpha,
78 Real z1,
79 Real z2,
80 Real z3,
81 Real zweight)
82 : r_(r), s_(s), alpha_(alpha), z1_(z1), z2_(z2), z3_(z3), zweight_(zweight)
83 {
84
85 QL_REQUIRE(r > 0, "sphere must have positive radius");
86
87 s = std::max(s, 0.0);
88 QL_REQUIRE(alpha > 0, "cylinder centre must have positive coordinate");
89
90 nonEmpty_ = std::fabs(alpha - s) <= r;
91
92 Real cylinderInside = r * r - (s + alpha) * (s + alpha);
93
94 if (cylinderInside > 0.0) {
95 topValue_ = alpha + s;
97 } else {
99 Real tmp = r * r - (s * s + alpha * alpha);
100
101 if (tmp <= 0) { // max to left of maximum
102 Real topValue2 = std::sqrt(s * s - tmp * tmp / (4 * alpha * alpha));
103 topValue_ = alpha - std::sqrt(s * s - topValue2 * topValue2);
104 } else {
105 topValue_ = alpha + tmp / (2.0 * alpha);
106 }
107 }
108 }
109
111 return nonEmpty_;
112 }
113
115 Real tolerance,
116 Real& y1,
117 Real& y2,
118 Real& y3) const
119 {
120 Real x1,x2,x3;
121 findByProjection(x1,x2,x3);
122
123 y1 = BrentMinimize(
124 bottomValue_, x1, topValue_,tolerance, maxIterations,
125 [&](Real x){ return objectiveFunction(x); });
126 y2 =std::sqrt(s_*s_ - (y1-alpha_)*(y1-alpha_));
127 y3= std::sqrt(r_*r_ - y1*y1-y2*y2);
128 }
129
131 {
132 // Real x1 = alpha_ - std::sqrt(s_*s_-x2*x2);
133
134 Real x2sq = s_*s_ - (x1-alpha_)*(x1-alpha_);
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;
137 Real x3= std::sqrt(r_*r_ - x1*x1-x2*x2);
138
139 Real err=0.0;
140 err+= (x1-z1_)*(x1-z1_);
141 err+= (x2-z2_)*(x2-z2_);
142 err+= (x3-z3_)*(x3-z3_)*zweight_;
143
144 return err;
145 }
146
148 Real& y2,
149 Real& y3) const {
150 Real z1moved = z1_-alpha_;
151 Real distance = std::sqrt( z1moved*z1moved + z2_*z2_);
152 Real scale = s_/distance;
153 Real y1moved = z1moved*scale;
154 y1 = alpha_+ y1moved;
155 y2 = scale*z2_;
156 Real residual = r_*r_ - y1*y1 -y2*y2;
157 if (residual >=0.0) {
158 y3 = std::sqrt(residual);
159 return true;
160 }
161 // we are outside the sphere
162 if (!isIntersectionNonEmpty()) {
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;
170 y1 = topValue_;
171 y2 = std::sqrt(r_*r_ -y1*y1);
172
173 return true;
174 }
175
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
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}
void findClosest(Size maxIterations, Real tolerance, Real &y1, Real &y2, Real &y3) const
SphereCylinderOptimizer(Real r, Real s, Real alpha, Real z1, Real z2, Real z3, Real zweight=1.0)
bool findByProjection(Real &y1, Real &y2, Real &y3) const
Classes and functions for error handling.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
QL_REAL Real
real number
Definition: types.hpp:50
unsigned QL_INTEGER Natural
positive integer
Definition: types.hpp:43
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
std::vector< Real > sphereCylinderOptimizerClosest(Real r, Real s, Real alpha, Real z1, Real z2, Real z3, Natural maxIterations, Real tolerance, Real zweight)
ext::shared_ptr< YieldTermStructure > r
const Real s_
Real F
Definition: sabr.cpp:200
Real alpha
Definition: sabr.cpp:200
Find closest point of the intersection of a sphere and cylinder to a given point.