1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273
|
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkPowellOptimizerv4_h
#define itkPowellOptimizerv4_h
#include "itkVector.h"
#include "itkMatrix.h"
#include "itkObjectToObjectOptimizerBase.h"
namespace itk
{
/**
* \class PowellOptimizerv4
* \brief Implements Powell optimization using Brent line search.
*
* The code in this class was adapted from the Wikipedia and the
* netlib.org zeroin function.
*
* https://www.netlib.org/go/zeroin.f
* https://en.wikipedia.org/wiki/Brent_method
* https://en.wikipedia.org/wiki/Golden_section_search
*
* This optimizer needs a cost function.
* Partial derivatives of that function are not required.
*
* For an N-dimensional parameter space, each iteration minimizes
* the function in N (initially orthogonal) directions. Typically only 2-5
* iterations are required. If gradients are available, consider a conjugate
* gradient line search strategy.
*
* The SetStepLength determines the initial distance to step in a line direction
* when bounding the minimum (using bracketing triple spaced using a golden
* search strategy).
*
* The StepTolerance terminates optimization when the parameter values are
* known to be within this (scaled) distance of the local extreme.
*
* The ValueTolerance terminates optimization when the cost function values at
* the current parameters and at the local extreme are likely (within a second
* order approximation) to be within this is tolerance.
*
* \ingroup ITKOptimizersv4
*/
template <typename TInternalComputationValueType>
class ITK_TEMPLATE_EXPORT PowellOptimizerv4 : public ObjectToObjectOptimizerBaseTemplate<TInternalComputationValueType>
{
public:
/** Standard "Self" type alias. */
using Self = PowellOptimizerv4;
using Superclass = ObjectToObjectOptimizerBaseTemplate<TInternalComputationValueType>;
using Pointer = SmartPointer<Self>;
using ConstPointer = SmartPointer<const Self>;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** \see LightObject::GetNameOfClass() */
itkOverrideGetNameOfClassMacro(PowellOptimizerv4);
using typename Superclass::ParametersType;
using typename Superclass::MeasureType;
using typename Superclass::ScalesType;
/** Set/Get maximum iteration limit. */
itkSetMacro(MaximumIteration, unsigned int);
itkGetConstReferenceMacro(MaximumIteration, unsigned int);
/** Set/Get the maximum number of line search iterations */
itkSetMacro(MaximumLineIteration, unsigned int);
itkGetConstMacro(MaximumLineIteration, unsigned int);
/** Set/Get StepLength for the (scaled) spacing of the sampling of
* parameter space while bracketing the extremum */
itkSetMacro(StepLength, double);
itkGetConstReferenceMacro(StepLength, double);
/** Set/Get StepTolerance. Once the local extreme is known to be within this
* distance of the current parameter values, optimization terminates */
itkSetMacro(StepTolerance, double);
itkGetConstReferenceMacro(StepTolerance, double);
/** Set/Get ValueTolerance. Once this current cost function value is known
* to be within this tolerance of the cost function value at the local
* extreme, optimization terminates */
itkSetMacro(ValueTolerance, double);
itkGetConstReferenceMacro(ValueTolerance, double);
/** Return Current Value */
itkGetConstReferenceMacro(CurrentCost, MeasureType);
const MeasureType &
GetValue() const override
{
return this->GetCurrentCost();
}
/** Get the current line search iteration */
itkGetConstReferenceMacro(CurrentLineIteration, unsigned int);
/** Start optimization. */
void
StartOptimization(bool doOnlyInitialization = false) override;
/** When users call StartOptimization, this value will be set false.
* By calling StopOptimization, this flag will be set true, and
* optimization will stop at the next iteration. */
void
StopOptimization()
{
m_Stop = true;
}
itkGetConstReferenceMacro(CatchGetValueException, bool);
itkSetMacro(CatchGetValueException, bool);
itkBooleanMacro(CatchGetValueException);
itkGetConstReferenceMacro(MetricWorstPossibleValue, double);
itkSetMacro(MetricWorstPossibleValue, double);
const std::string
GetStopConditionDescription() const override;
protected:
PowellOptimizerv4();
PowellOptimizerv4(const PowellOptimizerv4 &);
~PowellOptimizerv4() override = default;
void
PrintSelf(std::ostream & os, Indent indent) const override;
itkSetMacro(CurrentCost, double);
/** Used to specify the line direction through the n-dimensional parameter
* space the is currently being bracketed and optimized. */
void
SetLine(const ParametersType & origin, const vnl_vector<double> & direction);
/** Get the value of the n-dimensional cost function at this scalar step
* distance along the current line direction from the current line origin.
* Line origin and distances are set via SetLine */
double
GetLineValue(double x) const;
double
GetLineValue(double x, ParametersType & tempCoord) const;
/** Set the given scalar step distance (x) and function value (fx) as the
* "best-so-far" optimizer values. */
void
SetCurrentLinePoint(double x, double fx);
/** Used in bracketing the extreme along the current line.
* Adapted from NRC */
void
Swap(double * a, double * b) const;
/** Used in bracketing the extreme along the current line.
* Adapted from NRC */
void
Shift(double * a, double * b, double * c, double d) const;
/** The LineBracket routine from NRC. Later reimplemented from the description
* of the method available in the Wikipedia.
*
* Uses current origin and line direction (from SetLine) to find a triple of
* points (ax, bx, cx) that bracket the extreme "near" the origin. Search
* first considers the point StepLength distance from ax.
*
* IMPORTANT: The value of ax and the value of the function at ax (i.e., fa),
* must both be provided to this function. */
virtual void
LineBracket(double * x1, double * x2, double * x3, double * f1, double * f2, double * f3);
virtual void
LineBracket(double * x1, double * x2, double * x3, double * f1, double * f2, double * f3, ParametersType & tempCoord);
/** Given a bracketing triple of points and their function values, returns
* a bounded extreme. These values are in parameter space, along the
* current line and wrt the current origin set via SetLine. Optimization
* terminates based on MaximumIteration, StepTolerance, or ValueTolerance.
* Implemented as Brent line optimizers from NRC. */
virtual void
BracketedLineOptimize(double ax,
double bx,
double cx,
double fa,
double functionValueOfb,
double fc,
double * extX,
double * extVal);
virtual void
BracketedLineOptimize(double ax,
double bx,
double cx,
double fa,
double functionValueOfb,
double fc,
double * extX,
double * extVal,
ParametersType & tempCoord);
itkGetMacro(SpaceDimension, unsigned int);
void
SetSpaceDimension(unsigned int dim)
{
this->m_SpaceDimension = dim;
this->m_LineDirection.set_size(dim);
this->m_LineOrigin.set_size(dim);
this->m_CurrentPosition.set_size(dim);
this->Modified();
}
itkSetMacro(CurrentIteration, unsigned int);
itkGetMacro(Stop, bool);
itkSetMacro(Stop, bool);
private:
unsigned int m_SpaceDimension{ 0 };
/** Current iteration */
unsigned int m_CurrentLineIteration{ 0 };
/** Maximum iteration limit. */
unsigned int m_MaximumIteration{ 100 };
unsigned int m_MaximumLineIteration{ 100 };
bool m_CatchGetValueException{ false };
double m_MetricWorstPossibleValue{ 0 };
/** The minimal size of search */
double m_StepLength{ 0 };
double m_StepTolerance{ 0 };
ParametersType m_LineOrigin{};
vnl_vector<double> m_LineDirection{};
double m_ValueTolerance{ 0 };
/** Internal storage for the value type / used as a cache */
MeasureType m_CurrentCost{};
/** this is user-settable flag to stop optimization.
* when users call StartOptimization, this value will be set false.
* By calling StopOptimization, this flag will be set true, and
* optimization will stop at the next iteration. */
bool m_Stop{ false };
ParametersType m_CurrentPosition{};
std::ostringstream m_StopConditionDescription{};
}; // end of class
} // end of namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
# include "itkPowellOptimizerv4.hxx"
#endif
#endif
|