Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2004 Sandia Corporation and Argonne National
5 : : Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 : : with Sandia Corporation, the U.S. Government retains certain
7 : : rights in this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : [email protected], [email protected], [email protected],
24 : : [email protected], [email protected], [email protected]
25 : :
26 : : ***************************************************************** */
27 : : /*!
28 : : \file SteepestDescent.cpp
29 : : \brief
30 : :
31 : : Implements the SteepestDescent class member functions.
32 : :
33 : : \author Thomas Leurent
34 : : \date 2002-06-13
35 : : */
36 : :
37 : : #include "SteepestDescent.hpp"
38 : : #include "MsqFreeVertexIndexIterator.hpp"
39 : : #include "MsqTimer.hpp"
40 : : #include "MsqDebug.hpp"
41 : :
42 : : #include <memory>
43 : :
44 : : namespace MBMesquite
45 : : {
46 : :
47 : 0 : std::string SteepestDescent::get_name() const
48 : : {
49 [ # # ]: 0 : return "SteepestDescent";
50 : : }
51 : :
52 : 36 : PatchSet* SteepestDescent::get_patch_set()
53 : : {
54 : 36 : return PatchSetUser::get_patch_set();
55 : : }
56 : :
57 : 33 : SteepestDescent::SteepestDescent( ObjectiveFunction* of )
58 [ + - ]: 33 : : VertexMover( of ), PatchSetUser( true ), projectGradient( false ) //,
59 : : // cosineStep(false)
60 : : {
61 : 33 : }
62 : :
63 : 42 : void SteepestDescent::initialize( PatchData& /*pd*/, MsqError& /*err*/ ) {}
64 : :
65 : 57762 : void SteepestDescent::initialize_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
66 : :
67 : 57762 : void SteepestDescent::optimize_vertex_positions( PatchData& pd, MsqError& err )
68 : : {
69 : : MSQ_FUNCTION_TIMER( "SteepestDescent::optimize_vertex_positions" );
70 : :
71 : 57762 : const int SEARCH_MAX = 100;
72 : 57762 : const double c1 = 1e-4;
73 : : // std::vector<Vector3D> unprojected(pd.num_free_vertices());
74 [ + - ][ + - ]: 57762 : std::vector< Vector3D > gradient( pd.num_free_vertices() );
75 : 57762 : bool feasible = true; // bool for OF values
76 : : double min_edge_len, max_edge_len;
77 : 57762 : double step_size = 0, original_value = 0, new_value = 0;
78 : 57762 : double norm_squared = 0;
79 : : PatchDataVerticesMemento* pd_previous_coords;
80 [ + - ]: 57762 : TerminationCriterion* term_crit = get_inner_termination_criterion();
81 [ + - ]: 57762 : OFEvaluator& obj_func = get_objective_function_evaluator();
82 : :
83 : : // get vertex memento so we can restore vertex coordinates for bad steps.
84 [ + - ][ + - ]: 57762 : pd_previous_coords = pd.create_vertices_memento( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
85 : : // use auto_ptr to automatically delete memento when we exit this function
86 [ + + ]: 115524 : std::auto_ptr< PatchDataVerticesMemento > memento_deleter( pd_previous_coords );
87 : :
88 : : // Evaluate objective function.
89 : : //
90 : : // Always use 'update' version when beginning optimization so that
91 : : // if doing block coordinate descent the OF code knows the set of
92 : : // vertices we are modifying during the optimziation (the subset
93 : : // of the mesh contained in the current patch.) This has to be
94 : : // done up-front because typically an OF will just store the portion
95 : : // of the OF value (e.g. the numeric contribution to the sum for an
96 : : // averaging OF) for the initial patch.
97 [ + - ][ + - ]: 57762 : feasible = obj_func.update( pd, original_value, gradient, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
98 : : // calculate gradient dotted with itself
99 [ + - ]: 57762 : norm_squared = length_squared( gradient );
100 : :
101 : : // set an error if initial patch is invalid.
102 [ - + ]: 57762 : if( !feasible )
103 : : {
104 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "SteepestDescent passed invalid initial patch.", MsqError::INVALID_ARG );
105 : 0 : return;
106 : : }
107 : :
108 : : // use edge length as an initial guess for for step size
109 [ + - ]: 57762 : pd.get_minmax_edge_length( min_edge_len, max_edge_len );
110 : : // step_size = max_edge_len / std::sqrt(norm_squared);
111 : : // if (!finite(step_size)) // zero-length gradient
112 : : // return;
113 : : // if (norm_squared < DBL_EPSILON)
114 : : // return;
115 [ + + ][ + - ]: 57762 : if( norm_squared >= DBL_EPSILON ) step_size = max_edge_len / std::sqrt( norm_squared ) * pd.num_free_vertices();
116 : :
117 : : // The steepest descent loop...
118 : : // We loop until the user-specified termination criteria are met.
119 [ + - ][ + + ]: 116322 : while( !term_crit->terminate() )
[ + + ]
120 : : {
121 : : MSQ_DBGOUT( 3 ) << "Iteration " << term_crit->get_iteration_count() << std::endl;
122 : : MSQ_DBGOUT( 3 ) << " o original_value: " << original_value << std::endl;
123 : : MSQ_DBGOUT( 3 ) << " o grad norm suqared: " << norm_squared << std::endl;
124 : :
125 : : // Save current vertex coords so that they can be restored if
126 : : // the step was bad.
127 [ + - ][ + - ]: 58560 : pd.recreate_vertices_memento( pd_previous_coords, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
128 : :
129 : : // Reduce step size until it satisfies Armijo condition
130 : 58560 : int counter = 0;
131 : : for( ;; )
132 : : {
133 [ + - ][ + + ]: 234472 : if( ++counter > SEARCH_MAX || step_size < DBL_EPSILON )
[ + + ]
134 : : {
135 : : MSQ_DBGOUT( 3 ) << " o No valid step found. Giving Up." << std::endl;
136 : 20271 : return;
137 : : }
138 : :
139 : : // Move vertices to new positions.
140 : : // Note: step direction is -gradient so we pass +gradient and
141 : : // -step_size to achieve the same thing.
142 [ + - ][ + - ]: 214201 : pd.move_free_vertices_constrained( arrptr( gradient ), gradient.size(), -step_size, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
143 : : // Evaluate objective function for new vertices. We call the
144 : : // 'evaluate' form here because we aren't sure yet if we want to
145 : : // keep these vertices. Until we call 'update', we have the option
146 : : // of reverting a block coordinate decent objective function's state
147 : : // to that of the initial vertex coordinates. However, for block
148 : : // coordinate decent to work correctly, we will need to call an
149 : : // 'update' form if we decide to keep the new vertex coordinates.
150 [ + - ]: 214201 : feasible = obj_func.evaluate( pd, new_value, err );
151 [ + - ][ + + ]: 214201 : if( err.error_code() == err.BARRIER_VIOLATED )
152 [ + - ]: 24 : err.clear(); // barrier violated does not represent an actual error here
153 [ + - ][ - + ]: 214201 : MSQ_ERRRTN( err );
[ # # ][ # # ]
[ - + ]
154 : : MSQ_DBGOUT( 3 ) << " o step_size: " << step_size << std::endl;
155 : : MSQ_DBGOUT( 3 ) << " o new_value: " << new_value << std::endl;
156 : :
157 [ + + ]: 214201 : if( !feasible )
158 : : {
159 : : // OF value is invalid, decrease step_size a lot
160 : 25 : step_size *= 0.2;
161 : : }
162 [ + + ]: 214176 : else if( new_value > original_value - c1 * step_size * norm_squared )
163 : : {
164 : : // Armijo condition not met.
165 : 175887 : step_size *= 0.5;
166 : : }
167 : : else
168 : : {
169 : : // Armijo condition met, stop
170 : 38289 : break;
171 : : }
172 : :
173 : : // undo previous step : restore vertex coordinates
174 [ + - ][ + - ]: 175912 : pd.set_to_vertices_memento( pd_previous_coords, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
175 : : }
176 : :
177 : : // Re-evaluate objective function to get gradient.
178 : : // Calling the 'update' form here incorporates the new vertex
179 : : // positions into the 'accumulated' value if we are doing a
180 : : // block coordinate descent optimization.
181 [ + - ][ + - ]: 38289 : obj_func.update( pd, original_value, gradient, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
182 [ - + ]: 38289 : if( projectGradient )
183 : : {
184 : : // if (cosineStep) {
185 : : // unprojected = gradient;
186 : : // pd.project_gradient( gradient, err ); MSQ_ERRRTN(err);
187 : : // double dot = inner_product( arrptr(gradient), arrptr(unprojected), gradient.size()
188 : : // ); double lensqr1 = length_squared( gradient ); double lensqr2 = length_squared(
189 : : // unprojected ); double cossqr = dot * dot / lensqr1 / lensqr2; step_size *=
190 : : // sqrt(cossqr);
191 : : //}
192 : : // else {
193 [ # # ][ # # ]: 0 : pd.project_gradient( gradient, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
194 : : //}
195 : : }
196 : :
197 : : // Update terination criterion for next iteration.
198 : : // This is necessary for efficiency. Some values can be adjusted
199 : : // for each iteration so we don't need to re-caculate the value
200 : : // over the entire mesh.
201 [ + - ][ + - ]: 38289 : term_crit->accumulate_patch( pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
202 [ + - ][ + - ]: 38289 : term_crit->accumulate_inner( pd, original_value, arrptr( gradient ), err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
203 : :
204 : : // Calculate initial step size for next iteration using step size
205 : : // from this iteration
206 : 38289 : step_size *= norm_squared;
207 [ + - ]: 38289 : norm_squared = length_squared( gradient );
208 : : // if (norm_squared < DBL_EPSILON)
209 : : // break;
210 [ + + ]: 38289 : if( norm_squared >= DBL_EPSILON ) step_size /= norm_squared;
211 : 233674 : }
212 : : }
213 : :
214 : 342 : void SteepestDescent::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ )
215 : : {
216 : : // cout << "- Executing SteepestDescent::iteration_complete()\n";
217 : 342 : }
218 : :
219 : 36 : void SteepestDescent::cleanup()
220 : : {
221 : : // cout << "- Executing SteepestDescent::iteration_end()\n";
222 : 36 : }
223 : :
224 [ + - ][ + - ]: 36 : } // namespace MBMesquite
|
