C++ QualityMetric类代码示例

bool PMeanPTemplate::evaluate( EvalType type, 
                           PatchData& pd,
                           double& value_out,
                           bool free,
                           MsqError& err )
{
  QualityMetric* qm = get_quality_metric();
  if (type == ObjectiveFunction::ACCUMULATE)
    qm->get_single_pass( pd, qmHandles, free, err );
  else
    qm->get_evaluations( pd, qmHandles, free, err );  
  MSQ_ERRFALSE(err);
  
    // calculate OF value for just the patch
  std::vector<size_t>::const_iterator i;
  double value, working_sum = 0.0;
  for (i = qmHandles.begin(); i != qmHandles.end(); ++i)
  {
    bool result = qm->evaluate( pd, *i, value, err );
    if (MSQ_CHKERR(err) || !result)
      return false;
    
    working_sum += mPower.raise( value );
  }
  
    // 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 );
  return true;
}

bool LInfTemplate::evaluate( EvalType type, 
                             PatchData& pd,
                             double& value_out,
                             bool free,
                             MsqError& err )
{
  if (type != ObjectiveFunction::CALCULATE) {
    MSQ_SETERR(err)(
      "LInfTemplate does not support block coodinate descent algoritms",
      MsqError::INVALID_STATE );
    return false;
  }

  QualityMetric* qm = get_quality_metric();
  qm->get_evaluations( pd, qmHandles, free, err );  MSQ_ERRFALSE(err);
  const double negate = qm->get_negate_flag();
  
    // calculate OF value for just the patch
  std::vector<size_t>::const_iterator i;
  double value;
  value_out = -HUGE_VAL;
  for (i = qmHandles.begin(); i != qmHandles.end(); ++i)
  {
    bool result = qm->evaluate( pd, *i, value, err );
    if (MSQ_CHKERR(err) || !result)
      return false;

    value = negate * fabs(value);
    if (value > value_out)
      value_out = value;
  }
  
  return true;
}

bool LPtoPTemplate::evaluate_with_gradient( EvalType type, 
                                            PatchData& pd,
                                            double& OF_val,
                                            std::vector<Vector3D>& grad_out,
                                            MsqError& err )
{
  QualityMetric* qm = get_quality_metric();
  qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );  MSQ_ERRFALSE(err);
  
    // zero gradient
  grad_out.clear();
  grad_out.resize( pd.num_free_vertices(), Vector3D(0.0,0.0,0.0) );
  bool qm_bool=true;
  double QM_val;
  OF_val = 0.;
  int p1;
  
    // calculate OF value and gradient for just the patch
  std::vector<size_t>::const_iterator i;
  for (i = qmHandles.begin(); i != qmHandles.end(); ++i)
  {
    qm_bool = qm->evaluate_with_gradient( pd, *i, QM_val, mIndices, mGradient, err );
    if (MSQ_CHKERR(err) || !qm_bool)
      return false;
    
    QM_val = fabs(QM_val);
    double QM_pow = 1.0;
    double factor = qm->get_negate_flag();
    if (pVal == 1) 
      QM_pow = 1.0;
    else {
      QM_pow = QM_val;
      for (p1 = 2; p1 < pVal; ++p1)
        QM_pow *= QM_val;
      factor *= QM_pow * pVal;
    }
    
    OF_val += QM_pow * QM_val;
    for (size_t j = 0; j < mIndices.size(); ++j) {
      mGradient[j] *= factor;
      grad_out[mIndices[j]] += mGradient[j];
    }
  }
  
    // get overall OF value, update member data, etc.
  size_t global_count;
  OF_val = qm->get_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;
    std::vector<Vector3D>::iterator g;
    for (g = grad_out.begin(); g != grad_out.end(); ++g)
      *g *= inv_n;
  }  
  
  return true;
}

