MOAB: Mesh Oriented datABase  (version 5.4.0)
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,
00056                                    double sqr_sum,
00057                                    size_t count,
00058                                    EvalType type,
00059                                    double& result_sum,
00060                                    double& result_sqr,
00061                                    size_t& global_count )
00062 {
00063     switch( type )
00064     {
00065         case CALCULATE:
00066             result_sum   = sum;
00067             result_sqr   = sqr_sum;
00068             global_count = count;
00069             break;
00070 
00071         case ACCUMULATE:
00072             result_sum   = mSum += sum;
00073             result_sqr   = mSqrSum += sqr_sum;
00074             global_count = mCount += count;
00075             break;
00076 
00077         case SAVE:
00078             saveSum      = sum;
00079             saveSqrSum   = sqr_sum;
00080             saveCount    = count;
00081             result_sum   = mSum;
00082             result_sqr   = mSqrSum;
00083             global_count = mCount;
00084             break;
00085 
00086         case UPDATE:
00087             mSum -= saveSum;
00088             mSqrSum -= saveSqrSum;
00089             mCount -= saveCount;
00090             result_sum = mSum += saveSum = sum;
00091             result_sqr = mSqrSum += saveSqrSum = sqr_sum;
00092             global_count = mCount += saveCount = count;
00093             break;
00094 
00095         case TEMPORARY:
00096             result_sum   = mSum - saveSum + sum;
00097             result_sqr   = mSqrSum - saveSqrSum + sqr_sum;
00098             global_count = mCount + count - saveCount;
00099             break;
00100     }
00101 }
00102 
00103 bool VarianceTemplate::evaluate( EvalType type, PatchData& pd, double& value_out, bool free, MsqError& err )
00104 {
00105     QualityMetric* qm = get_quality_metric();
00106     if( type == ObjectiveFunction::ACCUMULATE )
00107         qm->get_single_pass( pd, qmHandles, free, err );
00108     else
00109         qm->get_evaluations( pd, qmHandles, free, err );
00110     MSQ_ERRFALSE( err );
00111 
00112     // calculate OF value for just the patch
00113     std::vector< size_t >::const_iterator i;
00114     double value, sum = 0.0, sqr = 0.0;
00115     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00116     {
00117         bool result = qm->evaluate( pd, *i, value, err );
00118         if( MSQ_CHKERR( err ) || !result ) return false;
00119 
00120         sum += value;
00121         sqr += value * value;
00122     }
00123 
00124     // get overall OF value, update member data, etc.
00125     size_t n;
00126     accumulate( sum, sqr, qmHandles.size(), type, sum, sqr, n );
00127 
00128     // don't divide by zero
00129     if( n < 1 )
00130     {
00131         value_out = 0.0;
00132         return true;
00133     }
00134 
00135     const double rms = sqr / n;
00136     const double avg = sum / n;
00137     value_out        = qm->get_negate_flag() * ( rms - avg * avg );
00138     return true;
00139 }
00140 
00141 bool VarianceTemplate::evaluate_with_gradient( EvalType type,
00142                                                PatchData& pd,
00143                                                double& value_out,
00144                                                std::vector< Vector3D >& grad_out,
00145                                                MsqError& err )
00146 {
00147     QualityMetric* qm = get_quality_metric();
00148     qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );
00149     MSQ_ERRFALSE( err );
00150 
00151     // zero gradient data
00152     grad_out.clear();  // store sum of metric * gradient of metric, and later OF gradient
00153     grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00154     gradSum.clear();  // store sum of gradients of metrics
00155     gradSum.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00156 
00157     // calculate OF value and gradient for just the patch
00158     Matrix3D op;
00159     std::vector< size_t >::const_iterator i;
00160     double value, sum = 0.0, sqr = 0.0;
00161     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00162     {
00163         bool result = qm->evaluate_with_gradient( pd, *i, value, mIndices, mGradient, err );
00164         if( MSQ_CHKERR( err ) || !result ) return false;
00165         if( fabs( value ) < DBL_EPSILON ) continue;
00166 
00167         sum += value;
00168         sqr += value * value;
00169 
00170         for( size_t j = 0; j < mIndices.size(); ++j )
00171         {
00172             const size_t r = mIndices[j];
00173             gradSum[r] += mGradient[j];
00174             mGradient[j] *= value;
00175             grad_out[r] += mGradient[j];
00176         }
00177     }
00178 
00179     // update accumulated values (if doing BCD)
00180     size_t n;
00181     accumulate( sum, sqr, qmHandles.size(), type, sum, sqr, n );
00182 
00183     // don't divide by zero
00184     if( n < 1 )
00185     {
00186         value_out = 0.0;
00187         grad_out.clear();
00188         grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00189         return true;
00190     }
00191 
00192     // calculate OF value
00193     const double rms = sqr / n;
00194     const double avg = sum / n;
00195     value_out        = qm->get_negate_flag() * ( rms - avg * avg );
00196 
00197     // Finish calculation of gradient and Hessian
00198     const double f = qm->get_negate_flag() * 2.0 / n;
00199     for( size_t k = 0; k < pd.num_free_vertices(); ++k )
00200     {
00201         gradSum[k] *= avg;
00202         grad_out[k] -= gradSum[k];
00203         grad_out[k] *= f;
00204     }
00205 
00206     return true;
00207 }
00208 
00209 bool VarianceTemplate::evaluate_with_Hessian_diagonal( EvalType type,
00210                                                        PatchData& pd,
00211                                                        double& value_out,
00212                                                        std::vector< Vector3D >& grad_out,
00213                                                        std::vector< SymMatrix3D >& hess_diag_out,
00214                                                        MsqError& err )
00215 {
00216     QualityMetric* qm = get_quality_metric();
00217     qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );
00218     MSQ_ERRFALSE( err );
00219 
00220     // zero gradient and Hessian data
00221     grad_out.clear();  // store sum of metric * gradient of metric, and later OF gradient
00222     grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00223     gradSum.clear();  // store sum of gradients of metrics
00224     gradSum.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00225     hessSum.clear();  // store sum of Hessians of metrics
00226     hessSum.resize( pd.num_free_vertices(), SymMatrix3D( 0.0 ) );
00227     hess_diag_out.clear();  // store sum of metric * outer_product(metric gradient), and later OF Hessian
00228     hess_diag_out.resize( pd.num_free_vertices(), SymMatrix3D( 0.0 ) );
00229 
00230     // calculate OF value and gradient for just the patch
00231     Matrix3D op;
00232     std::vector< size_t >::const_iterator i;
00233     double value, sum = 0.0, sqr = 0.0;
00234     for( i = qmHandles.begin(); i != qmHandles.end(); ++i )
00235     {
00236         bool result = qm->evaluate_with_Hessian_diagonal( pd, *i, value, mIndices, mGradient, mHessDiag, err );
00237         if( MSQ_CHKERR( err ) || !result ) return false;
00238         if( fabs( value ) < DBL_EPSILON ) continue;
00239 
00240         sum += value;
00241         sqr += value * value;
00242 
00243         for( size_t j = 0; j < mIndices.size(); ++j )
00244         {
00245             const size_t r = mIndices[j];
00246 
00247             hessSum[r] += mHessDiag[j];
00248             hess_diag_out[r] += outer( mGradient[j] );
00249             mHessDiag[j] *= value;
00250             hess_diag_out[r] += mHessDiag[j];
00251 
00252             gradSum[r] += mGradient[j];
00253             mGradient[j] *= value;
00254             grad_out[r] += mGradient[j];
00255         }
00256     }
00257 
00258     // update accumulated values (if doing BCD)
00259     size_t n;
00260     accumulate( sum, sqr, qmHandles.size(), type, sum, sqr, n );
00261 
00262     // don't divide by zero
00263     if( n < 1 )
00264     {
00265         value_out = 0.0;
00266         grad_out.clear();
00267         grad_out.resize( pd.num_free_vertices(), Vector3D( 0.0, 0.0, 0.0 ) );
00268         hess_diag_out.clear();
00269         hess_diag_out.resize( pd.num_free_vertices(), SymMatrix3D( 0.0 ) );
00270         return true;
00271     }
00272 
00273     // calculate OF value
00274     const double rms = sqr / n;
00275     const double avg = sum / n;
00276     value_out        = qm->get_negate_flag() * ( rms - avg * avg );
00277 
00278     // Finish calculation of gradient and Hessian
00279     const double inv_n = 1.0 / n;
00280     const double f     = qm->get_negate_flag() * 2.0 / n;
00281     for( size_t k = 0; k < pd.num_free_vertices(); ++k )
00282     {
00283         hessSum[k] *= avg;
00284         hess_diag_out[k] -= hessSum[k];
00285         hess_diag_out[k] -= inv_n * outer( gradSum[k] );
00286         hess_diag_out[k] *= f;
00287 
00288         gradSum[k] *= avg;
00289         grad_out[k] -= gradSum[k];
00290         grad_out[k] *= f;
00291     }
00292 
00293     return true;
00294 }
00295 
00296 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines