C++ VectorXd::maxCoeff方法代码示例

int CLogistic::PredictOnevsAll(const VectorXd& point)
{

	// compute prediction for one vs. all
	VectorXd prob = (point * classifier.transpose());

	// pick the most expressive
	return prob.maxCoeff();
}

double CGppe::get_fbest(int N)
{
    VectorXd ftest;
    double fbest;

    ftest = f.segment(f.rows() - N, N - 1);
    fbest = ftest.maxCoeff();
    if (fbest != fbest)
        fbest = f.maxCoeff();
    return fbest;
}

  // Updates the inference result.
  // Hand position is in mm.
  std::string InfEngine::Update(float* raw_feature) {
    using Eigen::Map;
    using Eigen::VectorXf;
    using Eigen::VectorXd;
    using std::string;

    int gesture_index = 0;
    int handpose_index = 0;
    string stage = "Unknown";

    json_spirit::mObject result;

    if (raw_feature != NULL) {
      int motion_feature_len = feature_len_ - n_principal_comps_;

      Map<VectorXf> des(raw_feature + motion_feature_len, descriptor_len_);
      VectorXf res(n_principal_comps_);
      res.noalias() = principal_comp_ * (des - pca_mean_);

      Map<VectorXf> motion_feature(raw_feature, motion_feature_len);
      VectorXf full_feature(feature_len_);
      full_feature << motion_feature, res;
      // Normalize feature.
      full_feature = (full_feature - std_mu_).cwiseQuotient(std_sigma_);

      if (svm_classifier_) {
        VectorXd prob = svm_classifier_->Predict(full_feature);
        VectorXd::Index max_index;
        prob.maxCoeff(&max_index);
        handpose_index = (int) max_index;
      }

      if (hmm_) {
        hmm_->Fwdback(full_feature);
        gesture_index = hmm_->MostLikelyLabelIndex();
        stage = hmm_->MostLikelyStage();
      }

      result["rightX"] = (int) (motion_feature[0] * 1000);
      result["rightY"] = (int) (motion_feature[1] * 1000);
    }

    const string& gesture_label = gesture_labels_[gesture_index];
    gesture_event_detector_.Detect(gesture_label, stage, &result);
    std::string s = write(result, json_spirit::pretty_print | json_spirit::raw_utf8);
    return s;
  }

int main(int argc, char *argv[])
{
  using namespace Eigen;
  using namespace std;

  cout<<"Usage:"<<endl;
  cout<<"[space]  toggle showing surface."<<endl;
  cout<<"'.'/','  push back/pull forward slicing plane."<<endl;
  cout<<endl;

  // Load mesh: (V,T) tet-mesh of convex hull, F contains original surface
  // triangles
  igl::readMESH("../shared/bunny.mesh",V,T,F);


  // Encapsulated call to point_mesh_squared_distance to determine bounds
  {
    VectorXd sqrD;
    VectorXi I;
    MatrixXd C;
    igl::point_mesh_squared_distance(V,V,F,sqrD,I,C);
    max_distance = sqrt(sqrD.maxCoeff());
  }

  // Precompute signed distance AABB tree
  tree.init(V,F);
  // Precompute vertex,edge and face normals
  igl::per_face_normals(V,F,FN);
  igl::per_vertex_normals(
    V,F,igl::PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,FN,VN);
  igl::per_edge_normals(
    V,F,igl::PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,FN,EN,E,EMAP);

  // Plot the generated mesh
  igl::viewer::Viewer viewer;
  update_visualization(viewer);
  viewer.callback_key_down = &key_down;
  viewer.core.show_lines = false;
  viewer.launch();
}

double CGppe::maximum_expected_improvement(const VectorXd & theta_t, const VectorXd& theta_x, const double& sigma,
        const MatrixXd& t, const MatrixXd & x, const VectorXd& idx_global, const VectorXd& ind_t, const VectorXd& ind_x, MatrixXd tstar, int N, double fbest)
{
    VectorXd idx_xstar=Nfirst(N);
    int Kt_ss = 1;
    double  mei;
    MatrixXd Kx_star, Kx_star_star, kstar, Kss, Css;
    MatrixXd Kt_star = covfunc_t->Compute(t, tstar);
	//dsp(GetKinv(),"Kinv");


    Kx_star = GetMatRow(Kx, idx_xstar.transpose()); //maybe need some transpose?

    Kx_star_star = GetMat(Kx, idx_xstar.transpose(), idx_xstar.transpose()); // test to test
    kstar = Kron(Kt_star, Kx_star);

    kstar = GetMatRow(kstar, idx_global);
    Kss = Kt_ss * Kx_star_star;


    mustar = kstar.transpose() * Kinv * GetVec(f, idx_global);
    Css    = Kss - kstar.transpose() * W * llt.solve(Kinv * kstar);
    varstar = Css.diagonal();


    VectorXd sigmastar = sqrt(varstar.array());
    VectorXd z = (fbest - mustar.array()) / sigmastar.array();
    VectorXd pdfval = normpdf(z);
    VectorXd cdfval = normcdf(z);
    VectorXd inter = z.array() * (1 - cdfval.array());
    VectorXd el = sigmastar.cwiseProduct(inter - pdfval);

	el=-1*el;
    mei = el.maxCoeff();
    //dsp(mei,"mei");
    return mei;
}