bool LPtoPTemplate::concrete_evaluate(PatchData &pd, double &fval,
                                      MsqError &err){
  size_t index=0;
  MsqMeshEntity* elems=pd.get_element_array(err);
  bool obj_bool=true;
    //double check for pVal=0;
  if(pVal==0){
    MSQ_SETERR(err)("pVal equal zero not allowed.  L_0 is not a valid norm.",
                    MsqError::INVALID_STATE);
    return false;
  }
  
    //Michael:  this may not do what we want
    //Set currentQM to be the first quality metric* in the list
  QualityMetric* currentQM = get_quality_metric();
  if(currentQM==NULL)
    currentQM=get_quality_metric_list().front();
  if(currentQM==NULL) {
    MSQ_SETERR(err)("NULL QualityMetric pointer in LPtoPTemplate",
                    MsqError::INVALID_STATE);
    return false;
  }
  size_t num_elements=pd.num_elements();
  size_t num_vertices=pd.num_vertices();
  size_t total_num=0;
  if(currentQM->get_metric_type()==QualityMetric::ELEMENT_BASED)   
    total_num=num_elements;
  else if (currentQM->get_metric_type()==QualityMetric::VERTEX_BASED)
    total_num=num_vertices;
  else {
    MSQ_SETERR(err)("Make sure MetricType is initialised in concrete "
                    "QualityMetric constructor.", MsqError::INVALID_STATE);
    return false;
  }
  
  msq_std::vector<double> metric_values(total_num);
  if(currentQM->get_metric_type()==QualityMetric::ELEMENT_BASED)
  {
    for (index=0; index<num_elements;index++)
    {
        //if invalid return false after clean-up
      obj_bool=currentQM->evaluate_element(pd, (&elems[index]),
                                           metric_values[index], err);
      if(MSQ_CHKERR(err) || !obj_bool){
        fval=0.0;
        return false;
      }
      
      metric_values[index]=fabs(metric_values[index]);
      MSQ_DBGOUT(3) << "      o  Quality metric value for element "
                    << index << "\t: " << metric_values[index] << "\n";
    }
  }
  else if(currentQM->get_metric_type()==QualityMetric::VERTEX_BASED)
  {
    MsqVertex* vertices=pd.get_vertex_array(err); MSQ_ERRZERO(err);
    for (index=0; index<num_vertices;index++)
    {
        //evaluate metric for this vertex
      obj_bool=currentQM->evaluate_vertex(pd, (&vertices[index]),
                                          metric_values[index], err);
        //if invalid return false after clean-up
      if(MSQ_CHKERR(err) || !obj_bool){
        fval=0.0;
        return false;
      }
      
      metric_values[index]=fabs(metric_values[index]);
    }
  }
  fval=compute_function(&metric_values[0], total_num, err);
  return !MSQ_CHKERR(err);
}

/* virtual function reimplemented from QualityMetric. No doxygen doc needed. */
bool LPtoPTemplate::compute_analytical_gradient(PatchData &pd,
                                              Vector3D *const &grad,
                                              double &OF_val,
                                              MsqError &err, size_t array_size)
{
  MSQ_FUNCTION_TIMER( "LPtoPTemplate::compute_analytical_gradient" );
  
    //initialize the scaling value
  double scaling_value=1.0;
 
  size_t num_elements=pd.num_elements();
  size_t num_vertices=pd.num_vertices();
  if( num_vertices!=array_size && array_size>0)
  {
    MSQ_SETERR(err)("Incorrect array size.", MsqError::INVALID_ARG);
    return false;
  }
  
  MsqMeshEntity* elems=pd.get_element_array(err); MSQ_ERRZERO(err);
  MsqVertex* vertices=pd.get_vertex_array(err); MSQ_ERRZERO(err);
  bool qm_bool=true;
  double QM_val;
  OF_val = 0.;
  size_t i;
  int p1;
  
  //Set currentQM to be quality metric (possibly composite) associated with the objective function
  QualityMetric* currentQM = get_quality_metric();
  if(currentQM==NULL) {
    MSQ_SETERR(err)("LPtoPTemplate has NULL QualityMetric pointer.",MsqError::INVALID_STATE);
    return false;
  }
  enum QualityMetric::MetricType qm_type=currentQM->get_metric_type();
  
  if (qm_type!=QualityMetric::ELEMENT_BASED &&
      qm_type!=QualityMetric::VERTEX_BASED) {
    MSQ_SETERR(err)("Make sure MetricType is initialised"
                    "in concrete QualityMetric constructor.",
                    MsqError::INVALID_STATE);
    return false;
  }


  // zeros out objective function gradient
  for (i=0; i<num_vertices; ++i)
    grad[i] =0;
  
  // Computes objective function gradient for an element based metric
  if(qm_type==QualityMetric::ELEMENT_BASED){
      //if scaling, divid by num_elements
    if(dividingByN){
      if(num_elements<=0) {
        MSQ_SETERR(err)("The number of elements should not be zero.",MsqError::INVALID_MESH);
        return false;
      }
      
      scaling_value/=num_elements;
    }
    size_t e, ve;
    size_t nfve; // num free vtx in element
    size_t nve; // num vtx in element
    MsqVertex* ele_free_vtces[MSQ_MAX_NUM_VERT_PER_ENT];
    const size_t *ele_vtces_ind;

    // loops over all elements.
    for (e=0; e<num_elements; ++e) {
      // stores the pointers to the free vertices within the element
      // (using pointer arithmetic).
      nfve = 0;
      nve = elems[e].vertex_count();
      ele_vtces_ind = elems[e].get_vertex_index_array();
      for (ve=0; ve<nve; ++ve) {
        if (vertices[ele_vtces_ind[ve]].is_free_vertex()) {
          ele_free_vtces[nfve] = vertices + ele_vtces_ind[ve];
          ++nfve;
        }
      }

      // Computes q and grad(q)
      Vector3D grad_vec[MSQ_MAX_NUM_VERT_PER_ENT];
      qm_bool = currentQM->compute_element_gradient(
                                     pd, &elems[e],
                                     ele_free_vtces,
                                     grad_vec, nfve, QM_val, err);
      if(MSQ_CHKERR(err) || !qm_bool) return false;

      // computes p*|Q(e)|^{p-1}
      QM_val = fabs(QM_val);
      double QM_pow=1.0;
      double factor;
      if (pVal==1) factor=1;
      else {
        QM_pow=QM_val;
        for (p1=1; p1<pVal-1; ++p1)
          QM_pow*=QM_val;
        factor = QM_pow * pVal;
      }
        //this scales the gradient
      factor *= (scaling_value * get_negate_flag());
//.........这里部分代码省略.........

/*\ For each element, each entry to be accumulated in the Hessian for
    this objective function (\f$ \sum_{e \in E} Q(e)^p \f$ where \f$ E \f$
    is the set of all elements in the patch) has the form:
    \f$ pQ(e)^{p-1} \nabla^2 Q(e) + p(p-1)Q(e)^{p-2} \nabla Q(e) [\nabla Q(e)]^T \f$.

    For \f$ p=2 \f$, this simplifies to
    \f$ 2Q(e) \nabla^2 Q(e) + 2 \nabla Q(e) [\nabla Q(e)]^T \f$.

    For \f$ p=1 \f$, this simplifies to \f$ \nabla^2 Q(e) \f$.

    The \f$ p=1 \f$ simplified version is implemented directly
    to speed up computation. 

    This function does not support vertex-based metrics.
    
    \param pd The PatchData object for which the objective function
           hessian is computed.
    \param hessian this object must have been previously initialized.
*/
bool LPtoPTemplate::evaluate_with_Hessian( EvalType type, 
                                           PatchData& pd,
                                           double& OF_val,
                                           std::vector<Vector3D>& grad,
                                           MsqHessian& hessian,
                                           MsqError& err )
{
  QualityMetric* qm = get_quality_metric();
  qm->get_evaluations( pd, qmHandles, OF_FREE_EVALS_ONLY, err );  MSQ_ERRFALSE(err);
  double negate_flag = qm->get_negate_flag();
  
    // zero gradient and hessian
  grad.clear();
  grad.resize( pd.num_free_vertices(), 0.0 );
  hessian.zero_out();
  
  double QM_val, QM_pow = 1.0;
  double fac1, fac2;
  Matrix3D elem_outer_product;
  bool qm_bool;
  size_t i, j, n;
  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( pd, *k, QM_val, mIndices, mGradient, mHessian, 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;
      n=0;
      for (i=0; i<nve; ++i) {
        for (j=i; j<nve; ++j) {
            //negate if necessary
          mHessian[n] *= negate_flag;
          hessian.add( mIndices[i], mIndices[j], mHessian[n], err ); MSQ_ERRFALSE(err);
          ++n;
        }
      }
      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();

      n=0;
      for (i=0; i<nve; ++i) {
        for (j=i; j<nve; ++j) {
          elem_outer_product.outer_product(mGradient[i], mGradient[j]);
          elem_outer_product *= fac2;
          mHessian[n] *= fac1;
          mHessian[n] += elem_outer_product;
          mHessian[n] *= negate_flag;
          hessian.add( mIndices[i], mIndices[j], mHessian[n], err ); MSQ_ERRFALSE(err);
          ++n;
        }
      }
    } else {
      MSQ_SETERR(err)(" invalid P value.", MsqError::INVALID_STATE);
      return false;
    }


//.........这里部分代码省略.........

QualityMetric::MetricType NumericQM::get_metric_type() const
  { return realMetric->get_metric_type(); }

/*! \fn LPtoPTemplate::compute_analytical_hessian(PatchData &pd, MsqHessian &hessian, MsqError &err)

    For each element, each entry to be accumulated in the Hessian for
    this objective function (\f$ \sum_{e \in E} Q(e)^p \f$ where \f$ E \f$
    is the set of all elements in the patch) has the form:
    \f$ pQ(e)^{p-1} \nabla^2 Q(e) + p(p-1)Q(e)^{p-2} \nabla Q(e) [\nabla Q(e)]^T \f$.

    For \f$ p=2 \f$, this simplifies to
    \f$ 2Q(e) \nabla^2 Q(e) + 2 \nabla Q(e) [\nabla Q(e)]^T \f$.

    For \f$ p=1 \f$, this simplifies to \f$ \nabla^2 Q(e) \f$.

    The \f$ p=1 \f$ simplified version is implemented directly
    to speed up computation. 

    This function does not support vertex-based metrics.
    
    \param pd The PatchData object for which the objective function
           hessian is computed.
    \param hessian: this object must have been previously initialized.
*/
bool LPtoPTemplate::compute_analytical_hessian(PatchData &pd,
                                               MsqHessian &hessian,
                                               Vector3D *const &grad,
                                               double &OF_val,
                                               MsqError &err)
{
  double scaling_value=1.0;
  
  MSQ_FUNCTION_TIMER( "LPtoPTemplate::compute_analytical_hessian" );

  MsqMeshEntity* elements = pd.get_element_array(err); MSQ_ERRZERO(err);
  MsqVertex* vertices = pd.get_vertex_array(err); MSQ_ERRZERO(err);
  size_t num_elems = pd.num_elements();
    //if scaling divide by the number of elements.
  if(dividingByN){
    if(num_elems<=0) {
      MSQ_SETERR(err)("LPtoP is attempting to divide by zero in analytical Hessian.",
                      MsqError::INVALID_MESH);
      return false;
    }
    scaling_value/=num_elems;
  }
  
  size_t num_vertices = pd.num_vertices();
  Matrix3D elem_hessian[MSQ_MAX_NUM_VERT_PER_ENT*(MSQ_MAX_NUM_VERT_PER_ENT+1)/2];
  Matrix3D elem_outer_product;
  Vector3D grad_vec[MSQ_MAX_NUM_VERT_PER_ENT];
  double QM_val;
  double fac1, fac2;
  Matrix3D grad_outprod;
  bool qm_bool;
  QualityMetric* currentQM = get_quality_metric();
  
  MsqVertex* ele_free_vtces[MSQ_MAX_NUM_VERT_PER_ENT];
  short i;
  for (i=0; i<MSQ_MAX_NUM_VERT_PER_ENT; ++i) ele_free_vtces[i]=NULL;
  
  const size_t* vtx_indices;
    
  size_t e, v;
  size_t nfve; // number of free vertices in element
  short j,n;

  hessian.zero_out();
  for (v=0; v<num_vertices; ++v) grad[v] = 0.;
  OF_val = 0.;
  
  // Loops over all elements in the patch.
  for (e=0; e<num_elems; ++e) {
    short nve = elements[e].vertex_count();
    
    // Gets a list of free vertices in the element.
    vtx_indices = elements[e].get_vertex_index_array();
    nfve=0;
    for (i=0; i<nve; ++i) {
      if ( vertices[vtx_indices[i]].is_free_vertex() ) {
        ele_free_vtces[nfve] = vertices + vtx_indices[i];
        ++nfve;
      }
    }
    
    // Computes \nabla^2 Q(e). Only the free vertices will have non-zero entries. 
    qm_bool = currentQM->compute_element_hessian(pd,
                                    elements+e, ele_free_vtces,
                                    grad_vec, elem_hessian,
                                    nfve, QM_val, err);
    if (MSQ_CHKERR(err) || !qm_bool) return false;


    // **** Computes Hessian ****
    double QM_pow=1.;
    if (pVal == 1) {
      n=0;
      for (i=0; i<nve; ++i) {
        for (j=i; j<nve; ++j) {
            //negate if necessary
          elem_hessian[n] *= (scaling_value * get_negate_flag());
          ++n;
        }
//.........这里部分代码省略.........

bool CompareMetric::evaluate_with_indices( PatchData& pd,
                 size_t handle,
                 double& value,
                 std::vector<size_t>& indices,
                 MsqError& err )
{
  double m2val;
  bool r1, r2;
  m2Handles.clear();
  r1 = metric1->evaluate_with_indices( pd, handle, value, indices, err ); MSQ_ERRZERO(err);
  r2 = metric2->evaluate_with_indices( pd, handle, m2val, m2Handles, 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_indices:\n"
                    "\t%s %f vs. %s %f\n", (unsigned long)handle,
                    r1?"true":"false",value,r2?"true":"false",m2val);
  }
  else {
    bool same = (indices.size() == m2Handles.size());
    std::sort( m2Handles.begin(), m2Handles.end() );
    for (std::vector<size_t>::iterator i = indices.begin(); i != indices.end(); ++i) 
      if (!std::binary_search( m2Handles.begin(), m2Handles.end(), *i ))
        same = false;
    if (!same) {
      MSQ_SETERR(err)(MsqError::INVALID_STATE,
                      "Metrics returned incompatible lists of vertex indices"
                      " for handle %lu in evaluate_with_indices\n.", 
                      (unsigned long)handle );
    }
  }
  
  return r1 && !err;
}

std::string CompareMetric::get_name() const
{
  std::string n = metric1->get_name();
  n += " =? ";
  n += metric2->get_name();
  return n;
}

bool CompareMetric::evaluate_with_gradient( PatchData& pd,
                 size_t handle,
                 double& value,
                 std::vector<size_t>& indices,
                 std::vector<Vector3D>& gradient,
                 MsqError& err )
{
  double m2val;
  bool r1, r2;
  m2Handles.clear();
  m2Grad.clear();
  r1 = metric1->evaluate_with_gradient( pd, handle, value, indices, gradient, err ); MSQ_ERRZERO(err);
  r2 = metric2->evaluate_with_gradient( pd, handle, m2val, m2Handles, m2Grad, 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_gradient:\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;
    int grad_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;
    }
      
    if (!same) {
      MSQ_SETERR(err)(MsqError::INVALID_STATE,
                      "Metrics returned incompatible lists of vertex indices"
                      " for handle %lu in evaluate_with_gradient\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_gradient\n.", 
                      grad_diff, (unsigned)gradient.size(), 
                      (unsigned long)handle );
    }
  }
  
  return r1 && !err;
}

bool CompareMetric::evaluate( PatchData& pd, 
                 size_t handle, 
                 double& value, 
                 MsqError& err )
{
  double m2val;
  bool r1, r2;
  r1 = metric1->evaluate( pd, handle, value, err ); MSQ_ERRZERO(err);
  r2 = metric2->evaluate( pd, handle, m2val, 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:\n"
                    "\t%s %f vs. %s %f\n", (unsigned long)handle,
                    r1?"true":"false",value,r2?"true":"false",m2val);
  }
  
  return r1 && !err;
}

void CompareMetric::get_evaluations( PatchData& pd, 
                                     std::vector<size_t>& handles, 
                                     bool free_vertices_only,
                                     MsqError& err )
{
  if (maskPlane)
    get_mask_axis(pd);

  m2Handles.clear();
  metric1->get_evaluations( pd, handles, free_vertices_only, err ); MSQ_ERRRTN(err);
  metric2->get_evaluations( pd, m2Handles, free_vertices_only, err ); MSQ_ERRRTN(err);
  bool same = (handles.size() == m2Handles.size());
  std::sort( m2Handles.begin(), m2Handles.end() );
  for (std::vector<size_t>::iterator i = handles.begin(); i != handles.end(); ++i) 
    if (!std::binary_search( m2Handles.begin(), m2Handles.end(), *i ))
      same = false;
  if (!same) {
    MSQ_SETERR(err)("Metrics have incompatible lists of evaluation handles.\n",
                    MsqError::INVALID_STATE);
  }
}

bool CompareMetric::evaluate_with_Hessian( PatchData& pd,
                 size_t handle,
                 double& value,
                 std::vector<size_t>& indices,
                 std::vector<Vector3D>& gradient,
                 std::vector<Matrix3D>& Hessian,
                 MsqError& err )
{
  double m2val;
  bool r1, r2;
  m2Handles.clear();
  m2Grad.clear();
  m2Hess.clear();
  r1 = metric1->evaluate_with_Hessian( pd, handle, value, indices, gradient, Hessian, err ); MSQ_ERRZERO(err);
  r2 = metric2->evaluate_with_Hessian( pd, handle, m2val, m2Handles, m2Grad, m2Hess, err ); MSQ_ERRZERO(err);
  if (r1 != r2 || fabs(value - m2val) > epsilon) {
    MSQ_SETERR(err)(MsqError::INVALID_STATE,
                    "Metrics returned different values for handle %lu in "
                    "evaluate_with_Hessian:\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;
    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;
          // call again for so debugger can step into it after failure is found
        std::vector<size_t> i2;
        std::vector<Vector3D> g2;
        std::vector<Matrix3D> h2;
        metric2->evaluate_with_Hessian(pd, handle, m2val, i2, g2, h2, err );
      }
    }  
     
    if (!same) {
      MSQ_SETERR(err)(MsqError::INVALID_STATE,
                      "Metrics returned incompatible lists of vertex indices"
                      " for handle %lu in evaluate_with_Hessian\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\n.", 
                      grad_diff, (unsigned)gradient.size(), 
                      (unsigned long)handle );
    }
    else {
      size_t row, col, row2, col2, idx, idx2;
      for (row = idx = 0; row < indices.size(); ++row) {
        row2 = std::lower_bound( m2Handles.begin(), m2Handles.end(), indices[row] ) - m2Handles.begin();
        for (col = row; col < indices.size(); ++col, ++idx) {
          col2 = std::lower_bound( m2Handles.begin(), m2Handles.end(), indices[col] ) - m2Handles.begin();
          if (row2 <= col2) {
            idx2 = indices.size()*row2 - row2*(row2+1)/2 + col2;
            if (!equal(Hessian[idx], m2Hess[idx2]))
              ++hess_diff;
          }
          else {
            idx2 = indices.size()*col2 - col2*(col2+1)/2 + row2;
            if (!equal(Hessian[idx], transpose(m2Hess[idx2])))
              ++hess_diff;
          }
        }
      }

      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\n.", 
                        hess_diff, (unsigned)Hessian.size(), 
                        (unsigned long)handle );
      }
    }
  }
  
  return r1 && !err;
}

int NumericQM::get_negate_flag() const
{ return realMetric->get_negate_flag(); }