MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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