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

void UnOptimizer::optimize_vertex_positions( PatchData &pd, MsqError &err) {
  assert( pd.num_free_vertices() == 1 && pd.vertex_by_index(0).is_free_vertex() );
  std::vector<Vector3D> grad(1);
  double val, junk, coeff;
  bool state;
  
  state = objectiveFunction->evaluate_with_gradient( ObjectiveFunction::CALCULATE,
                                                     pd, val, grad, err );
  MSQ_ERRRTN(err);
  if (!state) {
    MSQ_SETERR(err)(MsqError::INVALID_MESH);
    return;
  }
  grad[0] /= grad[0].length();
  
  PatchDataVerticesMemento* memento = pd.create_vertices_memento( err ); MSQ_ERRRTN(err);
  std::auto_ptr<PatchDataVerticesMemento> deleter( memento );
  pd.get_minmax_edge_length( junk, coeff );
  
  for (int i = 0; i < 100; ++i) {
    pd.set_free_vertices_constrained( memento, &grad[0], 1, coeff, err ); MSQ_ERRRTN(err);
    state = objectiveFunction->evaluate( ObjectiveFunction::CALCULATE, pd, val, true, err );
    MSQ_ERRRTN(err);
    if (state)
      break;
    coeff *= 0.5;
  }
  if (!state) {
    pd.set_to_vertices_memento( memento, err );
  }
}

void NonGradient::initialize_mesh_iteration(PatchData &pd, MsqError &err)
{
  int elementDimension = getPatchDimension( pd, err );  // to do: react to error
  int dimension = elementDimension * pd.num_free_vertices();
  //printPatch( pd, err );
  setDimension(dimension);
  int maxNumEval = 100*dimension;  // 1. Make this a user parameter
  setMaxNumEval(maxNumEval);
  double threshold = 1.e-10; // avoid division by zero
  setThreshold(threshold);
  double minEdgeLen = 0.0;
  double maxEdgeLen = 0.0;
//  double ftol = 0.;
  if( dimension > 0 )
  {
    pd.get_minmax_edge_length( minEdgeLen, maxEdgeLen );
    //ftol = minEdgeLen * 1.e-4; // Turn off Amoeba convergence criterion
    if( mNonGradDebug >= 1 ) 
    {      
         std::cout << "minimum edge length " << minEdgeLen << " maximum edge length " << maxEdgeLen << std::endl;
    }      
    MSQ_PRINT(3)("minimum edge length %e    maximum edge length %e\n", minEdgeLen,  maxEdgeLen);
  }
//  setTolerance(ftol);
  int numRow = dimension;
  int numCol = numRow+1;  
  if( numRow*numCol <= simplex.max_size() )
  { 
    simplex.assign(numRow*numCol, 0.);  // guard against previous simplex value
    double scale = minEdgeLen * mScaleDiameter;; 
    const MsqVertex* coord = pd.get_vertex_array(err);
    if( pd.num_free_vertices() > 1 )
    {
      MSQ_SETERR(err)("Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED);
    }
    size_t index = 0;
    for( int col = 0; col < numCol; col++ )
    {
      for (int row=0;row<numRow;row++)
      {
        simplex[ row + col*numRow ] = coord[index][row];
        if( row == col-1 )
        {
          simplex[ row + col*numRow ] += scale/ static_cast<double>(numCol);
        }
      }
    }
  }
  else
  {
    MSQ_SETERR(err)("Use patch with fewer free vertices", MsqError::OUT_OF_MEMORY);
    if( mNonGradDebug >= 1 ) 
    {      
      std::cout << "ERROR: Too many free vertices in patch" << std::endl;
    }      
    //MSQ_PRINT(1)("ERROR: Too many free vertices in patch\n");
  }
}

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