void  UpdaterMean::costsToWeights(const VectorXd& costs, string weighting_method, double eliteness, VectorXd& weights) const
{
  weights.resize(costs.size());
  if (weighting_method.compare("PI-BB")==0)
  {
    // PI^2 style weighting: continuous, cost exponention
    double h = eliteness; // In PI^2, eliteness parameter is known as "h"
    double range = costs.maxCoeff()-costs.minCoeff();
    if (range==0)
      weights.fill(1);
    else
      weights = (-h*(costs.array()-costs.minCoeff())/range).exp();
  } 
  else if (weighting_method.compare("CMA-ES")==0 || weighting_method.compare("CEM")==0 )
  {
    // CMA-ES and CEM are rank-based, so we must first sort the costs, and the assign a weight to 
    // each rank.
    VectorXd costs_sorted = costs; 
    std::sort(costs_sorted.data(), costs_sorted.data()+costs_sorted.size());
    // In Python this is more elegant because we have argsort.
    // indices = np.argsort(costs)
    // It is possible to do this with fancy lambda functions or std::pair in C++ too, but  I don't
    // mind writing two for loops instead ;-)
    
    weights.fill(0.0);
    int mu = eliteness; // In CMA-ES, eliteness parameter is known as "mu"
    assert(mu<costs.size());
    for (int ii=0; ii<mu; ii++)
    {
      double cur_cost = costs_sorted[ii];
      for (int jj=0; jj<costs.size(); jj++)
      {
        if (costs[jj] == cur_cost)
        {
          if (weighting_method.compare("CEM")==0)
            weights[jj] = 1.0/mu; // CEM
          else
            weights[jj] = log(mu+0.5) - log(ii+1); // CMA-ES
          break;
        }
      }
      
    }
    // For debugging
    //MatrixXd print_mat(3,costs.size());
    //print_mat.row(0) = costs_sorted;
    //print_mat.row(1) = costs;
    //print_mat.row(2) = weights;
    //cout << print_mat << endl;
  }
  else
  {
    cout << __FILE__ << ":" << __LINE__ << ":WARNING: Unknown weighting method '" << weighting_method << "'. Calling with PI-BB weighting." << endl; 
    costsToWeights(costs, "PI-BB", eliteness, weights);
    return;
  }
  
  // Relative standard deviation of total costs
  double mean = weights.mean();
  double std = sqrt((weights.array()-mean).pow(2).mean());
  double rel_std = std/mean;
  if (rel_std<1e-10)
  {
      // Special case: all costs are the same
      // Set same weights for all.
      weights.fill(1);
  }

  // Normalize weights
  weights = weights/weights.sum();

}

