Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2006 Sandia National Laboratories. Developed at the
5 : : University of Wisconsin--Madison under SNL contract number
6 : : 624796. The U.S. Government and the University of Wisconsin
7 : : retain certain rights to 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 : : (2006) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file PMeanPMetric.cpp
28 : : * \brief
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "PMeanPMetric.hpp"
34 : : #include "MsqError.hpp"
35 : : #include "QualityMetric.hpp"
36 : : #include "Vector3D.hpp"
37 : : #include "Matrix3D.hpp"
38 : : #include "SymMatrix3D.hpp"
39 : : #include "PatchData.hpp"
40 : :
41 : : namespace MBMesquite
42 : : {
43 : :
44 : 4360452 : bool PMeanPMetric::average( PatchData& pd, QualityMetric* metric, const std::vector< size_t >& qm_handles,
45 : : double& value, MsqError& err )
46 : : {
47 : 4360452 : bool rval = true;
48 : 4360452 : value = 0;
49 [ + - ][ + - ]: 13168572 : for( std::vector< size_t >::const_iterator i = qm_handles.begin(); i != qm_handles.end(); ++i )
[ + + ]
50 : : {
51 : : double mval;
52 [ + - ][ + - ]: 8808152 : if( !metric->evaluate( pd, *i, mval, err ) ) rval = false;
[ + + ]
53 [ + - ][ + + ]: 8808152 : MSQ_ERRZERO( err );
[ + - ][ + - ]
[ + + ]
54 [ + - ]: 8808120 : value += P.raise( mval );
55 : : }
56 : 4360420 : value /= qm_handles.size();
57 : 4360452 : return rval;
58 : : }
59 : :
60 : 0 : bool PMeanPMetric::average_with_indices( PatchData& pd, QualityMetric* metric, const std::vector< size_t >& qm_handles,
61 : : double& value, std::vector< size_t >& indices, MsqError& err )
62 : : {
63 : 0 : indices.clear();
64 : :
65 : 0 : bool rval = true;
66 : 0 : value = 0;
67 [ # # ][ # # ]: 0 : for( std::vector< size_t >::const_iterator i = qm_handles.begin(); i != qm_handles.end(); ++i )
[ # # ]
68 : : {
69 : : double mval;
70 : 0 : mIndices.clear();
71 [ # # ][ # # ]: 0 : if( !metric->evaluate_with_indices( pd, *i, mval, mIndices, err ) ) rval = false;
[ # # ]
72 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
73 [ # # ]: 0 : value += P.raise( mval );
74 : :
75 [ # # ][ # # ]: 0 : std::copy( mIndices.begin(), mIndices.end(), std::back_inserter( indices ) );
76 : : }
77 : 0 : std::sort( indices.begin(), indices.end() );
78 : 0 : indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
79 : :
80 : 0 : value /= qm_handles.size();
81 : 0 : return rval;
82 : : }
83 : :
84 : 741906 : bool PMeanPMetric::average_with_gradient( PatchData& pd, QualityMetric* metric, const std::vector< size_t >& qm_handles,
85 : : double& value, std::vector< size_t >& indices,
86 : : std::vector< Vector3D >& gradient, MsqError& err )
87 : : {
88 : 741906 : indices.clear();
89 : 741906 : gradient.clear();
90 : :
91 : 741906 : std::vector< Vector3D >::iterator g;
92 : 741906 : std::vector< size_t >::iterator j, k;
93 : 741906 : std::vector< size_t >::const_iterator i;
94 : :
95 : 741906 : bool rval = true;
96 : 741906 : value = 0;
97 [ + - ][ + - ]: 2806311 : for( i = qm_handles.begin(); i != qm_handles.end(); ++i )
[ + + ]
98 : : {
99 : :
100 : : double mval;
101 : 2064405 : mIndices.clear();
102 : 2064405 : mGrad.clear();
103 [ + - ][ + - ]: 2064405 : if( !metric->evaluate_with_gradient( pd, *i, mval, mIndices, mGrad, err ) ) rval = false;
[ - + ]
104 [ + - ][ - + ]: 2064405 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
105 [ + - ]: 2064405 : value += P.raise( mval );
106 : :
107 [ + - ][ + - ]: 2064405 : double p1val = P.value() * P1.raise( mval );
108 [ + - ][ + - ]: 4684911 : for( j = mIndices.begin(), g = mGrad.begin(); j != mIndices.end(); ++j, ++g )
[ + - ][ + + ]
109 : : {
110 : :
111 [ + - ][ + - ]: 2620506 : *g *= p1val;
112 [ + - ][ + - ]: 2620506 : k = std::lower_bound( indices.begin(), indices.end(), *j );
113 [ + - ][ + + ]: 2620506 : if( k == indices.end() || *k != *j )
[ + - ][ + - ]
[ + + ][ + - ]
[ + + # # ]
114 : : {
115 [ + - ][ + - ]: 1427016 : k = indices.insert( k, *j );
116 [ + - ]: 1427016 : size_t idx = k - indices.begin();
117 [ + - ][ + - ]: 1427016 : gradient.insert( gradient.begin() + idx, *g );
[ + - ]
118 : : }
119 : : else
120 : : {
121 [ + - ]: 1193490 : size_t idx = k - indices.begin();
122 [ + - ][ + - ]: 1193490 : gradient[idx] += *g;
[ + - ]
123 : : }
124 : : }
125 : : }
126 : :
127 : 741906 : double inv_n = 1.0 / qm_handles.size();
128 : 741906 : value *= inv_n;
129 [ + - ][ + - ]: 2168922 : for( g = gradient.begin(); g != gradient.end(); ++g )
[ + + ]
130 [ + - ][ + - ]: 1427016 : *g *= inv_n;
131 : :
132 : 741906 : return rval;
133 : : }
134 : :
135 : 0 : bool PMeanPMetric::average_with_Hessian_diagonal( PatchData& pd, QualityMetric* metric,
136 : : const std::vector< size_t >& qm_handles, double& value,
137 : : std::vector< size_t >& indices, std::vector< Vector3D >& gradient,
138 : : std::vector< SymMatrix3D >& diagonal, MsqError& err )
139 : : {
140 : : // clear temporary storage
141 : 0 : mIndices.clear();
142 : 0 : mGrad.clear();
143 : 0 : mDiag.clear();
144 : 0 : mOffsets.clear();
145 : 0 : mValues.clear();
146 : :
147 : : // Evaluate metric for all sample points,
148 : : // accumulating indices, gradients, and Hessians
149 : 0 : bool rval = true;
150 : 0 : std::vector< size_t >::const_iterator q;
151 [ # # ][ # # ]: 0 : for( q = qm_handles.begin(); q != qm_handles.end(); ++q )
[ # # ]
152 : : {
153 : : double mval;
154 : 0 : indices.clear();
155 : 0 : gradient.clear();
156 : 0 : diagonal.clear();
157 [ # # ][ # # ]: 0 : if( !metric->evaluate_with_Hessian_diagonal( pd, *q, mval, indices, gradient, diagonal, err ) ) rval = false;
[ # # ]
158 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
159 : :
160 [ # # ]: 0 : mValues.push_back( mval );
161 [ # # ]: 0 : mOffsets.push_back( mIndices.size() );
162 [ # # ][ # # ]: 0 : std::copy( indices.begin(), indices.end(), std::back_inserter( mIndices ) );
163 [ # # ][ # # ]: 0 : std::copy( gradient.begin(), gradient.end(), std::back_inserter( mGrad ) );
164 [ # # ][ # # ]: 0 : std::copy( diagonal.begin(), diagonal.end(), std::back_inserter( mDiag ) );
165 : : }
166 [ # # ]: 0 : mOffsets.push_back( mIndices.size() );
167 : :
168 : : // Combine lists of free vertex indices, and update indices
169 : : // in per-evaluation lists to point into the combined gradient
170 : : // and Hessian lists.
171 [ # # ]: 0 : indices = mIndices;
172 [ # # ]: 0 : std::sort( indices.begin(), indices.end() );
173 [ # # ][ # # ]: 0 : indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
174 : 0 : std::vector< size_t >::iterator i, j;
175 [ # # ][ # # ]: 0 : for( i = mIndices.begin(); i != mIndices.end(); ++i )
[ # # ]
176 : : {
177 [ # # ][ # # ]: 0 : j = std::lower_bound( indices.begin(), indices.end(), *i );
178 [ # # ][ # # ]: 0 : assert( *j == *i );
[ # # ]
179 [ # # ][ # # ]: 0 : *i = j - indices.begin();
180 : : }
181 : :
182 : : // Allocate space and zero output gradient and Hessian lists
183 : 0 : const size_t n = indices.size();
184 : 0 : const size_t m = mValues.size();
185 : 0 : const double inv_n = 1.0 / m;
186 : : double v;
187 : 0 : gradient.clear();
188 [ # # ][ # # ]: 0 : gradient.resize( n, Vector3D( 0, 0, 0 ) );
189 : 0 : diagonal.clear();
190 [ # # ][ # # ]: 0 : diagonal.resize( n, SymMatrix3D( 0.0 ) );
191 : :
192 : : // Average values, gradients, Hessians
193 : 0 : value = 0.0;
194 [ # # ]: 0 : for( size_t k = 0; k < m; ++k )
195 : : { // for each metric evaluate
196 [ # # ][ # # ]: 0 : if( P.value() == 1.0 )
197 : : {
198 [ # # ][ # # ]: 0 : v = P.raise( mValues[k] );
199 : : // for each gradient (or Hessian) for the local metric evaluation
200 [ # # ][ # # ]: 0 : for( size_t r = mOffsets[k]; r < mOffsets[k + 1]; ++r )
[ # # ]
201 : : {
202 [ # # ]: 0 : const size_t nr = mIndices[r];
203 : :
204 [ # # ][ # # ]: 0 : mDiag[r] *= inv_n;
205 [ # # ][ # # ]: 0 : diagonal[nr] += mDiag[r];
[ # # ]
206 [ # # ][ # # ]: 0 : mGrad[r] *= inv_n;
207 [ # # ][ # # ]: 0 : gradient[nr] += mGrad[r];
[ # # ]
208 : : }
209 : : }
210 : : else
211 : : {
212 [ # # ][ # # ]: 0 : const double r2 = P2.raise( mValues[k] );
213 [ # # ]: 0 : const double r1 = r2 * mValues[k];
214 [ # # ]: 0 : v = r1 * mValues[k];
215 [ # # ][ # # ]: 0 : const double h = r2 * P.value() * ( P.value() - 1 ) * inv_n;
216 [ # # ]: 0 : const double g = r1 * P.value() * inv_n;
217 : : // for each gradient (or Hessian) for the local metric evaluation
218 [ # # ][ # # ]: 0 : for( size_t r = mOffsets[k]; r < mOffsets[k + 1]; ++r )
[ # # ]
219 : : {
220 [ # # ]: 0 : const size_t nr = mIndices[r];
221 : :
222 [ # # ][ # # ]: 0 : mDiag[r] *= g;
223 [ # # ][ # # ]: 0 : diagonal[nr] += mDiag[r];
[ # # ]
224 [ # # ][ # # ]: 0 : diagonal[nr] += h * outer( mGrad[r] );
[ # # ][ # # ]
[ # # ]
225 [ # # ][ # # ]: 0 : mGrad[r] *= g;
226 [ # # ][ # # ]: 0 : gradient[nr] += mGrad[r];
[ # # ]
227 : : }
228 : : }
229 : 0 : value += v;
230 : : }
231 : :
232 : 0 : value *= inv_n;
233 : :
234 : 0 : return rval;
235 : : }
236 : :
237 : 124856 : bool PMeanPMetric::average_with_Hessian( PatchData& pd, QualityMetric* metric, const std::vector< size_t >& qm_handles,
238 : : double& value, std::vector< size_t >& indices,
239 : : std::vector< Vector3D >& gradient, std::vector< Matrix3D >& Hessian,
240 : : MsqError& err )
241 : : {
242 : : // clear temporary storage
243 : 124856 : mIndices.clear();
244 : 124856 : mGrad.clear();
245 : 124856 : mHess.clear();
246 : 124856 : mOffsets.clear();
247 : 124856 : mValues.clear();
248 : :
249 : : // Evaluate metric for all sample points,
250 : : // accumulating indices, gradients, and Hessians
251 : 124856 : bool rval = true;
252 : 124856 : std::vector< size_t >::const_iterator q;
253 [ + - ][ + - ]: 624280 : for( q = qm_handles.begin(); q != qm_handles.end(); ++q )
[ + + ]
254 : : {
255 : : double mval;
256 : 499424 : indices.clear();
257 : 499424 : gradient.clear();
258 : 499424 : Hessian.clear();
259 [ + - ][ + - ]: 499424 : if( !metric->evaluate_with_Hessian( pd, *q, mval, indices, gradient, Hessian, err ) ) rval = false;
[ - + ]
260 [ + - ][ - + ]: 499424 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
261 : :
262 [ + - ]: 499424 : mValues.push_back( mval );
263 [ + - ]: 499424 : mOffsets.push_back( mIndices.size() );
264 [ + - ][ + - ]: 499424 : std::copy( indices.begin(), indices.end(), std::back_inserter( mIndices ) );
265 [ + - ][ + - ]: 499424 : std::copy( gradient.begin(), gradient.end(), std::back_inserter( mGrad ) );
266 [ + - ][ + - ]: 499424 : std::copy( Hessian.begin(), Hessian.end(), std::back_inserter( mHess ) );
267 : : }
268 [ + - ]: 124856 : mOffsets.push_back( mIndices.size() );
269 : :
270 : : // Combine lists of free vertex indices, and update indices
271 : : // in per-evaluation lists to point into the combined gradient
272 : : // and Hessian lists.
273 [ + - ]: 124856 : indices = mIndices;
274 [ + - ]: 124856 : std::sort( indices.begin(), indices.end() );
275 [ + - ][ + - ]: 124856 : indices.erase( std::unique( indices.begin(), indices.end() ), indices.end() );
276 : 124856 : std::vector< size_t >::iterator i, j;
277 [ + - ][ + - ]: 1613564 : for( i = mIndices.begin(); i != mIndices.end(); ++i )
[ + + ]
278 : : {
279 [ + - ][ + - ]: 1488708 : j = std::lower_bound( indices.begin(), indices.end(), *i );
280 [ + - ][ + - ]: 1488708 : assert( *j == *i );
[ - + ]
281 [ + - ][ + - ]: 1488708 : *i = j - indices.begin();
282 : : }
283 : :
284 : : // Allocate space and zero output gradient and Hessian lists
285 : 124856 : const size_t N = indices.size();
286 : 124856 : const size_t m = mValues.size();
287 : 124856 : const double inv_n = 1.0 / m;
288 : : double v;
289 : 124856 : gradient.clear();
290 [ + - ][ + - ]: 124856 : gradient.resize( N, Vector3D( 0, 0, 0 ) );
291 : 124856 : Hessian.clear();
292 [ + - ][ + - ]: 124856 : Hessian.resize( N * ( N + 1 ) / 2, Matrix3D( 0.0 ) );
293 : :
294 : : // Average values, gradients, Hessians
295 [ + - ]: 124856 : Matrix3D outer;
296 : 124856 : value = 0.0;
297 : 124856 : std::vector< Matrix3D >::iterator met_hess_iter = mHess.begin();
298 [ + + ]: 624280 : for( size_t k = 0; k < m; ++k )
299 : : { // for each metric evaluate
300 [ + - ][ + - ]: 499424 : if( P.value() == 1.0 )
301 : : {
302 [ + - ][ + - ]: 499424 : v = P.raise( mValues[k] );
303 : : // for each gradient (or Hessian row) for the local metric evaluation
304 [ + - ][ + - ]: 1988132 : for( size_t r = mOffsets[k]; r < mOffsets[k + 1]; ++r )
[ + + ]
305 : : {
306 [ + - ]: 1488708 : const size_t nr = mIndices[r];
307 : : // for each column of the local metric Hessian
308 [ + - ][ + + ]: 4459904 : for( size_t c = r; c < mOffsets[k + 1]; ++c )
309 : : {
310 [ + - ]: 2971196 : const size_t nc = mIndices[c];
311 [ + - ][ + - ]: 2971196 : *met_hess_iter *= inv_n;
312 [ + + ]: 2971196 : if( nr <= nc )
313 [ + - ][ + - ]: 2229952 : Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += *met_hess_iter;
[ + - ]
314 : : else
315 [ + - ][ + - ]: 741244 : Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( *met_hess_iter );
[ + - ]
316 [ + - ]: 2971196 : ++met_hess_iter;
317 : : }
318 [ + - ][ + - ]: 1488708 : mGrad[r] *= inv_n;
319 [ + - ][ + - ]: 1488708 : gradient[nr] += mGrad[r];
[ + - ]
320 : : }
321 : : }
322 : : else
323 : : {
324 [ # # ][ # # ]: 0 : const double r2 = P2.raise( mValues[k] );
325 [ # # ]: 0 : const double r1 = r2 * mValues[k];
326 [ # # ]: 0 : v = r1 * mValues[k];
327 [ # # ][ # # ]: 0 : const double h = r2 * P.value() * ( P.value() - 1 ) * inv_n;
328 [ # # ]: 0 : const double g = r1 * P.value() * inv_n;
329 : : // for each gradient (or Hessian row) for the local metric evaluation
330 [ # # ][ # # ]: 0 : for( size_t r = mOffsets[k]; r < mOffsets[k + 1]; ++r )
[ # # ]
331 : : {
332 [ # # ]: 0 : const size_t nr = mIndices[r];
333 : : // for each column of the local metric Hessian
334 [ # # ][ # # ]: 0 : for( size_t c = r; c < mOffsets[k + 1]; ++c )
335 : : {
336 [ # # ]: 0 : const size_t nc = mIndices[c];
337 [ # # ][ # # ]: 0 : *met_hess_iter *= g;
338 [ # # ][ # # ]: 0 : outer.outer_product( mGrad[r], mGrad[c] );
[ # # ]
339 [ # # ]: 0 : outer *= h;
340 [ # # ][ # # ]: 0 : outer += *met_hess_iter;
341 [ # # ]: 0 : if( nr <= nc )
342 [ # # ][ # # ]: 0 : Hessian[N * nr - nr * ( nr + 1 ) / 2 + nc] += outer;
343 : : else
344 [ # # ][ # # ]: 0 : Hessian[N * nc - nc * ( nc + 1 ) / 2 + nr].plus_transpose_equal( outer );
345 [ # # ]: 0 : ++met_hess_iter;
346 : : }
347 [ # # ][ # # ]: 0 : mGrad[r] *= g;
348 [ # # ][ # # ]: 0 : gradient[nr] += mGrad[r];
[ # # ]
349 : : }
350 : : }
351 : 499424 : value += v;
352 : : }
353 : :
354 : 124856 : value *= inv_n;
355 : :
356 : 124856 : return rval;
357 : : }
358 : :
359 [ + - ][ + - ]: 120 : } // namespace MBMesquite
|
