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

//.........这里部分代码省略.........
  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;
      }
      
      // undo previous step : restore vertex coordinates
      pd.set_to_vertices_memento( pd_previous_coords, err );  MSQ_ERRRTN(err);
    }
   
      // Re-evaluate objective function to get gradient.
      // Calling the 'update' form here incorporates the new vertex 
      // positions into the 'accumulated' value if we are doing a 
      // block coordinate descent optimization.
    obj_func.update(pd, original_value, gradient, err ); MSQ_ERRRTN(err);
    if (projectGradient) {
      //if (cosineStep) {
      //  unprojected = gradient;
      //  pd.project_gradient( gradient, err ); MSQ_ERRRTN(err);
      //  double dot = inner_product( arrptr(gradient), arrptr(unprojected), gradient.size() );
      //  double lensqr1 = length_squared( gradient );
      //  double lensqr2 = length_squared( unprojected );
      //  double cossqr = dot * dot / lensqr1 / lensqr2;
      //  step_size *= sqrt(cossqr);
      //}
      //else {
        pd.project_gradient( gradient, err ); MSQ_ERRRTN(err);
      //}      
    }
      
      // Update terination criterion for next iteration.
      // This is necessary for efficiency.  Some values can be adjusted
      // for each iteration so we don't need to re-caculate the value
      // over the entire mesh.
    term_crit->accumulate_patch( pd, err );  MSQ_ERRRTN(err);
    term_crit->accumulate_inner( pd, original_value, arrptr(gradient), err ); MSQ_ERRRTN(err); 
      
      // Calculate initial step size for next iteration using step size 
      // from this iteration
    step_size *= norm_squared;
    norm_squared = length_squared( gradient );
//    if (norm_squared < DBL_EPSILON)
//      break;
  if (norm_squared >= DBL_EPSILON)
    step_size /= norm_squared;
  }
}