C++ MsqHessian类代码示例

/** Internal helper function for test_eval_type */
double ObjectiveFunctionTests::evaluate_internal( 
                                 ObjectiveFunction::EvalType type, 
                                 OFTestMode test_mode,
                                 ObjectiveFunction* of )
{
  MsqPrintError err(cout);
  vector<Vector3D> grad;
  vector<SymMatrix3D> diag;
  MsqHessian hess;
  bool valid = false;
  double result;
  
  switch (test_mode) {
    case EVAL:
      valid = of->evaluate( type, patch(), result, OF_FREE_EVALS_ONLY, err );
      break;
    case GRAD:
      valid = of->evaluate_with_gradient( type, patch(), result, grad, err );
      break;
    case DIAG:
      valid = of->evaluate_with_Hessian_diagonal( type, patch(), result, grad, diag, err );
      break;
    case HESS:
      hess.initialize( patch(), err );
      ASSERT_NO_ERROR( err );
      valid = of->evaluate_with_Hessian( type, patch(), result, grad, hess, err );
      break;
    default:
      CPPUNIT_ASSERT(false);
  }
  
  ASSERT_NO_ERROR( err );
  CPPUNIT_ASSERT(valid);
  return result;
}

void CompositeOFTest::get_hessians( MsqHessian& LP1_hess,
                                    MsqHessian& LP2_hess,
                                    ObjectiveFunction& OF,
                                    MsqHessian& OF_hess )
{
    MsqPrintError err(cout);
    PatchData pd;
    create_twelve_hex_patch( pd, err );
    ASSERT_NO_ERROR( err );

    LP1_hess.initialize( pd, err );
    ASSERT_NO_ERROR(err);
    LP2_hess.initialize( pd, err );
    ASSERT_NO_ERROR(err);
    OF_hess .initialize( pd, err );
    ASSERT_NO_ERROR(err);

    std::vector<Vector3D> grad;
    bool rval;
    double value;
    rval = LP1.evaluate_with_Hessian( ObjectiveFunction::CALCULATE, pd, value, grad, LP1_hess, err );
    ASSERT_NO_ERROR(err);
    CPPUNIT_ASSERT(rval);
    rval = LP2.evaluate_with_Hessian( ObjectiveFunction::CALCULATE, pd, value, grad, LP2_hess, err );
    ASSERT_NO_ERROR(err);
    CPPUNIT_ASSERT(rval);
    rval = OF .evaluate_with_Hessian( ObjectiveFunction::CALCULATE, pd, value, grad, OF_hess , err );
    ASSERT_NO_ERROR(err);
    CPPUNIT_ASSERT(rval);
}

void StdDevTemplateTest::test_hessian_fails_sqr()
{
  MsqError err;
  double value;
  bool rval;
  vector<Vector3D> grad;
  MsqHessian Hess;
  Hess.initialize( patch(), err );
  CPPUNIT_ASSERT(!MSQ_CHKERR(err));
  
  OFTestQM metric( &value, 1 );
  VarianceTemplate func( &metric );
  rval = func.evaluate_with_Hessian( ObjectiveFunction::CALCULATE, patch(), value, grad, Hess, err );
  CPPUNIT_ASSERT(err);
}

bool ObjectiveFunction::evaluate_with_Hessian_diagonal( EvalType type, 
                                        PatchData& pd,
                                        double& value_out,
                                        std::vector<Vector3D>& grad_out,
                                        std::vector<SymMatrix3D>& hess_diag_out,
                                        MsqError& err )
{
  MsqHessian hess;
  hess.initialize( pd, err ); MSQ_ERRZERO(err);
  bool val = evaluate_with_Hessian( type, pd, value_out, grad_out, hess, err );
  MSQ_ERRZERO(err);
  hess_diag_out.resize( hess.size() );
  for (size_t i = 0; i < hess.size(); ++i)
    hess_diag_out[i] = hess.get_block(i,i)->upper();
  return val;
}

bool CompositeOFAdd::evaluate_with_Hessian( EvalType type, 
                                            PatchData& pd,
                                            double& value_out,
                                            std::vector<Vector3D>& grad_out,
                                            MsqHessian& Hessian_out,
                                            MsqError& err )
{
  double value_2;
  bool ok;
  
  mHessian.initialize( Hessian_out );
  
  ok = objFunc1->evaluate_with_Hessian( type, pd, value_out, grad_out, Hessian_out, err );
  if (MSQ_CHKERR(err) || !ok) return false;
  ok = objFunc2->evaluate_with_Hessian( type, pd, value_2, mGradient, mHessian, err );
  if (MSQ_CHKERR(err) || !ok) return false;
  
  value_out += value_2;
  
  assert( grad_out.size() == pd.num_free_vertices() );
  assert( mGradient.size() == pd.num_free_vertices() );
  
  for (size_t i = 0; i < pd.num_free_vertices(); ++i) 
    grad_out[i] += mGradient[i];
  Hessian_out.add( mHessian );
  return true;
}

void CompositeOFTest::test_multiply_hessian()
{
    MsqError err;
    PatchData pd;
    create_twelve_hex_patch( pd, err );
    ASSERT_NO_ERROR( err );

    // this should always fail because the Hessian is not sparse
    CompositeOFMultiply OF( &LP1, &LP2 );
    double value;
    MsqHessian hess;
    hess.initialize( pd, err );
    ASSERT_NO_ERROR(err);
    std::vector<Vector3D> grad;
    OF.evaluate_with_Hessian( ObjectiveFunction::CALCULATE, pd, value, grad, hess, err );
    CPPUNIT_ASSERT(err);
}

void ObjectiveFunctionTest::test_compute_ana_hessian_tet_scaled()
{
    MsqPrintError err(cout);
    PatchData tetPatch;
    create_qm_two_tet_patch(tetPatch,err);
    ASSERT_NO_ERROR(err);

    // creates a mean ratio quality metric ...
    IdealWeightInverseMeanRatio* mean_ratio = new IdealWeightInverseMeanRatio(err);
    CPPUNIT_ASSERT(!err);
    mean_ratio->set_averaging_method(QualityMetric::SUM);

    // ... and builds an objective function with it
    LPtoPTemplate LP2(mean_ratio, 2, err);
    LP2.set_dividing_by_n(true);

    MsqHessian H;
    std::vector<Vector3D> g;
    double dummy;
    H.initialize(tetPatch, err);
    CPPUNIT_ASSERT(!err);
    LP2.evaluate_with_Hessian(ObjectiveFunction::CALCULATE, tetPatch, dummy, g, H, err);
    CPPUNIT_ASSERT(!err);
    CPPUNIT_ASSERT_EQUAL( tetPatch.num_free_vertices(), g.size() );

    Matrix3D mat00(" 2.44444  0.2566   0.181444 "
                   " 0.2566   2.14815  0.104757 "
                   " 0.181444 0.104757 2.07407 ");

    mat00*=.5;

    Matrix3D mat13(" 5.47514 3.16659    9.83479 "
                   " -1.11704 -5.29718 -3.67406 "
                   " 10.3635 -13.5358  -15.5638 ");

    mat13*=.5;

    CPPUNIT_ASSERT_MATRICES_EQUAL( mat00, *H.get_block(0,0), 1e-4 );
    CPPUNIT_ASSERT_MATRICES_EQUAL( mat13, *H.get_block(1,3), 1e-4 );

//    cout << H <<endl;
    delete mean_ratio;
}

void ObjectiveFunctionTests::test_handles_qm_error( 
                                   OFTestMode test_mode,
                                   ObjectiveFunctionTemplate* of )
  
{
  OFTestBadQM metric(true);
  of->set_quality_metric( &metric );
  
  MsqError err;
  vector<Vector3D> grad;
  vector<SymMatrix3D> diag;
  MsqHessian hess;
  double result;
  bool valid;
  
  switch (test_mode) {
    case EVAL:
      valid = of->evaluate( ObjectiveFunction::CALCULATE, patch(), result, OF_FREE_EVALS_ONLY, err );
      break;
    case GRAD:
      valid = of->evaluate_with_gradient( ObjectiveFunction::CALCULATE, patch(), result, grad, err );
      break;
    case DIAG:
      valid = of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, patch(), result, grad, diag, err );
      break;
    case HESS:
      hess.initialize( patch(), err );
      ASSERT_NO_ERROR( err );
      valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, patch(), result, grad, hess, err );
      break;
    default:
      CPPUNIT_ASSERT(false);
  }
  
  CPPUNIT_ASSERT(err);
}

bool PMeanPTemplate::evaluate_with_Hessian( EvalType type, 
                                        PatchData& pd,
                                        double& value_out,
                                        std::vector<Vector3D>& grad_out,
                                        MsqHessian& Hessian_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
  grad_out.clear();
  grad_out.resize( pd.num_free_vertices(), 0.0 );
  Hessian_out.zero_out();
  
    // calculate OF value and gradient for just the patch
  std::vector<size_t>::const_iterator i;
  size_t j, k, n;
  double value, working_sum = 0.0;
  const double f1 = qm->get_negate_flag() * mPower.value();
  const double f2 = f1 * (mPower.value() - 1);
  Matrix3D m;
  for (i = qmHandles.begin(); i != qmHandles.end(); ++i)
  {
    bool result = qm->evaluate_with_Hessian( pd, *i, value, mIndices, mGradient, mHessian, err );
    if (MSQ_CHKERR(err) || !result)
      return false;
    if (fabs(value) < DBL_EPSILON)
      continue;
    
    const size_t nfree = mIndices.size();
    n = 0;
    if (mPower.value() == 1.0) {
      working_sum += mPower.raise( value );
      for (j = 0; j < nfree; ++j) {
        mGradient[j] *= f1;
        grad_out[mIndices[j]] += mGradient[j];
        for (k = j; k < nfree; ++k) {
          mHessian[n] *= f1;
          Hessian_out.add( mIndices[j], mIndices[k], mHessian[n], err );  MSQ_ERRFALSE(err);
          ++n;
        }
     }
    }
    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) {
        for (k = j; k < nfree; ++k) {
          m.outer_product( mGradient[j], mGradient[k] );
          m *= hf;
          mHessian[n] *= gf;
          m += mHessian[n];
          Hessian_out.add( mIndices[j], mIndices[k], m, err );  MSQ_ERRFALSE(err);
          ++n;
        }
      }
      for (j = 0; j < nfree; ++j) {
        mGradient[j] *= gf;
        grad_out[mIndices[j]] += 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;
  std::vector<Vector3D>::iterator g;
  for (g = grad_out.begin(); g != grad_out.end(); ++g)
    *g *= inv_n;
  Hessian_out.scale( inv_n );
  return true;
}

/*\ 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;
    }


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

void ObjectiveFunctionTests::compare_numerical_hessian( ObjectiveFunction* of,
                                                        bool diagonal_only )
{
  const double delta = 0.0001;

  MsqPrintError err(std::cout);
  PatchData pd;
  create_qm_two_tet_patch( pd, err ); 
  ASSERT_NO_ERROR( err );
  CPPUNIT_ASSERT( pd.num_free_vertices() != 0 );
  
    // get analytical Hessian from objective function
  std::vector<Vector3D> grad;
  std::vector<SymMatrix3D> diag;
  MsqHessian hess;
  hess.initialize( pd, err );
  ASSERT_NO_ERROR( err );
  double value;
  bool valid;
  if (diagonal_only)
    valid = of->evaluate_with_Hessian_diagonal( ObjectiveFunction::CALCULATE, pd, value, grad, diag, err );
  else
    valid = of->evaluate_with_Hessian( ObjectiveFunction::CALCULATE, pd, value, grad, hess, err );
  ASSERT_NO_ERROR(err); CPPUNIT_ASSERT(valid);

  
    // do numerical approximation of each block and compare to analytical value
  for (size_t i = 0; i < pd.num_free_vertices(); ++i) {
    const size_t j_end = diagonal_only ? i+1 : pd.num_free_vertices();
    for (size_t j = i; j < j_end; ++j) {
        // do numerical approximation for block corresponding to
        // coorindates for ith and jth vertices.
      Matrix3D block;    
      for (int k = 0; k < 3; ++k) {
        for (int m = 0; m < 3; ++m) {
          double dk, dm, dkm;
          Vector3D ik = pd.vertex_by_index(i);
          Vector3D im = pd.vertex_by_index(j);
          
          Vector3D delta_k(0.0); delta_k[k] = delta;
          pd.move_vertex( delta_k, i, err ); ASSERT_NO_ERROR(err);
          valid = of->evaluate( ObjectiveFunction::CALCULATE, pd, dk, true, err );
          ASSERT_NO_ERROR(err); CPPUNIT_ASSERT(valid);
          
          Vector3D delta_m(0.0); delta_m[m] = delta;
          pd.move_vertex( delta_m, j, err ); ASSERT_NO_ERROR(err);
          valid = of->evaluate( ObjectiveFunction::CALCULATE, pd, dkm, true, err );
          ASSERT_NO_ERROR(err); CPPUNIT_ASSERT(valid);
          
            // be careful here that we do the right thing if i==j
          pd.set_vertex_coordinates( ik, i, err ); ASSERT_NO_ERROR(err);
          pd.set_vertex_coordinates( im, j, err ); ASSERT_NO_ERROR(err);
          pd.move_vertex( delta_m, j, err ); ASSERT_NO_ERROR(err);
          valid = of->evaluate( ObjectiveFunction::CALCULATE, pd, dm, true, err );
          ASSERT_NO_ERROR(err); CPPUNIT_ASSERT(valid);
          
          pd.set_vertex_coordinates( ik, i, err ); ASSERT_NO_ERROR(err);
          pd.set_vertex_coordinates( im, j, err ); ASSERT_NO_ERROR(err);
          
          block[k][m] = (dkm - dk - dm + value)/(delta*delta);
        }
      }
        // compare to analytical value
      if (diagonal_only) {
        CPPUNIT_ASSERT(i == j); // see j_end above
        CPPUNIT_ASSERT(i < diag.size());
        CHECK_EQUAL_MATRICES( block, Matrix3D(diag[i]) );
      }
      else {
        Matrix3D* m = hess.get_block( i, j );
        Matrix3D* mt = hess.get_block( j, i );
        if (NULL != m) {
          CHECK_EQUAL_MATRICES( block, *m );
        }
        if (NULL != mt) {
          CHECK_EQUAL_MATRICES( transpose(block), *m );
        }
        if (NULL == mt && NULL == m) {
          CHECK_EQUAL_MATRICES( Matrix3D(0.0), block );
        }
      }
    }
  }
}

/*! \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;
        }
//.........这里部分代码省略.........