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 #include "itkExhaustiveOptimizer.h"
19 
20 namespace itk
21 {
22 /**
23  * Constructor
24  */
25 ExhaustiveOptimizer
ExhaustiveOptimizer()26 ::ExhaustiveOptimizer()
27 {
28   itkDebugMacro("Constructor");
29 
30   m_StepLength = 1.0;
31   m_CurrentIteration   =   0;
32   m_CurrentValue = 0;
33   m_CurrentParameter = 0;
34   m_CurrentIndex.Fill(0);
35   m_Stop = false;
36   m_NumberOfSteps.Fill(0);
37 
38   m_StopConditionDescription.str("");
39 }
40 
41 /**
42  * Start walking
43  */
44 
StartOptimization()45 void ExhaustiveOptimizer::StartOptimization()
46 {
47   this->StartWalking();
48 }
49 
50 void
51 ExhaustiveOptimizer
StartWalking()52 ::StartWalking()
53 {
54   itkDebugMacro("StartWalking");
55   this->InvokeEvent( StartEvent() );
56   m_StopConditionDescription.str("");
57   m_StopConditionDescription << this->GetNameOfClass() << ": Running";
58 
59   ParametersType initialPos = this->GetInitialPosition();
60   m_MinimumMetricValuePosition = initialPos;
61   m_MaximumMetricValuePosition = initialPos;
62 
63   MeasureType initialValue = this->GetValue( this->GetInitialPosition() );
64   m_MaximumMetricValue = initialValue;
65   m_MinimumMetricValue = initialValue;
66 
67   m_CurrentIteration          = 0;
68   m_MaximumNumberOfIterations = 1;
69 
70   const unsigned int spaceDimension = this->GetInitialPosition().GetSize();
71 
72   for ( unsigned int i = 0; i < spaceDimension; i++ )
73     {
74     m_MaximumNumberOfIterations *= ( 2 * m_NumberOfSteps[i] + 1 );
75     }
76 
77   m_CurrentIndex.SetSize(spaceDimension);
78   m_CurrentIndex.Fill(0);
79 
80   const ScalesType & scales = this->GetScales();
81   // Make sure the scales have been set properly
82   if ( scales.size() != spaceDimension )
83     {
84     itkExceptionMacro(<< "The size of Scales is "
85                       << scales.size()
86                       << ", but the NumberOfParameters is "
87                       << spaceDimension
88                       << ".");
89     }
90 
91   // Setup first grid position.
92   ParametersType position(spaceDimension);
93   for ( unsigned int i = 0; i < spaceDimension; i++ )
94     {
95     position[i] = this->GetInitialPosition()[i] - m_NumberOfSteps[i] * m_StepLength * scales[i];
96     }
97   this->SetCurrentPosition(position);
98 
99   itkDebugMacro("Calling ResumeWalking");
100 
101   this->ResumeWalking();
102 }
103 
104 /**
105  * Resume the optimization
106  */
107 void
108 ExhaustiveOptimizer
ResumeWalking()109 ::ResumeWalking()
110 {
111   itkDebugMacro("ResumeWalk");
112   m_Stop = false;
113 
114   while ( !m_Stop )
115     {
116     ParametersType currentPosition = this->GetCurrentPosition();
117 
118     if ( m_Stop )
119       {
120       StopWalking();
121       break;
122       }
123 
124     m_CurrentValue = this->GetValue(currentPosition);
125 
126     if ( m_CurrentValue > m_MaximumMetricValue )
127       {
128       m_MaximumMetricValue = m_CurrentValue;
129       m_MaximumMetricValuePosition = currentPosition;
130       }
131     if ( m_CurrentValue < m_MinimumMetricValue )
132       {
133       m_MinimumMetricValue = m_CurrentValue;
134       m_MinimumMetricValuePosition = currentPosition;
135       }
136 
137     if ( m_Stop )
138       {
139       this->StopWalking();
140       break;
141       }
142 
143     m_StopConditionDescription.str("");
144     m_StopConditionDescription << this->GetNameOfClass() << ": Running. ";
145     m_StopConditionDescription << "@ index " << this->GetCurrentIndex() << " value is " << this->GetCurrentValue();
146 
147     this->InvokeEvent( IterationEvent() );
148     this->AdvanceOneStep();
149     m_CurrentIteration++;
150     }
151 }
152 
153 void
154 ExhaustiveOptimizer
StopWalking()155 ::StopWalking()
156 {
157   itkDebugMacro("StopWalking");
158 
159   m_Stop = true;
160   this->InvokeEvent( EndEvent() );
161 }
162 
163 void
164 ExhaustiveOptimizer
AdvanceOneStep()165 ::AdvanceOneStep()
166 {
167   itkDebugMacro("AdvanceOneStep");
168 
169   const unsigned int spaceDimension = m_CostFunction->GetNumberOfParameters();
170 
171   ParametersType newPosition(spaceDimension);
172   IncrementIndex(newPosition);
173 
174   itkDebugMacro(<< "new position = " << newPosition);
175 
176   this->SetCurrentPosition(newPosition);
177 }
178 
179 void
180 ExhaustiveOptimizer
IncrementIndex(ParametersType & newPosition)181 ::IncrementIndex(ParametersType & newPosition)
182 {
183   unsigned int       idx = 0;
184   const unsigned int spaceDimension = m_CostFunction->GetNumberOfParameters();
185 
186   while ( idx < spaceDimension )
187     {
188     m_CurrentIndex[idx]++;
189 
190     if ( m_CurrentIndex[idx] > ( 2 * m_NumberOfSteps[idx] ) )
191       {
192       m_CurrentIndex[idx] = 0;
193       idx++;
194       }
195     else
196       {
197       break;
198       }
199     }
200 
201   if ( idx == spaceDimension )
202     {
203     m_Stop = true;
204     m_StopConditionDescription.str("");
205     m_StopConditionDescription << this->GetNameOfClass() << ": ";
206     m_StopConditionDescription << "Completed sampling of parametric space of size " << spaceDimension;
207     }
208 
209   const ScalesType & scales = this->GetScales();
210   for ( unsigned int i = 0; i < spaceDimension; i++ )
211     {
212     newPosition[i] = ( m_CurrentIndex[i] - m_NumberOfSteps[i] )
213                      * m_StepLength * scales[i]
214                      + this->GetInitialPosition()[i];
215     }
216 }
217 
218 const std::string
219 ExhaustiveOptimizer
GetStopConditionDescription() const220 ::GetStopConditionDescription() const
221 {
222   return m_StopConditionDescription.str();
223 }
224 
225 void
226 ExhaustiveOptimizer
PrintSelf(std::ostream & os,Indent indent) const227 ::PrintSelf(std::ostream & os, Indent indent) const
228 {
229   Superclass::PrintSelf(os, indent);
230 
231   os << indent << "CurrentValue = " << m_CurrentValue << std::endl;
232   os << indent << "NumberOfSteps = " << m_NumberOfSteps << std::endl;
233   os << indent << "CurrentIteration = " << m_CurrentIteration << std::endl;
234   os << indent << "Stop = " << m_Stop << std::endl;
235   os << indent << "CurrentParameter = " << m_CurrentParameter << std::endl;
236   os << indent << "StepLength = " << m_StepLength << std::endl;
237   os << indent << "CurrentIndex = " << m_CurrentIndex << std::endl;
238   os << indent << "MaximumNumberOfIterations = " << m_MaximumNumberOfIterations << std::endl;
239   os << indent << "MaximumMetricValue = " << m_MaximumMetricValue << std::endl;
240   os << indent << "MinimumMetricValue = " << m_MinimumMetricValue << std::endl;
241   os << indent << "MinimumMetricValuePosition = " << m_MinimumMetricValuePosition << std::endl;
242   os << indent << "MaximumMetricValuePosition = " << m_MaximumMetricValuePosition << std::endl;
243 }
244 } // end namespace itk
245