bool ClassLayouter::computeNormally()
{
    FuzzyDependAttr::Ptr fuzzyAttr = m_parent->getAttr<FuzzyDependAttr>();
    if (!fuzzyAttr)
        return false;

    Graph G;
    GraphAttributes GA(G,
                       GraphAttributes::nodeGraphics | GraphAttributes::edgeGraphics |
                       GraphAttributes::nodeLabel    | GraphAttributes::nodeTemplate |
                       GraphAttributes::edgeDoubleWeight);

    SparseMatrix& veMat = fuzzyAttr->vtxEdgeMatrix();
    VectorXd	  edgeWeight	= fuzzyAttr->edgeWeightVector();
    if (edgeWeight.size() != veMat.cols())
    {
        m_status |= WARNING_USE_DEFAULT_EDGE_WEIGHT;
        edgeWeight.setOnes(veMat.cols());
    }
    const int nNodes = veMat.rows();
    const int nEdges = veMat.cols();
    if (nNodes <= 0 || nNodes != m_childList.size() || nEdges < 1)
    {
        m_status |= WARNING_NO_EDGE;
        return false;
    }

    bool isUseOrthoLayout = nEdges < 50;

    vector<node> nodeArray;
    vector<edge> edgeArray;
    NodeArray<float> nodeSize(G);
    EdgeArray<double> edgeLength(G);
    for (int i = 0; i < nNodes; ++i)
    {
        nodeArray.push_back(G.newNode());
        float r = m_nodeRadius[i];
        GA.width(nodeArray.back()) = r*2;
        GA.height(nodeArray.back()) = r*2;
        nodeSize[nodeArray.back()] = r * 2;
    }

    float maxEdgeWeight = edgeWeight.maxCoeff();
    float minEdgeWeight = edgeWeight.minCoeff();
    for (int ithEdge = 0; ithEdge < nEdges; ++ithEdge)
    {
        int src, dst;
        GraphUtility::getVtxFromEdge(veMat, ithEdge, src, dst);
        edgeArray.push_back(G.newEdge(nodeArray[src], nodeArray[dst]));
        //GA.doubleWeight(edgeArray.back()) = edgeWeight[ithEdge];
        edgeLength[edgeArray.back()] = 1;//log(edgeWeight[ithEdge] + 1);
    }


    // add dummy vertices and edges in order to merge parallel edge segments attached to the same node
    VectorXi compVec;
    int nComp = 1;
    const float dummyNodeRadius = 0.01;
    if(m_isMergeEdgeEnd && isUseOrthoLayout && GraphUtility::findConnectedComponentsVE(	veMat, compVec, nComp ))
    {
        bool* isCompSet = new bool[nComp];
        for (int i = 0; i < nComp; ++i)
            isCompSet[i] = false;

        // add a dummy node and edge for each connect component
        for (int ithVtx = 0; ithVtx < compVec.size(); ++ithVtx)
        {
            int ithCmp = compVec[ithVtx];
            if (isCompSet[ithCmp] == false)
            {
                // add dummy node and set its radius
                nodeArray.push_back(G.newNode());
                GA.width(nodeArray.back()) = dummyNodeRadius;
                GA.height(nodeArray.back()) = dummyNodeRadius;
                // add dummy edge
                edgeArray.push_back(G.newEdge(nodeArray[ithVtx], nodeArray.back()));
                isCompSet[ithCmp] = true;
            }
        }
        delete[] isCompSet;
    }


    MatrixXd pos;
    pos.resize(nNodes, 2);
    try
    {
        if (isUseOrthoLayout)
        {
            PlanarizationLayout layouter;
            OrthoLayout *ol = new OrthoLayout;
            float sep = max(m_nodeRadius.sum() / nNodes, LayoutSetting::s_baseRadius * 12);
            ol->separation(sep);
            ol->cOverhang(0.1);
            ol->setOptions(1+4);
            layouter.setPlanarLayouter(ol);
            layouter.call(GA);
            for (int v = 0; v < nNodes; ++v)
            {
                double x = GA.x(nodeArray[v]);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    seed = (sal_->getIndxSalient())->at(si);
    if (!isSeedValid())
      continue;
    seeds.push_back(seed);
    ns++;

    // Step 4: fetch and validate neighborhood
    //double last_nn = pcl::getTime ();
    search.radiusSearch(filtered_cloud_->points[seed], radius_, nn_indices, nn_sqr_distances);
    nn_cloud_->points.resize (0);
    for (int i=0; i<nn_indices.size(); i++)
      nn_cloud_->points.push_back(filtered_cloud_->points[nn_indices[i]]);
    //double now_nn = pcl::getTime ();
    //cerr << "Time for nn: " << now_nn-last_nn << endl;
   
    if (do_vis_ && show_nn_ && viewer_)
      AutoSeg::showNN();

    // Step 5: fit the patch
    //double last_pf = pcl::getTime ();
    PatchFit *patchFit;
    patchFit = new PatchFit(nn_cloud_);
    patchFit->setSSmax (50);
    patchFit->fit();
    //double now_pf = pcl::getTime ();
    //cerr << "Time for fitting " << nn_cloud_->points.size() <<
    //        " points: " << now_pf-last_pf << endl;
  

    // Step 6: validate patch
    //double last_val = pcl::getTime ();
    if (do_validation_)
    {
      (*patchFit->p_).computeResidual (nn_cloud_);
      if ((*patchFit->p_).getResidual() > t_residual_)
        continue;

      VectorXd k;
      k = (*patchFit->p_).getK ();
      if (k.minCoeff () < t_min_curv_ || k.maxCoeff ()>t_max_curv_)
        continue;
    }
    num_patches++;
    //double now_val = pcl::getTime ();
    //cerr << "Time for validation: " << now_val-last_val << endl;

    // Visualize patch
    //double last_pp = pcl::getTime ();
    if (do_vis_ && show_patches_ && viewer_)
    {
      PatchPlot *patch_plot;
      (*patchFit->p_).setID (ns);
      (*patchFit->p_).setSS (0.025);
      (*patchFit->p_).infer_params();
      (*patchFit->p_).gs();

      //PatchPlot patch_plot (*patchFit->p_);
      patch_plot = new PatchPlot (*patchFit->p_);
    
      patch_plot->showPatch (viewer_, 0, 1, 0, boost::to_string(ns) + "_patch");
    
      delete (patch_plot);
    }
    //double now_pp = pcl::getTime ();
    //cerr << "Time for ploting patch: " << now_pp-last_pp << endl;

    if (!no_stats_)
    {
      // Save normal vector
      (*patchFit->p_).setCnn (sal_->getNormalNormalAngle (seed));
      if (use_gravity_)
        (*patchFit->p_).setCng (sal_->getNormalGravityAngle (seed));

      // Step 7: save statistics
      AutoSeg::saveStat (*patchFit->p_);
    
      // Step 8: print statistics
      Vector2d sk = (*patchFit->p_).getSK();
      cout << "patch stats" << " "
           << sk (0) << " "
           << sk (1) << " "
           << (*patchFit->p_).getCnn () << " "
           << (*patchFit->p_).getCng () << endl;
    }
    // delete objects
    delete (patchFit);

  } // for each seed

  if(do_vis_ && show_patches_) max_patch_plot_id_ = ns-1;
  else max_patch_plot_id_ = 0;

  // Total timing
  double now = getTime();
  cerr << "Total autoseg for " << num_patches << " patches out of " << ns <<
          " seed points: " << now-last << " sec." << endl;

  // destroy objects
  delete (sal_);
}

//.........这里部分代码省略.........
      beta(j) = (maxStep - rescaledLambda) * betaLS(0) / maxStep;
    }

  } else {

    // further initializations for iterative steps
    VectorXi active;  // active predictors
  	int k = 0;        // number of active predictors
    // previous and current regression coefficients
  	VectorXd previousBeta = VectorXd::Zero(p+s), currentBeta = VectorXd::Zero(p+s);
  	// previous and current penalty parameter
    double previousLambda = 0, currentLambda = 0;
  	// indicates variables to be dropped
    VectorXi drops;
  	// keep track of sign of correlations for the active variables 
    // (double precision is necessary for solve())
    VectorXd signs;
  	// Cholesky L of Gram matrix of active variables
    MatrixXd L;
  	int rank = 0;		// rank of Cholesky L
    // maximum number of variables to be sequenced
  	int maxActive = findMaxActive(n, p, useIntercept);

  	// modified LARS algorithm for lasso solution
  	while(k < maxActive) {

  		// extract current correlations of inactive variables
  		VectorXd corInactiveY(m);
  		for(int j = 0; j < m; j++) {
  			corInactiveY(j) = corY(inactive(j));
  		}
  		// compute absolute values of correlations and find maximum
  		VectorXd absCorInactiveY = corInactiveY.cwiseAbs();
  		double maxCor = absCorInactiveY.maxCoeff();
  		// update current lambda
  		if(k == 0) {	// no active variables
  			previousLambda = maxCor;
  		} else {
  			previousLambda = currentLambda;
  		}
  		currentLambda = maxCor;
  		if(currentLambda <= rescaledLambda) break;

  		if(drops.size() == 0) {
  			// new active variables
  			VectorXi newActive = findNewActive(absCorInactiveY, maxCor - eps);
  			// do calculations for new active variables
  			for(int j = 0; j < newActive.size(); j++) {
  				// update Cholesky L of Gram matrix of active variables
  				// TODO: put this into void function
  				int newJ = inactive(newActive(j));
  				VectorXd xNewJ;
  				double newX;
  				if(useGram) {
  					newX = Gram(newJ, newJ);
  				} else {
  					xNewJ = xs.col(newJ);
  					newX = xNewJ.squaredNorm();
  				}
  				double normNewX = sqrt(newX);
  				if(k == 0) {	// no active variables, L is empty
  					L.resize(1,1);
  					L(0, 0) = normNewX;
  					rank = 1;
  				} else {
  					VectorXd oldX(k);