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

void TerminationCriterion::reset_patch(PatchData &pd, MsqError &err)
{
  const unsigned long totalFlag = terminationCriterionFlag | cullingMethodFlag;
  if (totalFlag & MOVEMENT_FLAGS)
  {
    if (previousVerticesMemento)
      pd.recreate_vertices_memento(previousVerticesMemento,err); 
    else
      previousVerticesMemento = pd.create_vertices_memento(err);
    MSQ_ERRRTN(err);
  }

  if (totalFlag & UNTANGLED_MESH) {
    patchInvertedCount = count_inverted( pd, err );
    //MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o Num Patch Inverted: " << " " << patchInvertedCount << std::endl;
    MSQ_ERRRTN(err);
  }
}

void TerminationCriterion::accumulate_patch( PatchData& pd, MsqError& err )
{
  if (terminationCriterionFlag & MOVEMENT_FLAGS)
  {
    double patch_max_dist = pd.get_max_vertex_movement_squared( previousVerticesMemento, err );
    if (patch_max_dist > maxSquaredMovement)
      maxSquaredMovement = patch_max_dist;
    pd.recreate_vertices_memento( previousVerticesMemento, err );  MSQ_ERRRTN(err);
  }
    
    //if terminating on bounded vertex movement (a bounding box for the mesh)
  if(terminationCriterionFlag & BOUNDED_VERTEX_MOVEMENT)
  {
    const MsqVertex* vert = pd.get_vertex_array(err);
    int num_vert = pd.num_free_vertices();
    int i=0;
      //for each vertex
    for(i=0;i<num_vert;++i)
    {
        //if any of the coordinates are greater than eps
      if( (vert[i][0]>boundedVertexMovementEps) ||
          (vert[i][1]>boundedVertexMovementEps) ||
          (vert[i][2]>boundedVertexMovementEps) )
      {
        ++vertexMovementExceedsBound;
      }
    }
  }


  if ((terminationCriterionFlag|cullingMethodFlag) & UNTANGLED_MESH) {
    size_t new_count = count_inverted( pd, err );
      // be careful here because size_t is unsigned
    globalInvertedCount += new_count;
    globalInvertedCount -= patchInvertedCount;
    patchInvertedCount = new_count;
    //if (innerOuterType==TYPE_OUTER) 
    //  MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o Num Patch Inverted: " << " " << patchInvertedCount << " globalInvertedCount= " << globalInvertedCount << std::endl;
      
    MSQ_ERRRTN(err);
  }
}

