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