MOAB: Mesh Oriented datABase  (version 5.2.1)
VarianceTemplate.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retain certain rights to 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     (2006) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file VarianceTemplate.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "VarianceTemplate.hpp"
00034 #include "QualityMetric.hpp"
00035 #include "MsqError.hpp"
00036 #include "MsqHessian.hpp"
00037 #include "PatchData.hpp"
00038 
00039 namespace MBMesquite
00040 {
00041 
00042 ObjectiveFunction* VarianceTemplate::clone() const
00043 {
00044     return new VarianceTemplate( *this );
00045 }
00046 
00047 void VarianceTemplate::clear()
00048 {
00049     mCount = 0;
00050     mSum = mSqrSum = 0;
00051     saveCount      = 0;
00052     saveSum = saveSqrSum = 0;
00053 }
00054 
00055 void VarianceTemplate::accumulate( double sum, double sqr_sum, size_t count, EvalType type, double& result_sum,
00056                                    double& result_sqr, size_t& global_count )
00057 {
00058     switch( type )
00059     {
00060         case CALCULATE:
00061             result_sum   = sum;
00062             result_sqr   = sqr_sum;
00063             global_count = count;
00064             break;
00065 
00066         case ACCUMULATE:
00067             result_sum   = mSum += sum;
00068             result_sqr   = mSqrSum += sqr_sum;
00069             global_count = mCount += count;
00070             break;
00071 
00072         case SAVE:
00073             saveSum      = sum;
00074             saveSqrSum   = sqr_sum;
00075             saveCount    = count;
00076             result_sum   = mSum;
00077             result_sqr   = mSqrSum;
00078             global_count = mCount;
00079             break;
00080 
00081         case UPDATE:
00082             mSum -= saveSum;
00083             mSqrSum -= saveSqrSum;
00084             mCount -= saveCount;
00085             result_sum = mSum += saveSum = sum;
00086             result_sqr = mSqrSum += saveSqrSum = sqr_sum;
00087             global_count = mCount += saveCount = count;
00088             break;
00089 
00090         case TEMPORARY:
00091             result_sum   = mSum - saveSum + sum;
00092             result_sqr   = mSqrSum - saveSqrSum + sqr_sum;
00093             global_count = mCount + count - saveCount;
00094             break;
00095     }
00096 }
00097 
00098 bool VarianceTemplate::evaluate( EvalType type, PatchData& pd, double& value_out, bool free, MsqError& err )
00099 {
00100     QualityMetric* qm = get_quality_metric();
00101     if( type == ObjectiveFunction::ACCUMULATE )
00102         qm->get_single_pass( pd, qmHandles, free, err );
00103     else
00104         qm->get_evaluations( pd, qmHandles, free, err );
00105     MSQ_ERRFALSE( err );
00106 
00107     // calculate OF value for just the patch
00108     std::vector< size_t >::const_iterator i;
00109     double value, sum = 0.0, sqr = 0.0;
00110     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00111     {
00112         bool result = qm->evaluate( pd, *i, value, err );
00113         if( MSQ_CHKERR( err ) || !result ) return false;
00114 
00115         sum += value;
00116         sqr += value * value;
00117     }
00118 
00119     // get overall OF value, update member data, etc.
00120     size_t n;
00121     accumulate( sum, sqr, qmHandles.size(), type, sum, sqr, n );
00122 
00123     // don't divide by zero
00124     if( n < 1 )
00125     {
00126         value_out = 0.0;
00127         return true;
00128     }
00129 
00130     const double rms = sqr / n;
00131     const double avg = sum / n;
00132     value_out        = qm->get_negate_flag() * ( rms - avg * avg );
00133     return true;
00134 }
00135 
00136 bool VarianceTemplate::evaluate_with_gradient( EvalType type, PatchData& pd, double& value_out,
00137                                                std::vector< Vector3D >& grad_out, MsqError& err )
00138 {
00139     QualityMetric* qm = get_quality_metric();
00140     qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );
00141     MSQ_ERRFALSE( err );
00142 
00143     // zero gradient data
00144     grad_out.clear();  // store sum of metric * gradient of metric, and later OF gradient
00145     grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00146     gradSum.clear();  // store sum of gradients of metrics
00147     gradSum.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00148 
00149     // calculate OF value and gradient for just the patch
00150     Matrix3D op;
00151     std::vector< size_t >::const_iterator i;
00152     double value, sum = 0.0, sqr = 0.0;
00153     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00154     {
00155         bool result = qm->evaluate_with_gradient( pd, *i, value, mIndices, mGradient, err );
00156         if( MSQ_CHKERR( err ) || !result ) return false;
00157         if( fabs( value ) < DBL_EPSILON ) continue;
00158 
00159         sum += value;
00160         sqr += value * value;
00161 
00162         for( size_t j = 0; j < mIndices.size(); ++j )
00163         {
00164             const size_t r = mIndices[j];
00165             gradSum[r] += mGradient[j];
00166             mGradient[j] *= value;
00167             grad_out[r] += mGradient[j];
00168         }
00169     }
00170 
00171     // update accumulated values (if doing BCD)
00172     size_t n;
00173     accumulate( sum, sqr, qmHandles.size(), type, sum, sqr, n );
00174 
00175     // don't divide by zero
00176     if( n < 1 )
00177     {
00178         value_out = 0.0;
00179         grad_out.clear();
00180         grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00181         return true;
00182     }
00183 
00184     // calculate OF value
00185     const double rms = sqr / n;
00186     const double avg = sum / n;
00187     value_out        = qm->get_negate_flag() * ( rms - avg * avg );
00188 
00189     // Finish calculation of gradient and Hessian
00190     const double f = qm->get_negate_flag() * 2.0 / n;
00191     for( size_t k = 0; k < pd.num_free_vertices(); ++k )
00192     {
00193         gradSum[k] *= avg;
00194         grad_out[k] -= gradSum[k];
00195         grad_out[k] *= f;
00196     }
00197 
00198     return true;
00199 }
00200 
00201 bool VarianceTemplate::evaluate_with_Hessian_diagonal( EvalType type, PatchData& pd, double& value_out,
00202                                                        std::vector< Vector3D >& grad_out,
00203                                                        std::vector< SymMatrix3D >& hess_diag_out, MsqError& err )
00204 {
00205     QualityMetric* qm = get_quality_metric();
00206     qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );
00207     MSQ_ERRFALSE( err );
00208 
00209     // zero gradient and Hessian data
00210     grad_out.clear();  // store sum of metric * gradient of metric, and later OF gradient
00211     grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00212     gradSum.clear();  // store sum of gradients of metrics
00213     gradSum.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00214     hessSum.clear();  // store sum of Hessians of metrics
00215     hessSum.resize( pd.num_free_vertices(), SymMatrix3D( 0.0 ) );
00216     hess_diag_out.clear();  // store sum of metric * outer_product(metric gradient), and later OF Hessian
00217     hess_diag_out.resize( pd.num_free_vertices(), SymMatrix3D( 0.0 ) );
00218 
00219     // calculate OF value and gradient for just the patch
00220     Matrix3D op;
00221     std::vector< size_t >::const_iterator i;
00222     double value, sum = 0.0, sqr = 0.0;
00223     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00224     {
00225         bool result = qm->evaluate_with_Hessian_diagonal( pd, *i, value, mIndices, mGradient, mHessDiag, err );
00226         if( MSQ_CHKERR( err ) || !result ) return false;
00227         if( fabs( value ) < DBL_EPSILON ) continue;
00228 
00229         sum += value;
00230         sqr += value * value;
00231 
00232         for( size_t j = 0; j < mIndices.size(); ++j )
00233         {
00234             const size_t r = mIndices[j];
00235 
00236             hessSum[r] += mHessDiag[j];
00237             hess_diag_out[r] += outer( mGradient[j] );
00238             mHessDiag[j] *= value;
00239             hess_diag_out[r] += mHessDiag[j];
00240 
00241             gradSum[r] += mGradient[j];
00242             mGradient[j] *= value;
00243             grad_out[r] += mGradient[j];
00244         }
00245     }
00246 
00247     // update accumulated values (if doing BCD)
00248     size_t n;
00249     accumulate( sum, sqr, qmHandles.size(), type, sum, sqr, n );
00250 
00251     // don't divide by zero
00252     if( n < 1 )
00253     {
00254         value_out = 0.0;
00255         grad_out.clear();
00256         grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00257         hess_diag_out.clear();
00258         hess_diag_out.resize( pd.num_free_vertices(), SymMatrix3D( 0.0 ) );
00259         return true;
00260     }
00261 
00262     // calculate OF value
00263     const double rms = sqr / n;
00264     const double avg = sum / n;
00265     value_out        = qm->get_negate_flag() * ( rms - avg * avg );
00266 
00267     // Finish calculation of gradient and Hessian
00268     const double inv_n = 1.0 / n;
00269     const double f     = qm->get_negate_flag() * 2.0 / n;
00270     for( size_t k = 0; k < pd.num_free_vertices(); ++k )
00271     {
00272         hessSum[k] *= avg;
00273         hess_diag_out[k] -= hessSum[k];
00274         hess_diag_out[k] -= inv_n * outer( gradSum[k] );
00275         hess_diag_out[k] *= f;
00276 
00277         gradSum[k] *= avg;
00278         grad_out[k] -= gradSum[k];
00279         grad_out[k] *= f;
00280     }
00281 
00282     return true;
00283 }
00284 
00285 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines