C++ QualityMetric::evaluate_with_Hessian_diagonal方法代码示例

bool LPtoPTemplate::evaluate_with_Hessian_diagonal( EvalType type, 
                                        PatchData& pd,
                                        double& OF_val,
                                        std::vector<Vector3D>& grad,
                                        std::vector<SymMatrix3D>& hess_diag,
                                        MsqError& err )
{
  QualityMetric* qm = get_quality_metric();
  qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );  MSQ_ERRFALSE(err);
  
    // zero gradient and hessian
  grad.clear();
  grad.resize( pd.num_free_vertices(), 0.0 );
  hess_diag.clear();
  hess_diag.resize( pd.num_free_vertices(), 0.0 );
  
  double QM_val, QM_pow = 1.0;
  double fac1, fac2;
  const double negate_flag = qm->get_negate_flag();
  bool qm_bool;
  size_t i;
  short p;
   
  // Loops over all elements in the patch.
  OF_val = 0.0;
  std::vector<size_t>::const_iterator k;
  for (k = qmHandles.begin(); k != qmHandles.end(); ++k)
  {
    // Computes \nabla^2 Q(e). Only the free vertices will have non-zero entries. 
    qm_bool = qm->evaluate_with_Hessian_diagonal( pd, *k, QM_val, mIndices, mGradient, mDiag, err );
    if (MSQ_CHKERR(err) || !qm_bool) return false;
    QM_val = fabs(QM_val);

    // **** Computes Hessian ****
    const size_t nve = mIndices.size();
    if (pVal == 1) {
      QM_pow = 1.0;
      for (i=0; i<nve; ++i) {
        mDiag[i] *= negate_flag;
        hess_diag[mIndices[i]] += mDiag[i];
      }
      fac1 = 1;
    }
    else if (pVal >= 2) {
       // Computes the coefficients:
      QM_pow = 1.0;
      for (p=0; p<pVal-2; ++p)
        QM_pow *= QM_val;
      // 1 - computes p(p-1)Q(e)^{p-2}
      fac2 = pVal* (pVal-1) * QM_pow;
      // 2 - computes  pQ(e)^{p-1}
      QM_pow *= QM_val;
      fac1 = pVal * QM_pow;

        //fac1 *= qm->get_negate_flag();
        //fac2 *= qm->get_negate_flag();

      for (i=0; i<nve; ++i) {
        SymMatrix3D op(mGradient[i]);
        op *= fac2;
        mDiag[i] *= fac1;
        op += mDiag[i];
        op *= negate_flag;
        hess_diag[mIndices[i]] += op;
      }
    } else {
      MSQ_SETERR(err)(" invalid P value.", MsqError::INVALID_STATE);
      return false;
    }


    // **** Computes Gradient ****

    // For each vertex in the element ... 
    for (i=0; i<nve; ++i) {
      // ... computes p*q^{p-1}*grad(q) ...
      mGradient[i] *= fac1*negate_flag;
      // ... and accumulates it in the objective function gradient.
        //also scale the gradient by the scaling factor
      assert (mIndices[i] < pd.num_free_vertices());
      grad[mIndices[i]] += mGradient[i];
    }

    // **** computes Objective Function value \sum_{i=1}^{N_e} |q_i|^P ****
    OF_val += QM_pow * QM_val;
  }

  size_t global_count;
  OF_val = negate_flag
         * get_value( OF_val, qmHandles.size(), type, global_count, err );
//  if (!global_count)
//    return false;  // invalid mesh

  if (dividingByN && global_count) {
    const double inv_n = 1.0 / global_count;
    for (i = 0; i < pd.num_free_vertices(); ++i) {
      grad[i] *= inv_n;
      hess_diag[i] *= inv_n;
    }
  }
//.........这里部分代码省略.........

bool PMeanPTemplate::evaluate_with_Hessian_diagonal( EvalType type, 
                                        PatchData& pd,
                                        double& value_out,
                                        std::vector<Vector3D>& grad_out,
                                        std::vector<SymMatrix3D>& hess_diag_out,
                                        MsqError& err )
{
  QualityMetric* qm = get_quality_metric();
  qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );  MSQ_ERRFALSE(err);
  
    // zero gradient and hessian
  const size_t s = pd.num_free_vertices();
  grad_out.clear();
  grad_out.resize( s, 0.0 );
  hess_diag_out.clear();
  hess_diag_out.resize( s, 0.0 );
  
    // calculate OF value and gradient for just the patch
  std::vector<size_t>::const_iterator i;
  size_t j;
  double value, working_sum = 0.0;
  const double f1 = qm->get_negate_flag() * mPower.value();
  const double f2 = f1 * (mPower.value() - 1);
  for (i = qmHandles.begin(); i != qmHandles.end(); ++i)
  {
    bool result = qm->evaluate_with_Hessian_diagonal( pd, *i, value, mIndices, mGradient, mDiag, err );
    if (MSQ_CHKERR(err) || !result)
      return false;
    if (fabs(value) < DBL_EPSILON)
      continue;
    
    const size_t nfree = mIndices.size();
    if (mPower.value() == 1.0) {
      working_sum += mPower.raise( value );
      for (j = 0; j < nfree; ++j) {
        const size_t idx = mIndices[j];
        hess_diag_out[idx] += f1 * mDiag[j];
        mGradient[j] *= f1;
        grad_out[idx] += mGradient[j];
      }
    }
    else {
      const double r2 = mPowerMinus2.raise( value );
      const double r1 = r2 * value;
      working_sum += r1 * value;
      const double hf = f2 * r2;
      const double gf = f1 * r1;
      for (j = 0; j < nfree; ++j) {
        const size_t idx = mIndices[j];

        hess_diag_out[idx] += hf * outer( mGradient[j] );
        hess_diag_out[idx] += gf * mDiag[j];

        mGradient[j] *= gf;
        grad_out[idx] += mGradient[j];
      }
    }
  }
  
    // get overall OF value, update member data, etc.
  size_t global_count;
  value_out = qm->get_negate_flag() 
            * get_value( working_sum, qmHandles.size(), type, global_count );
  const double inv_n = 1.0 / global_count;
  for (j = 0; j < s; ++j) {
    grad_out[j] *= inv_n;
    hess_diag_out[j] *= inv_n;
  }
  return true;
}

bool CompareMetric::evaluate_with_Hessian_diagonal( PatchData& pd,
                 size_t handle,
                 double& value,
                 std::vector<size_t>& indices,
                 std::vector<Vector3D>& gradient,
                 std::vector<SymMatrix3D>& diagonal,
                 MsqError& err )
{
  double m2val;
  bool r1, r2;
  m2Handles.clear();
  m2Grad.clear();
  m2Diag.clear();
  r1 = metric1->evaluate_with_Hessian_diagonal( pd, handle, value, indices, gradient, diagonal, err ); MSQ_ERRZERO(err);
  r2 = metric2->evaluate_with_Hessian_diagonal( pd, handle, m2val, m2Handles, m2Grad, m2Diag, err ); MSQ_ERRZERO(err);
  if (r1 != r2 || (r1 && fabs(value - m2val) > epsilon)) {
    MSQ_SETERR(err)(MsqError::INVALID_STATE,
                    "Metrics returned different values for handle %lu in "
                    "evaluate_with_Hessian_diagonal:\n"
                    "\t%s %f vs. %s %f\n", (unsigned long)handle,
                    r1?"true":"false",value,r2?"true":"false",m2val);
  }
  else {
    std::vector<size_t>::const_iterator i, j;
    std::vector<Vector3D>::const_iterator r, s;
    std::vector<SymMatrix3D>::const_iterator u, v;
    int grad_diff = 0, hess_diff = 0;
    bool same = (indices.size() == m2Handles.size());
    std::sort( m2Handles.begin(), m2Handles.end() );
    for (i = indices.begin(); i != indices.end(); ++i) {
      j = std::lower_bound( m2Handles.begin(), m2Handles.end(), *i );
      if (j == m2Handles.end() || *j != *i) {
        same = false;
        continue;
      }
      
      r = gradient.begin() + (i - indices.begin());
      s = m2Grad.begin() + (j - m2Handles.begin());
      if (!equal(*r,*s))
        ++grad_diff;
      
      u = diagonal.begin() + (i - indices.begin());
      v = m2Diag.begin() + (j - m2Handles.begin());
      if (!equal(*u,*v))
        ++hess_diff;
    }
      
    if (!same) {
      MSQ_SETERR(err)(MsqError::INVALID_STATE,
                      "Metrics returned incompatible lists of vertex indices"
                      " for handle %lu in evaluate_with_Hessian_diagonal\n.", 
                      (unsigned long)handle );
    }
    else if (grad_diff) {
      MSQ_SETERR(err)(MsqError::INVALID_STATE,
                      "Metrics returned different gradient vectors for "
                      " %d of %u vertices for handle %lu in "
                      "evaluate_with_Hessian_diagonal\n.", 
                      grad_diff, (unsigned)gradient.size(), 
                      (unsigned long)handle );
    }
    else if (hess_diff) {
      MSQ_SETERR(err)(MsqError::INVALID_STATE,
                      "Metrics returned different Hessian blocks for "
                      " %d of %u vertices for handle %lu in "
                      "evaluate_with_Hessian_diagonal\n.", 
                      hess_diff, (unsigned)diagonal.size(), 
                      (unsigned long)handle );
    }
  }
  
  return r1 && !err;
}