double ConjugateGradient::get_step(PatchData &pd,double f0,int &j,
                                   MsqError &err)
{
  // get OF evaluator
  OFEvaluator& objFunc = get_objective_function_evaluator();

  size_t num_vertices=pd.num_free_vertices();
    //initial guess for alp
  double alp=1.0;
  int jmax=100;
  double rho=0.5;
    //feasible=false implies the mesh is not in the feasible region
  bool feasible=false;
  int found=0;
    //f and fnew hold the objective function value
  double f=0;
  double fnew=0;
    //Counter to avoid infinitly scaling alp
  j=0;
  //save memento
  pd.recreate_vertices_memento(pMemento, err);
    //if we must check feasiblility
    //while step takes mesh into infeasible region and ...
  while (j<jmax && !feasible && alp>MSQ_MIN) {
    ++j;
    pd.set_free_vertices_constrained(pMemento,arrptr(pGrad),num_vertices,alp,err);
    feasible=objFunc.evaluate(pd,f,err); MSQ_ERRZERO(err);
      //if not feasible, try a smaller alp (take smaller step)
    if(!feasible){
      alp*=rho;
    }
  }//end while ...
  
    //if above while ended due to j>=jmax, no valid step was found.
  if(j>=jmax){
    MSQ_PRINT(2)("\nFeasible Point Not Found");
    return 0.0;
  }
    //Message::print_info("\nOriginal f %f, first new f = %f, alp = %f",f0,f,alp);
    //if new f is larger than original, our step was too large
  if(f>=f0){
    j=0;
    while (j<jmax && found == 0){
      ++j;
      alp *= rho;
      pd.set_free_vertices_constrained(pMemento,arrptr(pGrad),num_vertices,alp,err);
        //Get new obj value
        //if patch is now invalid, then the feasible region is  convex or
        //we have an error.  For now, we assume an error.
      if(! objFunc.evaluate(pd,f,err) ){
        MSQ_SETERR(err)("Non-convex feasiblility region found.",MsqError::INVALID_MESH);
      }
      pd.set_to_vertices_memento(pMemento,err);MSQ_ERRZERO(err);
        //if our step has now improved the objective function value
      if(f<f0){
        found=1;
      }
    }//   end while j less than jmax
      //Message::print_info("\nj = %d found = %d f = %20.18f f0 = %20.18f\n",j,found,f,f0);
      //if above ended because of j>=jmax, take no step
    if(found==0){
        //Message::print_info("alp = %10.8f, but returning zero\n",alp);
      alp=0.0; 
      return alp;
    }

    j=0;
      //while shrinking the step improves the objFunc value further,
      //scale alp down.  Return alp, when scaling once more would
      //no longer improve the objFunc value.  
    while(j<jmax){
      ++j;
      alp*=rho;
      //step alp in search direction from original positions
      pd.set_free_vertices_constrained(pMemento,arrptr(pGrad),num_vertices,alp,err);MSQ_ERRZERO(err);

        //get new objective function value
      if (! objFunc.evaluate(pd,fnew,err))
        MSQ_SETERR(err)("Non-convex feasiblility region found while "
                        "computing new f.",MsqError::INVALID_MESH);
      if(fnew<f){
        f=fnew;
      }
      else{
	//Reset the vertices to original position
	pd.set_to_vertices_memento(pMemento,err);MSQ_ERRZERO(err);
	alp/=rho;
	return alp;
      }
    }
    //Reset the vertices to original position and return alp
    pd.set_to_vertices_memento(pMemento,err);MSQ_ERRZERO(err);
    return alp;
  }
    //else our new f was already smaller than our original
  else{
    j=0;
      //check to see how large of step we can take
    while (j<jmax && found == 0) {
      ++j;
//.........这里部分代码省略.........

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

//.........这里部分代码省略.........
      //         method though, unless the preconditioner is not positive
      //         definite.
      // If direction is positive, does a gradient (steepest descent) step.

      if (alpha > -epsilon) {
    
        MSQ_DBGOUT(3) << "  o  alpha = " << alpha << " (rejected)" << std::endl;
    
        if (!havePrintedDirectionMessage) {
          MSQ_PRINT(1)("Newton direction not guaranteed descent.  Ensure preconditioner is positive definite.\n");
          havePrintedDirectionMessage = true;
        }

        // TODD: removed performing gradient step here since we will use
        // gradient if step does not produce descent.  Instead we set
        // alpha to a small negative value.

        alpha = -epsilon;

        // alpha = inner(grad, grad, nv); // compute norm squared of gradient
        // if (alpha < 1) alpha = 1;	// take max with constant
        // for (i = 0; i < nv; ++i) {
        //   d[i] = -grad[i] / alpha; 	// compute scaled gradient
        // }
        // alpha = inner(grad, d, nv);  	// recompute alpha
        // 				// equal to one for large gradient
      }
      else {
        MSQ_DBGOUT(3) << "  o  alpha = " << alpha << std::endl;
      }
    
      alpha *= sigma;
      beta = 1.0;
      pd.recreate_vertices_memento(coordsMem, err); MSQ_ERRRTN(err);
    
      // TODD: Unrolling the linesearch loop.  We do a function and
      // gradient evaluation when beta = 1.  Otherwise, we end up
      // in the linesearch regime.  We expect that several
      // evaluations will be made, so we only do a function evaluation
      // and finish with a gradient evaluation.  When beta = 1, we also
      // check the gradient for stability.

      // TODD -- the Armijo linesearch is based on the objective function,
      //         so theoretically we only need to evaluate the objective
      //         function.  However, near a very accurate solution, say with
      //         the two norm of the gradient of the objective function less
      //         than 1e-5, the numerical error in the objective function
      //         calculation is enough that the Armijo linesearch will
      //         fail.  To correct this situation, the iterate is accepted
      //         when the norm of the gradient is also small.  If you need
      //         high accuracy and have a large mesh, talk with Todd about
      //         the numerical issues so that we can fix it.

      // TODD -- the Armijo linesearch here is only good for differentiable
      //         functions.  When projections are involved, you should change
      //	       to a form of the linesearch meant for nondifferentiable
      //         functions.

      pd.move_free_vertices_constrained(arrptr(d), nv, beta, err); MSQ_ERRRTN(err);
      fn_bool = objFunc.evaluate(pd, new_value, grad, err); MSQ_ERRRTN(err);
      if ((fn_bool && (original_value - new_value >= -alpha*beta - epsilon)) ||
          (fn_bool && (length(arrptr(grad), nv) < 100*convTol))) {
        // Armijo linesearch rules passed.
        MSQ_DBGOUT(3) << "  o  beta = " << beta << " (accepted without line search)" << std::endl;
      }
      else {

//.........这里部分代码省略.........
      //get the gradient norms
    if (totalFlag & (GRADIENT_INF_NORM_ABSOLUTE|GRADIENT_INF_NORM_RELATIVE))
    {
      currentGradInfNorm = initialGradInfNorm = Linf(mGrad);
      MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o Initial gradient Inf norm: " << " "
                                     << RPM(initialGradInfNorm) << std::endl;
    }  
      
    if (totalFlag & (GRADIENT_L2_NORM_ABSOLUTE|GRADIENT_L2_NORM_RELATIVE))
    {
      currentGradL2NormSquared = initialGradL2NormSquared = length_squared(mGrad);
      MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o Initial gradient L2 norm: " << " "
                                     << RPM(std::sqrt(initialGradL2NormSquared)) << std::endl;
    }  

      //the OFvalue comes for free, so save it
    previousOFValue=currentOFValue;
    initialOFValue=currentOFValue;
  }
  //find the initial objective function value if needed and not already
  //computed.  If we needed the gradient, we have the OF value for free.
  // Also, if possible, get initial OF value if writing plot file.  Solvers
  // often supply the OF value for subsequent iterations so by calculating
  // the initial value we can generate OF value plots.
  else if ((totalFlag & OF_FLAGS) || 
           (plotFile.is_open() && pd.num_free_vertices() && obj_eval.have_objective_function()))
  {
      //ensure the obj_ptr is not null
    if(!obj_eval.have_objective_function()){
      MSQ_SETERR(err)("Error termination criteria set which uses objective "
                      "functions, but no objective function is available.",
                      MsqError::INVALID_STATE);
      return;
    }
    
    bool b = obj_eval.evaluate(pd, currentOFValue, err); MSQ_ERRRTN(err);
    if (!b){
      MSQ_SETERR(err)("Initial patch is invalid for evaluation.",MsqError::INVALID_STATE);
      return;
    }
      //std::cout<<"\nReseting initial of value = "<<initialOFValue;
    previousOFValue=currentOFValue;
    initialOFValue=currentOFValue;
  }
  
  if (totalFlag & (GRAD_FLAGS|OF_FLAGS))
    MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o Initial OF value: " << " " << RPM(initialOFValue) << std::endl;
  
    // Store current vertex locations now, because we'll
    // need them later to compare the current movement with.
  if (totalFlag & VERTEX_MOVEMENT_RELATIVE)
  {
    if (initialVerticesMemento)
    {
      pd.recreate_vertices_memento( initialVerticesMemento, err );
    }
    else
    {
      initialVerticesMemento = pd.create_vertices_memento( err );
    }
    MSQ_ERRRTN(err);
    maxSquaredInitialMovement = DBL_MAX;
  }
  else {
    maxSquaredInitialMovement = 0;
  }
  
  if (terminationCriterionFlag & UNTANGLED_MESH) {
    globalInvertedCount = count_inverted( pd, err );
    //if (innerOuterType==TYPE_OUTER) MSQ_DBGOUT_P0_ONLY(debugLevel) << par_string() << "  o Num Inverted: " << " " << globalInvertedCount << std::endl;
    patchInvertedCount = 0;
    MSQ_ERRRTN(err);
  }

  if (timeStepFileType) {
      // If didn't already calculate gradient abive, calculate it now.
    if (!(totalFlag & GRAD_FLAGS)) {
      mGrad.resize( pd.num_free_vertices() );
      obj_eval.evaluate(pd, currentOFValue, mGrad, err);
      err.clear();
    }
    write_timestep( pd, mGrad.empty() ? 0 : arrptr(mGrad), err);
  }
    
  if (plotFile.is_open()) {
      // two newlines so GNU plot knows that we are starting a new data set
    plotFile << std::endl << std::endl;
      // write column headings as comment in data file
    plotFile << "#Iter\tCPU\tObjFunc\tGradL2\tGradInf\tMovement\tInverted" << std::endl;
      // write initial values
    plotFile << 0 
     << '\t' << mTimer.since_birth() 
     << '\t' << initialOFValue 
     << '\t' << std::sqrt( currentGradL2NormSquared ) 
     << '\t' << currentGradInfNorm 
     << '\t' << 0.0
     << '\t' << globalInvertedCount
     << std::endl;
  }
}

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