MOAB: Mesh Oriented datABase  (version 5.2.1)
AveragingQM.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
00025 
00026     (2006) kraftche@cae.wisc.edu
00027 
00028   ***************************************************************** */
00029 /*!
00030   \file   AveragingQM.cpp
00031   \brief
00032 
00033   \author Michael Brewer
00034   \author Thomas Leurent
00035   \author Jason Kraftcheck
00036   \date   2002-05-14
00037 */
00038 
00039 #include "AveragingQM.hpp"
00040 #include "MsqVertex.hpp"
00041 #include "MsqMeshEntity.hpp"
00042 #include "MsqDebug.hpp"
00043 #include "MsqTimer.hpp"
00044 #include "PatchData.hpp"
00045 
00046 namespace MBMesquite
00047 {
00048 
00049 double AveragingQM::average_corner_gradients( EntityTopology type, uint32_t fixed_vertices, unsigned num_corner,
00050                                               double corner_values[], const Vector3D corner_grads[],
00051                                               Vector3D vertex_grads[], MsqError& err )
00052 {
00053     const unsigned num_vertex = TopologyInfo::corners( type );
00054     const unsigned dim        = TopologyInfo::dimension( type );
00055     const unsigned per_vertex = dim + 1;
00056 
00057     unsigned i, j, num_adj;
00058     const unsigned *adj_idx, *rev_idx;
00059 
00060     // NOTE: This function changes the corner_values array such that
00061     //       it contains the gradient coefficients.
00062     double avg = average_metric_and_weights( corner_values, num_corner, err );
00063     MSQ_ERRZERO( err );
00064 
00065     for( i = 0; i < num_vertex; ++i )
00066     {
00067         if( fixed_vertices & ( 1 << i ) )  // skip fixed vertices
00068             continue;
00069 
00070         adj_idx = TopologyInfo::adjacent_vertices( type, i, num_adj );
00071         rev_idx = TopologyInfo::reverse_vertex_adjacency_offsets( type, i, num_adj );
00072         if( i < num_corner )  // not all vertices are corners (e.g. pyramid)
00073             vertex_grads[i] = corner_values[i] * corner_grads[per_vertex * i];
00074         else
00075             vertex_grads[i] = 0;
00076         for( j = 0; j < num_adj; ++j )
00077         {
00078             const unsigned v = adj_idx[j], c = rev_idx[j] + 1;
00079             if( v >= num_corner )  // if less corners than vertices (e.g. pyramid apex)
00080                 continue;
00081             vertex_grads[i] += corner_values[v] * corner_grads[per_vertex * v + c];
00082         }
00083     }
00084 
00085     return avg;
00086 }
00087 
00088 /**\brief Iterate over only diagonal blocks of element corner Hessian data
00089  *
00090  * Given concatenation of corner Hessian data for an element, iterate
00091  * over only the diagonal terms for each corner.  This class allows
00092  * common code to be used to generate Hessian diagonal blocks from either
00093  * the diagonal blocks for each corner or the full Hessian data for each
00094  * corner, where this class is used for the latter.
00095  */
00096 class CornerHessDiagIterator
00097 {
00098   private:
00099     const Matrix3D* cornerHess;     //!< Current location in concatenated Hessian data.
00100     const EntityTopology elemType;  //!< Element topology for Hessian data
00101     unsigned mCorner;               //!< The element corner for which cornerHess
00102                                     //!< is pointing into the corresponding Hessian data.
00103     unsigned mStep;                 //!< Amount to step to reach next diagonal block.
00104   public:
00105     CornerHessDiagIterator( const Matrix3D* corner_hessians, EntityTopology elem_type )
00106         : cornerHess( corner_hessians ), elemType( elem_type ), mCorner( 0 )
00107     {
00108         TopologyInfo::adjacent_vertices( elemType, mCorner, mStep );
00109         ++mStep;
00110     }
00111 
00112     SymMatrix3D operator*() const
00113     {
00114         return cornerHess->upper();
00115     }
00116 
00117     CornerHessDiagIterator& operator++()
00118     {
00119         cornerHess += mStep;
00120         if( !--mStep )
00121         {
00122             TopologyInfo::adjacent_vertices( elemType, ++mCorner, mStep );
00123             ++mStep;
00124         }
00125         return *this;
00126     }
00127 
00128     CornerHessDiagIterator operator++( int )
00129     {
00130         CornerHessDiagIterator copy( *this );
00131         operator++();
00132         return copy;
00133     }
00134 };
00135 
00136 template < typename HessIter >
00137 static inline double sum_corner_diagonals( EntityTopology type, unsigned num_corner, const double corner_values[],
00138                                            const Vector3D corner_grads[], HessIter corner_diag_blocks,
00139                                            Vector3D vertex_grads[], SymMatrix3D vertex_hessians[] )
00140 {
00141     unsigned i, n, r, R, idx[4];
00142     const unsigned* adj_list;
00143     double avg = 0.0;
00144 
00145     // calculate mean
00146     for( i = 0; i < num_corner; ++i )
00147         avg += corner_values[i];
00148 
00149     const Vector3D* grad = corner_grads;
00150     HessIter hess        = corner_diag_blocks;
00151     for( i = 0; i < num_corner; ++i )
00152     {
00153         adj_list = TopologyInfo::adjacent_vertices( type, i, n );
00154         idx[0]   = i;
00155         idx[1]   = adj_list[0];
00156         idx[2]   = adj_list[1];
00157         idx[3]   = adj_list[2 % n];  // %n so don't read off end if 2D
00158 
00159         for( r = 0; r <= n; ++r )
00160         {
00161             R = idx[r];
00162             vertex_grads[R] += *grad;
00163             vertex_hessians[R] += *hess;
00164             ++grad;
00165             ++hess;
00166         }
00167     }
00168     return avg;
00169 }
00170 
00171 template < typename HessIter >
00172 static inline double sum_sqr_corner_diagonals( EntityTopology type, unsigned num_corner, const double corner_values[],
00173                                                const Vector3D corner_grads[], HessIter corner_diag_blocks,
00174                                                Vector3D vertex_grads[], SymMatrix3D vertex_hessians[] )
00175 {
00176     unsigned i, n, r, R, idx[4];
00177     const unsigned* adj_list;
00178     double v, avg = 0.0;
00179 
00180     // calculate mean
00181     for( i = 0; i < num_corner; ++i )
00182         avg += corner_values[i] * corner_values[i];
00183 
00184     const Vector3D* grad = corner_grads;
00185     HessIter hess        = corner_diag_blocks;
00186     for( i = 0; i < num_corner; ++i )
00187     {
00188         adj_list = TopologyInfo::adjacent_vertices( type, i, n );
00189         idx[0]   = i;
00190         idx[1]   = adj_list[0];
00191         idx[2]   = adj_list[1];
00192         idx[3]   = adj_list[2 % n];  // %n so don't read off end if 2D
00193         ++n;
00194 
00195         v = 2.0 * corner_values[i];
00196         for( r = 0; r < n; ++r )
00197         {
00198             R = idx[r];
00199             vertex_grads[R] += v * *grad;
00200             vertex_hessians[R] += 2.0 * outer( *grad );
00201             vertex_hessians[R] += v * *hess;
00202             ++grad;
00203             ++hess;
00204         }
00205     }
00206     return avg;
00207 }
00208 
00209 template < typename HessIter >
00210 static inline double pmean_corner_diagonals( EntityTopology type, unsigned num_corner, const double corner_values[],
00211                                              const Vector3D corner_grads[], HessIter corner_diag_blocks,
00212                                              Vector3D vertex_grads[], SymMatrix3D vertex_hessians[], double p )
00213 {
00214     const unsigned N = TopologyInfo::corners( type );
00215     unsigned i, n, r, R, idx[4];
00216     const unsigned* adj_list;
00217     double m = 0.0, nm;
00218     double gf[8], hf[8];
00219     double inv = 1.0 / num_corner;
00220     assert( num_corner <= 8 );
00221 
00222     // calculate mean
00223     for( i = 0; i < num_corner; ++i )
00224     {
00225         nm = pow( corner_values[i], p );
00226         m += nm;
00227 
00228         gf[i] = inv * p * nm / corner_values[i];
00229         hf[i] = ( p - 1 ) * gf[i] / corner_values[i];
00230     }
00231     nm = inv * m;
00232 
00233     const Vector3D* grad = corner_grads;
00234     HessIter hess        = corner_diag_blocks;
00235     for( i = 0; i < num_corner; ++i )
00236     {
00237         adj_list = TopologyInfo::adjacent_vertices( type, i, n );
00238         idx[0]   = i;
00239         idx[1]   = adj_list[0];
00240         idx[2]   = adj_list[1];
00241         idx[3]   = adj_list[2 % n];  // %n so don't read off end if 2D
00242         ++n;
00243 
00244         for( r = 0; r < n; ++r )
00245         {
00246             R = idx[r];
00247             vertex_grads[R] += gf[i] * *grad;
00248             vertex_hessians[R] += hf[i] * outer( *grad );
00249             vertex_hessians[R] += gf[i] * *hess;
00250             ++grad;
00251             ++hess;
00252         }
00253     }
00254 
00255     m     = pow( nm, 1.0 / p );
00256     gf[0] = m / ( p * nm );
00257     hf[0] = ( 1.0 / p - 1 ) * gf[0] / nm;
00258     for( r = 0; r < N; ++r )
00259     {
00260         vertex_hessians[r] *= gf[0];
00261         vertex_hessians[r] += hf[0] * outer( vertex_grads[r] );
00262         vertex_grads[r] *= gf[0];
00263     }
00264 
00265     return m;
00266 }
00267 
00268 template < typename HessIter >
00269 static inline double average_corner_diagonals( EntityTopology type, QualityMetric::AveragingMethod method,
00270                                                unsigned num_corner, const double corner_values[],
00271                                                const Vector3D corner_grads[], HessIter corner_diag_blocks,
00272                                                Vector3D vertex_grads[], SymMatrix3D vertex_hessians[], MsqError& err )
00273 {
00274     unsigned i;
00275     double avg, inv;
00276 
00277     // Zero gradients and Hessians
00278     const unsigned num_vertex = TopologyInfo::corners( type );
00279     for( i = 0; i < num_vertex; ++i )
00280     {
00281         vertex_grads[i].set( 0.0 );
00282         vertex_hessians[i] = SymMatrix3D( 0.0 );
00283     }
00284 
00285     switch( method )
00286     {
00287         case QualityMetric::SUM:
00288             avg = sum_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks, vertex_grads,
00289                                         vertex_hessians );
00290             break;
00291 
00292         case QualityMetric::LINEAR:
00293             avg = sum_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks, vertex_grads,
00294                                         vertex_hessians );
00295             inv = 1.0 / num_corner;
00296             avg *= inv;
00297             for( i = 0; i < num_vertex; ++i )
00298             {
00299                 vertex_grads[i] *= inv;
00300                 vertex_hessians[i] *= inv;
00301             }
00302             break;
00303 
00304         case QualityMetric::SUM_SQUARED:
00305             avg = sum_sqr_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks,
00306                                             vertex_grads, vertex_hessians );
00307             break;
00308 
00309         case QualityMetric::RMS:
00310             avg = pmean_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks,
00311                                           vertex_grads, vertex_hessians, 2.0 );
00312             break;
00313 
00314         case QualityMetric::HARMONIC:
00315             avg = pmean_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks,
00316                                           vertex_grads, vertex_hessians, -1.0 );
00317             break;
00318 
00319         case QualityMetric::HMS:
00320             avg = pmean_corner_diagonals( type, num_corner, corner_values, corner_grads, corner_diag_blocks,
00321                                           vertex_grads, vertex_hessians, -2.0 );
00322             break;
00323 
00324         default:
00325             MSQ_SETERR( err )( "averaging method not available.", MsqError::INVALID_STATE );
00326             return 0.0;
00327     }
00328 
00329     return avg;
00330 }
00331 
00332 double AveragingQM::average_corner_hessian_diagonals( EntityTopology element_type, uint32_t, unsigned num_corners,
00333                                                       const double corner_values[], const Vector3D corner_grads[],
00334                                                       const Matrix3D corner_hessians[], Vector3D vertex_grads[],
00335                                                       SymMatrix3D vertex_hessians[], MsqError& err )
00336 {
00337     return average_corner_diagonals( element_type, avgMethod, num_corners, corner_values, corner_grads,
00338                                      CornerHessDiagIterator( corner_hessians, element_type ), vertex_grads,
00339                                      vertex_hessians, err );
00340 }
00341 
00342 double AveragingQM::average_corner_hessian_diagonals( EntityTopology element_type, uint32_t, unsigned num_corners,
00343                                                       const double corner_values[], const Vector3D corner_grads[],
00344                                                       const SymMatrix3D corner_hess_diag[], Vector3D vertex_grads[],
00345                                                       SymMatrix3D vertex_hessians[], MsqError& err )
00346 {
00347     return average_corner_diagonals( element_type, avgMethod, num_corners, corner_values, corner_grads,
00348                                      corner_hess_diag, vertex_grads, vertex_hessians, err );
00349 }
00350 
00351 static inline double sum_corner_hessians( EntityTopology type, unsigned num_corner, const double corner_values[],
00352                                           const Vector3D corner_grads[], const Matrix3D corner_hessians[],
00353                                           Vector3D vertex_grads[], Matrix3D vertex_hessians[] )
00354 {
00355     const unsigned N = TopologyInfo::corners( type );
00356     unsigned i, n, r, c, R, C, idx[4];
00357     const unsigned* adj_list;
00358     double avg = 0.0;
00359 
00360     // calculate mean
00361     for( i = 0; i < num_corner; ++i )
00362         avg += corner_values[i];
00363 
00364     const Vector3D* grad = corner_grads;
00365     const Matrix3D* hess = corner_hessians;
00366     for( i = 0; i < num_corner; ++i )
00367     {
00368         adj_list = TopologyInfo::adjacent_vertices( type, i, n );
00369         idx[0]   = i;
00370         idx[1]   = adj_list[0];
00371         idx[2]   = adj_list[1];
00372         idx[3]   = adj_list[2 % n];  // %n so don't read off end if 2D
00373 
00374         for( r = 0; r <= n; ++r )
00375         {
00376             R = idx[r];
00377             vertex_grads[R] += *grad;
00378             ++grad;
00379             for( c = r; c <= n; ++c )
00380             {
00381                 C = idx[c];
00382                 if( R <= C )
00383                     vertex_hessians[N * R - R * ( R + 1 ) / 2 + C] += *hess;
00384                 else
00385                     vertex_hessians[N * C - C * ( C + 1 ) / 2 + R].plus_transpose_equal( *hess );
00386                 ++hess;
00387             }
00388         }
00389     }
00390     return avg;
00391 }
00392 
00393 static inline double sum_sqr_corner_hessians( EntityTopology type, unsigned num_corner, const double corner_values[],
00394                                               const Vector3D corner_grads[], const Matrix3D corner_hessians[],
00395                                               Vector3D vertex_grads[], Matrix3D vertex_hessians[] )
00396 {
00397     const unsigned N = TopologyInfo::corners( type );
00398     unsigned i, n, r, c, R, C, idx[4];
00399     const unsigned* adj_list;
00400     double v, avg = 0.0;
00401     Matrix3D op;
00402 
00403     // calculate mean
00404     for( i = 0; i < num_corner; ++i )
00405         avg += corner_values[i] * corner_values[i];
00406 
00407     const Vector3D* grad = corner_grads;
00408     const Matrix3D* hess = corner_hessians;
00409     for( i = 0; i < num_corner; ++i )
00410     {
00411         adj_list = TopologyInfo::adjacent_vertices( type, i, n );
00412         idx[0]   = i;
00413         idx[1]   = adj_list[0];
00414         idx[2]   = adj_list[1];
00415         idx[3]   = adj_list[2 % n];  // %n so don't read off end if 2D
00416         ++n;
00417 
00418         v = 2.0 * corner_values[i];
00419         for( r = 0; r < n; ++r )
00420         {
00421             R = idx[r];
00422             vertex_grads[R] += v * grad[r];
00423             for( c = r; c < n; ++c )
00424             {
00425                 C = idx[c];
00426                 op.outer_product( 2.0 * grad[r], grad[c] );
00427                 op += v * *hess;
00428                 if( R <= C )
00429                     vertex_hessians[N * R - R * ( R + 1 ) / 2 + C] += op;
00430                 else
00431                     vertex_hessians[N * C - C * ( C + 1 ) / 2 + R].plus_transpose_equal( op );
00432                 ++hess;
00433             }
00434         }
00435         grad += n;
00436     }
00437     return avg;
00438 }
00439 
00440 static inline double pmean_corner_hessians( EntityTopology type, unsigned num_corner, const double corner_values[],
00441                                             const Vector3D corner_grads[], const Matrix3D corner_hessians[],
00442                                             Vector3D vertex_grads[], Matrix3D vertex_hessians[], double p )
00443 {
00444     const unsigned N = TopologyInfo::corners( type );
00445     unsigned i, n, r, c, R, C, idx[4];
00446     const unsigned* adj_list;
00447     double m = 0.0, nm;
00448     Matrix3D op;
00449     double gf[8], hf[8];
00450     double inv = 1.0 / num_corner;
00451     assert( num_corner <= 8 );
00452 
00453     // calculate mean
00454     for( i = 0; i < num_corner; ++i )
00455     {
00456         nm = pow( corner_values[i], p );
00457         m += nm;
00458 
00459         gf[i] = inv * p * nm / corner_values[i];
00460         hf[i] = ( p - 1 ) * gf[i] / corner_values[i];
00461     }
00462     nm = inv * m;
00463 
00464     const Vector3D* grad = corner_grads;
00465     const Matrix3D* hess = corner_hessians;
00466     for( i = 0; i < num_corner; ++i )
00467     {
00468         adj_list = TopologyInfo::adjacent_vertices( type, i, n );
00469         idx[0]   = i;
00470         idx[1]   = adj_list[0];
00471         idx[2]   = adj_list[1];
00472         idx[3]   = adj_list[2 % n];  // %n so don't read off end if 2D
00473         ++n;
00474 
00475         for( r = 0; r < n; ++r )
00476         {
00477             R = idx[r];
00478             vertex_grads[R] += gf[i] * grad[r];
00479             for( c = r; c < n; ++c )
00480             {
00481                 C = idx[c];
00482                 op.outer_product( grad[r], grad[c] );
00483                 op *= hf[i];
00484                 op += gf[i] * *hess;
00485                 if( R <= C )
00486                     vertex_hessians[N * R - R * ( R + 1 ) / 2 + C] += op;
00487                 else
00488                     vertex_hessians[N * C - C * ( C + 1 ) / 2 + R].plus_transpose_equal( op );
00489                 ++hess;
00490             }
00491         }
00492         grad += n;
00493     }
00494 
00495     m     = pow( nm, 1.0 / p );
00496     gf[0] = m / ( p * nm );
00497     hf[0] = ( 1.0 / p - 1 ) * gf[0] / nm;
00498     for( r = 0; r < N; ++r )
00499     {
00500         for( c = r; c < N; ++c )
00501         {
00502             op.outer_product( vertex_grads[r], vertex_grads[c] );
00503             op *= hf[0];
00504             *vertex_hessians *= gf[0];
00505             *vertex_hessians += op;
00506             ++vertex_hessians;
00507         }
00508         vertex_grads[r] *= gf[0];
00509     }
00510 
00511     return m;
00512 }
00513 
00514 double AveragingQM::average_corner_hessians( EntityTopology type, uint32_t, unsigned num_corner,
00515                                              const double corner_values[], const Vector3D corner_grads[],
00516                                              const Matrix3D corner_hessians[], Vector3D vertex_grads[],
00517                                              Matrix3D vertex_hessians[], MsqError& err )
00518 {
00519     unsigned i;
00520     double avg, inv;
00521 
00522     // Zero gradients and Hessians
00523     const unsigned num_vertex = TopologyInfo::corners( type );
00524     for( i = 0; i < num_vertex; ++i )
00525         vertex_grads[i].set( 0.0 );
00526     const unsigned num_hess = num_vertex * ( num_vertex + 1 ) / 2;
00527     for( i = 0; i < num_hess; ++i )
00528         vertex_hessians[i].zero();
00529 
00530     switch( avgMethod )
00531     {
00532         case QualityMetric::SUM:
00533             avg = sum_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads,
00534                                        vertex_hessians );
00535             break;
00536 
00537         case QualityMetric::LINEAR:
00538             avg = sum_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads,
00539                                        vertex_hessians );
00540             inv = 1.0 / num_corner;
00541             avg *= inv;
00542             for( i = 0; i < num_vertex; ++i )
00543                 vertex_grads[i] *= inv;
00544             for( i = 0; i < num_hess; ++i )
00545                 vertex_hessians[i] *= inv;
00546             break;
00547 
00548         case QualityMetric::SUM_SQUARED:
00549             avg = sum_sqr_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads,
00550                                            vertex_hessians );
00551             break;
00552 
00553         case QualityMetric::RMS:
00554             avg = pmean_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads,
00555                                          vertex_hessians, 2.0 );
00556             break;
00557 
00558         case QualityMetric::HARMONIC:
00559             avg = pmean_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads,
00560                                          vertex_hessians, -1.0 );
00561             break;
00562 
00563         case QualityMetric::HMS:
00564             avg = pmean_corner_hessians( type, num_corner, corner_values, corner_grads, corner_hessians, vertex_grads,
00565                                          vertex_hessians, -2.0 );
00566             break;
00567 
00568         default:
00569             MSQ_SETERR( err )( "averaging method not available.", MsqError::INVALID_STATE );
00570             return 0.0;
00571     }
00572 
00573     return avg;
00574 }
00575 
00576 double AveragingQM::average_metric_and_weights( double metrics[], int count, MsqError& err )
00577 {
00578     static bool min_max_warning = false;
00579     double avg                  = 0.0;
00580     int i, tmp_count;
00581     double f;
00582 
00583     switch( avgMethod )
00584     {
00585 
00586         case QualityMetric::MINIMUM:
00587             if( !min_max_warning )
00588             {
00589                 MSQ_DBGOUT( 1 ) << "Minimum and maximum not continuously differentiable.\n"
00590                                    "Element of subdifferential will be returned.\n";
00591                 min_max_warning = true;
00592             }
00593 
00594             avg = metrics[0];
00595             for( i = 1; i < count; ++i )
00596                 if( metrics[i] < avg ) avg = metrics[i];
00597 
00598             tmp_count = 0;
00599             for( i = 0; i < count; ++i )
00600             {
00601                 if( metrics[i] - avg <= MSQ_MIN )
00602                 {
00603                     metrics[i] = 1.0;
00604                     ++tmp_count;
00605                 }
00606                 else
00607                 {
00608                     metrics[i] = 0.0;
00609                 }
00610             }
00611 
00612             f = 1.0 / tmp_count;
00613             for( i = 0; i < count; ++i )
00614                 metrics[i] *= f;
00615 
00616             break;
00617 
00618         case QualityMetric::MAXIMUM:
00619             if( !min_max_warning )
00620             {
00621                 MSQ_DBGOUT( 1 ) << "Minimum and maximum not continuously differentiable.\n"
00622                                    "Element of subdifferential will be returned.\n";
00623                 min_max_warning = true;
00624             }
00625 
00626             avg = metrics[0];
00627             for( i = 1; i < count; ++i )
00628                 if( metrics[i] > avg ) avg = metrics[i];
00629 
00630             tmp_count = 0;
00631             for( i = 0; i < count; ++i )
00632             {
00633                 if( avg - metrics[i] <= MSQ_MIN )
00634                 {
00635                     metrics[i] = 1.0;
00636                     ++tmp_count;
00637                 }
00638                 else
00639                 {
00640                     metrics[i] = 0.0;
00641                 }
00642             }
00643 
00644             f = 1.0 / tmp_count;
00645             for( i = 0; i < count; ++i )
00646                 metrics[i] *= f;
00647 
00648             break;
00649 
00650         case QualityMetric::SUM:
00651             for( i = 0; i < count; ++i )
00652             {
00653                 avg += metrics[i];
00654                 metrics[i] = 1.0;
00655             }
00656 
00657             break;
00658 
00659         case QualityMetric::SUM_SQUARED:
00660             for( i = 0; i < count; ++i )
00661             {
00662                 avg += ( metrics[i] * metrics[i] );
00663                 metrics[i] *= 2;
00664             }
00665 
00666             break;
00667 
00668         case QualityMetric::LINEAR:
00669             f = 1.0 / count;
00670             for( i = 0; i < count; ++i )
00671             {
00672                 avg += metrics[i];
00673                 metrics[i] = f;
00674             }
00675             avg *= f;
00676 
00677             break;
00678 
00679         case QualityMetric::GEOMETRIC:
00680             avg = 1.0;
00681             for( i = 0; i < count; ++i )
00682                 avg *= metrics[i];
00683             avg = pow( avg, 1.0 / count );
00684 
00685             f = avg / count;
00686             for( i = 0; i < count; ++i )
00687                 metrics[i] = f / metrics[i];
00688 
00689             break;
00690 
00691         case QualityMetric::RMS:
00692             for( i = 0; i < count; ++i )
00693                 avg += metrics[i] * metrics[i];
00694             avg = sqrt( avg / count );
00695 
00696             f = 1. / ( avg * count );
00697             for( i = 0; i < count; ++i )
00698                 metrics[i] *= f;
00699 
00700             break;
00701 
00702         case QualityMetric::HARMONIC:
00703             for( i = 0; i < count; ++i )
00704                 avg += 1.0 / metrics[i];
00705             avg = count / avg;
00706 
00707             for( i = 0; i < count; ++i )
00708                 metrics[i] = ( avg * avg ) / ( count * metrics[i] * metrics[i] );
00709 
00710             break;
00711 
00712         case QualityMetric::HMS:
00713             for( i = 0; i < count; ++i )
00714                 avg += 1. / ( metrics[i] * metrics[i] );
00715             avg = sqrt( count / avg );
00716 
00717             f = avg * avg * avg / count;
00718             for( i = 0; i < count; ++i )
00719                 metrics[i] = f / ( metrics[i] * metrics[i] * metrics[i] );
00720 
00721             break;
00722 
00723         default:
00724             MSQ_SETERR( err )( "averaging method not available.", MsqError::INVALID_STATE );
00725     }
00726 
00727     return avg;
00728 }
00729 
00730 /*!
00731   average_metrics takes an array of length num_value and averages the
00732   contents using averaging method 'method'.
00733 */
00734 double AveragingQM::average_metrics( const double metric_values[], int num_values, MsqError& err )
00735 {
00736     // MSQ_MAX needs to be made global?
00737     // double MSQ_MAX=1e10;
00738     double total_value = 0.0;
00739     double temp_value  = 0.0;
00740     int i              = 0;
00741     int j              = 0;
00742     // if no values, return zero
00743     if( num_values <= 0 ) { return 0.0; }
00744 
00745     switch( avgMethod )
00746     {
00747         case QualityMetric::GEOMETRIC:
00748             total_value = 1.0;
00749             for( i = 0; i < num_values; ++i )
00750             {
00751                 total_value *= ( metric_values[i] );
00752             }
00753             total_value = pow( total_value, 1.0 / num_values );
00754             break;
00755 
00756         case QualityMetric::HARMONIC:
00757             // ensure no divide by zero, return zero
00758             for( i = 0; i < num_values; ++i )
00759             {
00760                 if( metric_values[i] < MSQ_MIN )
00761                 {
00762                     if( metric_values[i] > MSQ_MIN ) { return 0.0; }
00763                 }
00764                 total_value += ( 1 / metric_values[i] );
00765             }
00766             // ensure no divide by zero, return MSQ_MAX_CAP
00767             if( total_value < MSQ_MIN )
00768             {
00769                 if( total_value > MSQ_MIN ) { return MSQ_MAX_CAP; }
00770             }
00771             total_value = num_values / total_value;
00772             break;
00773 
00774         case QualityMetric::LINEAR:
00775             for( i = 0; i < num_values; ++i )
00776             {
00777                 total_value += metric_values[i];
00778             }
00779             total_value /= (double)num_values;
00780             break;
00781 
00782         case QualityMetric::MAXIMUM:
00783             total_value = metric_values[0];
00784             for( i = 1; i < num_values; ++i )
00785             {
00786                 if( metric_values[i] > total_value ) { total_value = metric_values[i]; }
00787             }
00788             break;
00789 
00790         case QualityMetric::MINIMUM:
00791             total_value = metric_values[0];
00792             for( i = 1; i < num_values; ++i )
00793             {
00794                 if( metric_values[i] < total_value ) { total_value = metric_values[i]; }
00795             }
00796             break;
00797 
00798         case QualityMetric::RMS:
00799             for( i = 0; i < num_values; ++i )
00800             {
00801                 total_value += ( metric_values[i] * metric_values[i] );
00802             }
00803             total_value /= (double)num_values;
00804             total_value = sqrt( total_value );
00805             break;
00806 
00807         case QualityMetric::HMS:
00808             // ensure no divide by zero, return zero
00809             for( i = 0; i < num_values; ++i )
00810             {
00811                 if( metric_values[i] * metric_values[i] < MSQ_MIN ) { return 0.0; }
00812                 total_value += ( 1.0 / ( metric_values[i] * metric_values[i] ) );
00813             }
00814 
00815             // ensure no divide by zero, return MSQ_MAX_CAP
00816             if( total_value < MSQ_MIN ) { return MSQ_MAX_CAP; }
00817             total_value = sqrt( num_values / total_value );
00818             break;
00819 
00820         case QualityMetric::STANDARD_DEVIATION:
00821             total_value = 0;
00822             temp_value  = 0;
00823             for( i = 0; i < num_values; ++i )
00824             {
00825                 temp_value += metric_values[i];
00826                 total_value += ( metric_values[i] * metric_values[i] );
00827             }
00828             temp_value /= (double)num_values;
00829             temp_value *= temp_value;
00830             total_value /= (double)num_values;
00831             total_value = total_value - temp_value;
00832             break;
00833 
00834         case QualityMetric::SUM:
00835             for( i = 0; i < num_values; ++i )
00836             {
00837                 total_value += metric_values[i];
00838             }
00839             break;
00840 
00841         case QualityMetric::SUM_SQUARED:
00842             for( i = 0; i < num_values; ++i )
00843             {
00844                 total_value += ( metric_values[i] * metric_values[i] );
00845             }
00846             break;
00847 
00848         case QualityMetric::MAX_MINUS_MIN:
00849             // total_value used to store the maximum
00850             // temp_value used to store the minimum
00851             temp_value = MSQ_MAX_CAP;
00852             for( i = 0; i < num_values; ++i )
00853             {
00854                 if( metric_values[i] < temp_value ) { temp_value = metric_values[i]; }
00855                 if( metric_values[i] > total_value ) { total_value = metric_values[i]; }
00856             }
00857 
00858             total_value -= temp_value;
00859             break;
00860 
00861         case QualityMetric::MAX_OVER_MIN:
00862             // total_value used to store the maximum
00863             // temp_value used to store the minimum
00864             temp_value = MSQ_MAX_CAP;
00865             for( i = 0; i < num_values; ++i )
00866             {
00867                 if( metric_values[i] < temp_value ) { temp_value = metric_values[i]; }
00868                 if( metric_values[i] > total_value ) { total_value = metric_values[i]; }
00869             }
00870 
00871             // ensure no divide by zero, return MSQ_MAX_CAP
00872             if( fabs( temp_value ) < MSQ_MIN ) { return MSQ_MAX_CAP; }
00873             total_value /= temp_value;
00874             break;
00875 
00876         case QualityMetric::SUM_OF_RATIOS_SQUARED:
00877             for( j = 0; j < num_values; ++j )
00878             {
00879                 // ensure no divide by zero, return MSQ_MAX_CAP
00880                 if( fabs( metric_values[j] ) < MSQ_MIN ) { return MSQ_MAX_CAP; }
00881                 for( i = 0; i < num_values; ++i )
00882                 {
00883                     total_value +=
00884                         ( ( metric_values[i] / metric_values[j] ) * ( metric_values[i] / metric_values[j] ) );
00885                 }
00886             }
00887             total_value /= (double)( num_values * num_values );
00888             break;
00889 
00890         default:
00891             // Return error saying Averaging Method mode not implemented
00892             MSQ_SETERR( err )
00893             ( "Requested Averaging Method Not Implemented", MsqError::NOT_IMPLEMENTED );
00894             return 0;
00895     }
00896     return total_value;
00897 }
00898 
00899 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines