C++ PatchData::move_free_vertices_constrained方法代码示例

/*!Performs Conjugate gradient minimization on the PatchData, pd.*/
void ConjugateGradient::optimize_vertex_positions(PatchData &pd, 
                                                MsqError &err){
  // pd.reorder();

  MSQ_FUNCTION_TIMER( "ConjugateGradient::optimize_vertex_positions" );

  Timer c_timer;
  size_t num_vert=pd.num_free_vertices();
  if(num_vert<1){
     MSQ_DBGOUT(1) << "\nEmpty free vertex list in ConjugateGradient\n";
     return;
  }
/*  
    //zero out arrays
  int zero_loop=0;
  while(zero_loop<arraySize){
    fGrad[zero_loop].set(0,0,0);
    pGrad[zero_loop].set(0,0,0);
    fNewGrad[zero_loop].set(0,0,0);
    ++zero_loop;
  }
*/
  
  // get OF evaluator
  OFEvaluator& objFunc = get_objective_function_evaluator();
  
  size_t ind;
    //Michael cull list:  possibly set soft_fixed flags here
  
   //MsqFreeVertexIndexIterator free_iter(pd, err);  MSQ_ERRRTN(err);
  
      
   double f=0;
   //Michael, this isn't equivalent to CUBIT because we only want to check
   //the objective function value of the 'bad' elements
   //if invalid initial patch set an error.
   bool temp_bool = objFunc.update(pd, f, fGrad, err);
   assert(fGrad.size() == num_vert);
   if(MSQ_CHKERR(err))
      return;
  if( ! temp_bool){
    MSQ_SETERR(err)("Conjugate Gradient not able to get valid gradient "
                    "and function values on intial patch.", 
                    MsqError::INVALID_MESH);
    return;
  }
  double grad_norm=MSQ_MAX_CAP;
  
  if(conjGradDebug>0){
    MSQ_PRINT(2)("\nCG's DEGUB LEVEL = %i \n",conjGradDebug);
    grad_norm=Linf(arrptr(fGrad),fGrad.size());
    MSQ_PRINT(2)("\nCG's FIRST VALUE = %f,grad_norm = %f",f,grad_norm);
    MSQ_PRINT(2)("\n   TIME %f",c_timer.since_birth());
    grad_norm=MSQ_MAX_CAP;
  }

    //Initializing pGrad (search direction).  
  pGrad.resize(fGrad.size());
  for (ind = 0; ind < num_vert; ++ind)
    pGrad[ind]=(-fGrad[ind]);

  int j=0; // total nb of step size changes ... not used much
  int i=0; // iteration counter
  unsigned m=0; // 
  double alp=MSQ_MAX_CAP; // alp: scale factor of search direction
    //we know inner_criterion is false because it was checked in
    //loop_over_mesh before being sent here.
  TerminationCriterion* term_crit=get_inner_termination_criterion();
  
    //while ((i<maxIteration && alp>stepBound && grad_norm>normGradientBound)
    //     && !inner_criterion){
  while(!term_crit->terminate()){
    ++i;
      //std::cout<<"\Michael delete i = "<<i;
    int k=0;
    alp=get_step(pd,f,k,err);
    j+=k;
    if(conjGradDebug>2){
      MSQ_PRINT(2)("\n  Alp initial, alp = %20.18f",alp);
    }

     // if alp == 0, revert to steepest descent search direction
    if(alp==0){
      for (m = 0; m < num_vert; ++m) {
        pGrad[m]=(-fGrad[m]);
      }
      alp=get_step(pd,f,k,err);
      j+=k;
      if(conjGradDebug>1){
        MSQ_PRINT(2)("\n CG's search direction reset.");
        if(conjGradDebug>2)
          MSQ_PRINT(2)("\n  Alp was zero, alp = %20.18f",alp);
      }
      
    }
    if(alp!=0){
      pd.move_free_vertices_constrained( arrptr(pGrad), num_vert, alp, err );
      MSQ_ERRRTN(err);
      
//.........这里部分代码省略.........

void QuasiNewton::optimize_vertex_positions( PatchData& pd, MsqError& err )
{
  TerminationCriterion& term = *get_inner_termination_criterion();
  OFEvaluator& func = get_objective_function_evaluator();
  
  const double sigma = 1e-4;
  const double beta0 = 0.25;
  const double beta1 = 0.80;
  const double tol1 = 1e-8;
  const double epsilon = 1e-10;

  double norm_r; //, norm_g;
  double alpha, beta;
  double obj, objn;

  size_t i;
  
    // Initialize stuff
  const size_t nn = pd.num_free_vertices();
  double a[QNVEC], b[QNVEC], r[QNVEC];
  for (i = 0; i < QNVEC; ++i)
    r[i] = 0;
  for (i = 0; i <= QNVEC; ++i) {
    v[i].clear();
    v[i].resize( nn, Vector3D(0.0) );
    w[i].clear();
    w[i].resize( nn, Vector3D(0.0) );
  }
  d.resize( nn );
  mHess.resize( nn );  //hMesh(mesh);

  bool valid = func.update( pd, obj, v[QNVEC], mHess, err ); MSQ_ERRRTN(err);
  if (!valid) {
    MSQ_SETERR(err)("Initial objective function is not valid", MsqError::INVALID_MESH);
    return;
  }

  while (!term.terminate()) {
    pd.recreate_vertices_memento( mMemento, err ); MSQ_ERRRTN(err);
    pd.get_free_vertex_coordinates( w[QNVEC] );

    x = v[QNVEC];
    for (i = QNVEC; i--; ) {
      a[i] = r[i] * inner( &(w[i][0]), arrptr(x), nn );
      plus_eq_scaled( arrptr(x), -a[i], &v[i][0], nn );
    }
     
    solve( arrptr(d), arrptr(x) );
  
    for (i = QNVEC; i--; ) {
      b[i] = r[i] * inner( &(v[i][0]), arrptr(d), nn );
      plus_eq_scaled( arrptr(d), a[i]-b[i], &(w[i][0]), nn );
    }
    
    alpha = -inner( &(v[QNVEC][0]), arrptr(d), nn );  /* direction is negated */
    if (alpha > 0.0) {
      MSQ_SETERR(err)("No descent.", MsqError::INVALID_MESH);
      return;
    }
   
    alpha *= sigma;
    beta = 1.0;
    
    pd.move_free_vertices_constrained( arrptr(d), nn, -beta, err ); MSQ_ERRRTN(err);
    valid = func.evaluate( pd, objn, v[QNVEC], err ); 
    if (err.error_code() == err.BARRIER_VIOLATED)             
      err.clear();  // barrier violated does not represent an actual error here
    MSQ_ERRRTN(err);
    if (!valid ||
        (obj - objn < -alpha*beta - epsilon &&
         length( &(v[QNVEC][0]), nn ) >= tol1)) {
      
      if (!valid)  // function not defined at trial point
        beta *= beta0;
      else  // unacceptable iterate
        beta *= beta1;
      
      for (;;) {
        if (beta < tol1) {
          pd.set_to_vertices_memento( mMemento, err ); MSQ_ERRRTN(err);
          MSQ_SETERR(err)("Newton step not good", MsqError::INTERNAL_ERROR);
          return;
        }
      
        pd.set_free_vertices_constrained( mMemento, arrptr(d), nn, -beta, err ); MSQ_ERRRTN(err);
        valid = func.evaluate( pd, objn, err );
        if (err.error_code() == err.BARRIER_VIOLATED)             
          err.clear();  // barrier violated does not represent an actual error here
        MSQ_ERRRTN(err);
        if (!valid) // function undefined at trial point
          beta *= beta0;
        else if (obj - objn < -alpha*beta - epsilon) // unacceptlable iterate
          beta *= beta1;
        else
          break;
      }
    }
    
    for (i = 0; i < QNVEC-1; ++i) {
      r[i] = r[i+1];
//.........这里部分代码省略.........

void SteepestDescent::optimize_vertex_positions(PatchData &pd, 
                                                MsqError &err)
{
  MSQ_FUNCTION_TIMER( "SteepestDescent::optimize_vertex_positions" );

  const int SEARCH_MAX = 100;
  const double c1 = 1e-4;
  //std::vector<Vector3D> unprojected(pd.num_free_vertices()); 
  std::vector<Vector3D> gradient(pd.num_free_vertices()); 
  bool feasible=true;//bool for OF values
  double min_edge_len, max_edge_len;
  double step_size=0, original_value=0, new_value=0;
  double norm_squared=0;
  PatchDataVerticesMemento* pd_previous_coords;
  TerminationCriterion* term_crit=get_inner_termination_criterion();
  OFEvaluator& obj_func = get_objective_function_evaluator();
  
    // get vertex memento so we can restore vertex coordinates for bad steps.
  pd_previous_coords = pd.create_vertices_memento( err ); MSQ_ERRRTN(err);
    // use auto_ptr to automatically delete memento when we exit this function
  std::auto_ptr<PatchDataVerticesMemento> memento_deleter( pd_previous_coords );

    // Evaluate objective function.
    //
    // Always use 'update' version when beginning optimization so that
    // if doing block coordinate descent the OF code knows the set of
    // vertices we are modifying during the optimziation (the subset
    // of the mesh contained in the current patch.)  This has to be
    // done up-front because typically an OF will just store the portion
    // of the OF value (e.g. the numeric contribution to the sum for an
    // averaging OF) for the initial patch.
  feasible = obj_func.update( pd, original_value, gradient, err ); MSQ_ERRRTN(err);
    // calculate gradient dotted with itself
  norm_squared = length_squared( gradient );
  
    //set an error if initial patch is invalid.
  if(!feasible){
    MSQ_SETERR(err)("SteepestDescent passed invalid initial patch.",
                    MsqError::INVALID_ARG);
    return;
  }

    // use edge length as an initial guess for for step size
  pd.get_minmax_edge_length( min_edge_len, max_edge_len );
  //step_size = max_edge_len / std::sqrt(norm_squared);
  //if (!finite(step_size))  // zero-length gradient
  //  return;
//  if (norm_squared < DBL_EPSILON)
//    return;
  if (norm_squared >= DBL_EPSILON)
    step_size = max_edge_len / std::sqrt(norm_squared) * pd.num_free_vertices();

    // The steepest descent loop...
    // We loop until the user-specified termination criteria are met.
  while (!term_crit->terminate()) {
    MSQ_DBGOUT(3) << "Iteration " << term_crit->get_iteration_count() << std::endl;
    MSQ_DBGOUT(3) << "  o  original_value: " << original_value << std::endl;
    MSQ_DBGOUT(3) << "  o  grad norm suqared: " << norm_squared << std::endl;

      // Save current vertex coords so that they can be restored if
      // the step was bad.
    pd.recreate_vertices_memento( pd_previous_coords, err ); MSQ_ERRRTN(err);

      // Reduce step size until it satisfies Armijo condition
    int counter = 0;
    for (;;) {
      if (++counter > SEARCH_MAX || step_size < DBL_EPSILON) {
        MSQ_DBGOUT(3) << "    o  No valid step found.  Giving Up." << std::endl;
        return;
      }
      
      // Move vertices to new positions.
      // Note: step direction is -gradient so we pass +gradient and 
      //       -step_size to achieve the same thing.
      pd.move_free_vertices_constrained( arrptr(gradient), gradient.size(), -step_size, err ); MSQ_ERRRTN(err);
      // Evaluate objective function for new vertices.  We call the
      // 'evaluate' form here because we aren't sure yet if we want to
      // keep these vertices.  Until we call 'update', we have the option
      // of reverting a block coordinate decent objective function's state
      // to that of the initial vertex coordinates.  However, for block
      // coordinate decent to work correctly, we will need to call an
      // 'update' form if we decide to keep the new vertex coordinates.
      feasible = obj_func.evaluate( pd, new_value, err ); 
      if (err.error_code() == err.BARRIER_VIOLATED) 
        err.clear();  // barrier violated does not represent an actual error here
      MSQ_ERRRTN(err);
      MSQ_DBGOUT(3) << "    o  step_size: " << step_size << std::endl;
      MSQ_DBGOUT(3) << "    o  new_value: " << new_value << std::endl;

      if (!feasible) {
        // OF value is invalid, decrease step_size a lot
        step_size *= 0.2;
      }
      else if (new_value > original_value - c1 * step_size * norm_squared) {
        // Armijo condition not met.
        step_size *= 0.5;
      }
      else {
        // Armijo condition met, stop
        break;
//.........这里部分代码省略.........