1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18 #ifndef itkOnePlusOneEvolutionaryOptimizerv4_hxx
19 #define itkOnePlusOneEvolutionaryOptimizerv4_hxx
20
21 #include "itkMath.h"
22 #include "itkOnePlusOneEvolutionaryOptimizerv4.h"
23 #include "vnl/vnl_matrix.h"
24 namespace itk
25 {
26 template<typename TInternalComputationValueType>
27 OnePlusOneEvolutionaryOptimizerv4<TInternalComputationValueType>
OnePlusOneEvolutionaryOptimizerv4()28 ::OnePlusOneEvolutionaryOptimizerv4()
29 {
30 m_CatchGetValueException = false;
31 m_MetricWorstPossibleValue = 0;
32
33 m_Epsilon = (double)1.5e-4;
34 m_RandomGenerator = nullptr;
35
36 m_Initialized = false;
37 m_GrowthFactor = 1.05;
38 m_ShrinkFactor = std::pow(m_GrowthFactor, -0.25);
39 m_InitialRadius = 1.01;
40 m_MaximumIteration = 100;
41 m_Stop = false;
42 m_StopConditionDescription.str("");
43 m_CurrentCost = 0;
44 m_FrobeniusNorm = 0.0;
45 }
46
47 template<typename TInternalComputationValueType>
48 void
49 OnePlusOneEvolutionaryOptimizerv4<TInternalComputationValueType>
SetNormalVariateGenerator(NormalVariateGeneratorType * generator)50 ::SetNormalVariateGenerator(NormalVariateGeneratorType *generator)
51 {
52 if ( m_RandomGenerator != generator )
53 {
54 m_RandomGenerator = generator;
55 this->Modified();
56 }
57 }
58
59 template<typename TInternalComputationValueType>
60 void
61 OnePlusOneEvolutionaryOptimizerv4<TInternalComputationValueType>
Initialize(double initialRadius,double grow,double shrink)62 ::Initialize(double initialRadius, double grow, double shrink)
63 {
64 m_InitialRadius = initialRadius;
65
66 if ( Math::AlmostEquals( grow, -1 ) )
67 {
68 m_GrowthFactor = 1.05;
69 }
70 else
71 {
72 m_GrowthFactor = grow;
73 }
74 if ( Math::AlmostEquals( shrink, -1 ) )
75 {
76 m_ShrinkFactor = std::pow(m_GrowthFactor, -0.25);
77 }
78 else
79 {
80 m_ShrinkFactor = shrink;
81 }
82 }
83
84 template<typename TInternalComputationValueType>
85 void
86 OnePlusOneEvolutionaryOptimizerv4<TInternalComputationValueType>
StartOptimization(bool)87 ::StartOptimization(bool /* doOnlyInitialization */)
88 {
89 if ( this->m_Metric.IsNull() )
90 {
91 return;
92 }
93
94 Superclass::StartOptimization();
95
96 this->InvokeEvent( StartEvent() );
97 m_Stop = false;
98
99 unsigned int spaceDimension = this->m_Metric->GetNumberOfParameters();
100 vnl_matrix< double > A(spaceDimension, spaceDimension);
101 vnl_vector< double > parent( this->m_Metric->GetParameters() );
102 vnl_vector< double > f_norm(spaceDimension);
103 vnl_vector< double > child(spaceDimension);
104 vnl_vector< double > delta(spaceDimension);
105
106 ParametersType parentPosition(spaceDimension);
107 ParametersType childPosition(spaceDimension);
108
109 for ( unsigned int i = 0; i < spaceDimension; i++ )
110 {
111 parentPosition[i] = parent[i];
112 }
113 this->m_Metric->SetParameters( parentPosition );
114
115 double pvalue = m_MetricWorstPossibleValue;
116 try
117 {
118 pvalue = this->m_Metric->GetValue();
119 }
120 catch ( ... )
121 {
122 if ( m_CatchGetValueException )
123 {
124 pvalue = m_MetricWorstPossibleValue;
125 }
126 else
127 {
128 throw;
129 }
130 }
131
132 itkDebugMacro(<< ": initial position: " << parentPosition);
133 itkDebugMacro(<< ": initial fitness: " << pvalue);
134
135 this->m_Metric->SetParameters(parentPosition);
136 const ScalesType & scales = this->GetScales();
137
138 // Make sure the scales have been set properly
139 if ( scales.size() != spaceDimension )
140 {
141 itkExceptionMacro(<< "The size of Scales is "
142 << scales.size()
143 << ", but the NumberOfParameters for the CostFunction is "
144 << spaceDimension
145 << ".");
146 }
147
148 A.set_identity();
149 for ( unsigned int i = 0; i < spaceDimension; i++ )
150 {
151 A(i, i) = m_InitialRadius / scales[i];
152 }
153
154 for ( this->m_CurrentIteration = 0;
155 this->m_CurrentIteration < m_MaximumIteration;
156 this->m_CurrentIteration++ )
157 {
158 if ( m_Stop )
159 {
160 m_StopConditionDescription.str("");
161 m_StopConditionDescription << this->GetNameOfClass() << ": ";
162 m_StopConditionDescription << "StopOptimization() called";
163 break;
164 }
165
166 for ( unsigned int i = 0; i < spaceDimension; i++ )
167 {
168 if ( !m_RandomGenerator )
169 {
170 itkExceptionMacro(<< "Random Generator is not set!");
171 }
172 f_norm[i] = m_RandomGenerator->GetVariate();
173 }
174
175 delta = A * f_norm;
176 child = parent + delta;
177
178 for ( unsigned int i = 0; i < spaceDimension; i++ )
179 {
180 childPosition[i] = child[i];
181 }
182 // Update the metric so we can check the metric value in childPosition
183 this->m_Metric->SetParameters( childPosition );
184
185 double cvalue = m_MetricWorstPossibleValue;
186 try
187 {
188 cvalue = this->m_Metric->GetValue();
189 // While we got the metric value in childPosition,
190 // the metric parameteres are set back to parentPosition
191 this->m_Metric->SetParameters( parentPosition );
192 }
193 catch ( ... )
194 {
195 if ( m_CatchGetValueException )
196 {
197 cvalue = m_MetricWorstPossibleValue;
198 }
199 else
200 {
201 throw;
202 }
203 }
204
205 itkDebugMacro(<< "iter: " << this->m_CurrentIteration << ": parent position: "
206 << parentPosition);
207 itkDebugMacro(<< "iter: " << this->m_CurrentIteration << ": parent fitness: "
208 << pvalue);
209 itkDebugMacro(<< "iter: " << this->m_CurrentIteration << ": random vector: " << f_norm);
210 itkDebugMacro(<< "iter: " << this->m_CurrentIteration << ": A: " << std::endl << A);
211 itkDebugMacro(<< "iter: " << this->m_CurrentIteration << ": delta: " << delta);
212 itkDebugMacro(<< "iter: " << this->m_CurrentIteration << ": child position: "
213 << childPosition);
214 itkDebugMacro(<< "iter: " << this->m_CurrentIteration << ": child fitness: "
215 << cvalue);
216
217 double adjust = m_ShrinkFactor;
218
219 if ( cvalue < pvalue )
220 {
221 itkDebugMacro(<< "iter: " << this->m_CurrentIteration << ": increasing search radius");
222 pvalue = cvalue;
223 parent.swap(child);
224 adjust = m_GrowthFactor;
225 for ( unsigned int i = 0; i < spaceDimension; i++ )
226 {
227 parentPosition[i] = parent[i];
228 }
229 this->m_Metric->SetParameters(parentPosition);
230 }
231 else
232 {
233 itkDebugMacro(<< "iter: " << this->m_CurrentIteration << ": decreasing search radius");
234 }
235
236 m_CurrentCost = pvalue;
237 // convergence criterion: f-norm of A < epsilon_A
238 // Compute double precision sum of absolute values of
239 // a single precision vector
240 m_FrobeniusNorm = A.fro_norm();
241 itkDebugMacro(<< "A f-norm:" << m_FrobeniusNorm);
242 if ( m_FrobeniusNorm <= m_Epsilon )
243 {
244 itkDebugMacro(<< "converges at iteration = " << this->m_CurrentIteration);
245 m_StopConditionDescription.str("");
246 m_StopConditionDescription << this->GetNameOfClass() << ": ";
247 m_StopConditionDescription << "Fnorm (" << m_FrobeniusNorm
248 << ") is less than Epsilon (" << m_Epsilon
249 << " at iteration #" << this->m_CurrentIteration;
250 this->InvokeEvent( EndEvent() );
251 return;
252 }
253
254 // A += (adjust - 1)/ (f_norm * f_norm) * A * f_norm * f_norm;
255 // Blas_R1_Update(A, A * f_norm, f_norm,
256 // ((adjust - 1) / Blas_Dot_Prod(f_norm, f_norm)));
257 // = DGER(Fortran)
258 // performs the rank 1 operation
259 // A := alpha*x*y' + A,
260 // where y' = transpose(y)
261 // where alpha is a scalar, x is an m element vector, y is an n element
262 // vector and A is an m by n matrix.
263 // x = A * f_norm , y = f_norm, alpha = (adjust - 1) / Blas_Dot_Prod(
264 // f_norm, f_norm)
265
266 //A = A + (adjust - 1.0) * A;
267 double alpha = ( ( adjust - 1.0 ) / dot_product(f_norm, f_norm) );
268 for ( unsigned int c = 0; c < spaceDimension; c++ )
269 {
270 for ( unsigned int r = 0; r < spaceDimension; r++ )
271 {
272 A(r, c) += alpha * delta[r] * f_norm[c];
273 }
274 }
275
276 this->InvokeEvent( IterationEvent() );
277 itkDebugMacro( << "Current position: " << this->GetCurrentPosition() );
278 }
279 m_StopConditionDescription.str("");
280 m_StopConditionDescription << this->GetNameOfClass() << ": ";
281 m_StopConditionDescription << "Maximum number of iterations ("
282 << m_MaximumIteration
283 << ") exceeded. ";
284 this->InvokeEvent( EndEvent() );
285 }
286
287 template<typename TInternalComputationValueType>
288 const std::string
289 OnePlusOneEvolutionaryOptimizerv4<TInternalComputationValueType>
GetStopConditionDescription() const290 ::GetStopConditionDescription() const
291 {
292 return m_StopConditionDescription.str();
293 }
294
295 template<typename TInternalComputationValueType>
296 const typename OnePlusOneEvolutionaryOptimizerv4<TInternalComputationValueType>::MeasureType &
297 OnePlusOneEvolutionaryOptimizerv4<TInternalComputationValueType>
GetValue() const298 ::GetValue() const
299 {
300 return this->GetCurrentCost();
301 }
302
303 template<typename TInternalComputationValueType>
304 void
305 OnePlusOneEvolutionaryOptimizerv4<TInternalComputationValueType>
PrintSelf(std::ostream & os,Indent indent) const306 ::PrintSelf(std::ostream & os, Indent indent) const
307 {
308 Superclass::PrintSelf(os, indent);
309
310 if ( m_RandomGenerator )
311 {
312 os << indent << "Random Generator " << m_RandomGenerator.GetPointer() << std::endl;
313 }
314 else
315 {
316 os << indent << "Random Generator " << "(none)" << std::endl;
317 }
318 os << indent << "Maximum Iteration " << GetMaximumIteration() << std::endl;
319 os << indent << "Epsilon " << GetEpsilon() << std::endl;
320 os << indent << "Initial Radius " << GetInitialRadius() << std::endl;
321 os << indent << "Growth Fractor " << GetGrowthFactor() << std::endl;
322 os << indent << "Shrink Fractor " << GetShrinkFactor() << std::endl;
323 os << indent << "Initialized " << GetInitialized() << std::endl;
324 os << indent << "Current Cost " << GetCurrentCost() << std::endl;
325 os << indent << "Frobenius Norm " << GetFrobeniusNorm() << std::endl;
326 os << indent << "CatchGetValueException " << GetCatchGetValueException()
327 << std::endl;
328 os << indent << "MetricWorstPossibleValue " << GetMetricWorstPossibleValue()
329 << std::endl;
330 }
331 } // end of namespace itk
332 #endif
333