C++ MsqError::error_code方法代码示例

void TagVertexMesh::initialize( Mesh* mesh, std::string name, MsqError& err )
{
  MeshDecorator::set_mesh( mesh );
  tagName = name;

  tagHandle = get_mesh()->tag_get( tagName, err );
    // If tag isn't defined yet, we're done for now.
  if (err.error_code() == MsqError::TAG_NOT_FOUND) {
    err.clear();
    return;
  } 
  else if (MSQ_CHKERR(err))
    return;
  
    // If tag is already defined, make sure it is the correct type.
  std::string t_name;
  Mesh::TagType type;
  unsigned length;
  tag_properties( tagHandle, t_name, type, length, err ); 
  MSQ_ERRRTN(err);
  if (!(type == Mesh::DOUBLE && length == 3) &&
      !(type == Mesh::BYTE && length == 3*sizeof(double)))
    MSQ_SETERR(err)(MsqError::TAG_ALREADY_EXISTS,
                    "Tag \"%s\" has invalid type or size.",
                    tagName.c_str());
 
    // If tag is already defined and init was true, reset tag
    // values.
  haveTagHandle = true;
}

void CompositeOFTest::test_eval_fails( ObjectiveFunction& OF )
{
    MsqError err;
    double value;
    OF.evaluate( ObjectiveFunction::CALCULATE, patch(), value, false, err );
    CPPUNIT_ASSERT_EQUAL(MsqError::INTERNAL_ERROR, err.error_code());
}

double 
NonGradient::evaluate( double *point,  PatchData &pd, MsqError &err )
{
  if( pd.num_free_vertices() > 1 )
  {
    MSQ_SETERR(err)("Only one free vertex per patch implemented", MsqError::NOT_IMPLEMENTED);
  }
  const size_t vertexIndex = 0; 
  Vector3D originalVec = pd.vertex_by_index(vertexIndex);
  Vector3D pointVec;
  for( int index = 0; index<3; index++)
  {
    pointVec[index] = point[index];
  }
  pd.set_vertex_coordinates( pointVec, vertexIndex, err ); 
  pd.snap_vertex_to_domain( vertexIndex, err );  

  OFEvaluator& obj_func = get_objective_function_evaluator();
  TerminationCriterion* term_crit=get_inner_termination_criterion();

  double value;
  bool feasible = obj_func.evaluate( pd, value, err ); // MSQ_ERRZERO(err);
  term_crit->accumulate_inner( pd, value, 0, err );  //MSQ_CHKERR(err);
  if (err.error_code() == err.BARRIER_VIOLATED)
    err.clear();   // barrier violation not an error here
  MSQ_ERRZERO(err);
  pd.set_vertex_coordinates( originalVec, vertexIndex, err ); 
  if( !feasible && mUseExactPenaltyFunction )  
  { // "value" undefined btw
    double ensureFiniteRtol= .25;
    value = DBL_MAX * ensureFiniteRtol;
    //std::cout << " Infeasible patch: " << value << std::endl; printPatch( pd, err );
  }
  return value;
}

TagHandle TargetWriter::get_tag_handle( const std::string& base_name,
                                        unsigned num_dbl, Mesh* mesh, MsqError& err )
{
  std::ostringstream sstr;
  sstr << base_name << num_dbl;
  
  TagHandle handle = mesh->tag_get( sstr.str().c_str(), err );
  if (!MSQ_CHKERR(err))
  {
    std::string temp_name;
    Mesh::TagType temp_type;
    unsigned temp_length;
    mesh->tag_properties( handle, temp_name, temp_type, temp_length, err );
    MSQ_ERRZERO(err);
    
    if (temp_type != Mesh::DOUBLE || temp_length != num_dbl)
    {
      MSQ_SETERR(err)( MsqError::TAG_ALREADY_EXISTS,
                      "Mismatched type or length for existing tag \"%s\"",
                       sstr.str().c_str() );
    }
  }
  else if (err.error_code() == MsqError::TAG_NOT_FOUND)
  {
    err.clear();
    handle = mesh->tag_create( sstr.str().c_str(), Mesh::DOUBLE, num_dbl, 0, err );
    MSQ_ERRZERO(err);
  }
  return handle;
}

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

//.........这里部分代码省略.........
            norm_dp1 = norm_d + 2.0*alpha*dMp + alpha*alpha*norm_p;
            if (norm_dp1 >= radius) {
                alpha = (sqrt(dMp*dMp+norm_p*(radius-norm_d))-dMp)/norm_p;
                plus_eq_scaled( d, alpha, p, nn );
                break;
            }

            plus_eq_scaled( d, alpha, p, nn );
            plus_eq_scaled( r, alpha, w, nn );
            norm_r = length_squared(r, nn);

            apply_preconditioner( z, r, err);
            MSQ_ERRRTN(err); //prec->apply(z, r, prec, mesh);

            rzm1 = rz;
            rz = inner(r, z, nn);
            beta = rz / rzm1;
            times_eq_minus( p, beta, z, nn );

            dMp = beta*(dMp + alpha*norm_p);
            norm_p = rz + beta*beta*norm_p;
            norm_d = norm_dp1;
        }

#ifdef DO_STEEP_DESC
        if (norm_d <= tr_num_tol) {
            norm_g = length(arrptr(mGrad), nn);
            double ll = 1.0;
            if (norm_g < tr_num_tol)
                break;
            if (norm_g > radius)
                ll = radius / nurm_g;
            for (int i = 0; i < nn; ++i)
                d[i] = ll * mGrad[i];
        }
#endif

        alpha = inner( arrptr(mGrad), d, nn ); // inner(mesh->g, d, nn);

        memset(p, 0, 3*sizeof(double)*nn);
        //matmul(p, mHess, d); //matmul(p, mesh, d);
        mHess.product( p, d );
        beta = 0.5*inner(p, d, nn);
        kappa = alpha + beta;

        /* Put the new point into the locations */
        pd.move_free_vertices_constrained( d, nn, 1.0, 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 not defined at trial point */
            radius *= tr_decr_undef;
            pd.set_to_vertices_memento( mMemento, err );
            MSQ_ERRRTN(err);
            continue;
        }


        if ((fabs(kappa) <= tr_num_tol) && (fabs(objn - obj) <= tr_num_tol)) {
            kappa = 1;
        }
        else {
            kappa = (objn - obj) / kappa;
        }

        if (kappa < eta_1) {
            /* Iterate is unacceptable */
            radius *= tr_decr_def;
            pd.set_to_vertices_memento( mMemento, err );
            MSQ_ERRRTN(err);
            continue;
        }

        /* Iterate is acceptable */
        if (kappa >= eta_2) {
            /* Iterate is a very good step, increase radius */
            radius *= tr_incr;
            if (radius > 1e20) {
                radius = 1e20;
            }
        }

        func.update( pd, obj, mGrad, mHess, err );
        compute_preconditioner( err );
        MSQ_ERRRTN(err);
        pd.recreate_vertices_memento( mMemento, err );
        MSQ_ERRRTN(err);

        // checks stopping criterion
        term.accumulate_patch( pd, err );
        MSQ_ERRRTN(err);
        term.accumulate_inner( pd, objn, arrptr(mGrad), err );
        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);
    if(err.error_code() == err.BARRIER_VIOLATED)
      err.clear();  // barrier violation does not represent an actual error here
    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;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      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); 
      if (err.error_code() == err.BARRIER_VIOLATED)
        err.clear();  // barrier violated does not represent an actual error here
      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 {
        if (!fn_bool) {
  	// Function undefined.  Use the higher decrease rate.
          beta *= beta0;
          MSQ_DBGOUT(3) << "  o  beta = " << beta << " (invalid step)" << std::endl;
        }
        else {
          // Function defined, but not sufficient decrease
          // Use the lower decrease rate.
          beta *= beta1;
          MSQ_DBGOUT(3) << "  o  beta = " << beta << " (insufficient decrease)" << std::endl;
        }
        pd.set_to_vertices_memento(coordsMem, err); MSQ_ERRRTN(err);
      
        // Standard Armijo linesearch rules
 
        MSQ_DBGOUT(3) << "  o  Doing line search" << std::endl;
        while (beta >= tol1) {
          // 6. Search along the direction
          //    (a) trial = x + beta*d
          pd.move_free_vertices_constrained(arrptr(d), nv, beta, err); MSQ_ERRRTN(err);
          //    (b) function evaluation
          fn_bool = objFunc.evaluate(pd, new_value, err); 
          if (err.error_code() == err.BARRIER_VIOLATED)