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

double LogSumExp(const Eigen::VectorXd& x) {
  const double x_max = x.maxCoeff();
  return log((x.array() - x_max).exp().sum()) + x_max;
};

int Lemke(const Eigen::MatrixXd& _M, const Eigen::VectorXd& _q,
          Eigen::VectorXd* _z) {
  int n = _q.size();

  const double zer_tol = 1e-5;
  const double piv_tol = 1e-8;
  int maxiter = 1000;
  int err = 0;

  if (_q.minCoeff() > 0) {
    // LOG(INFO) << "Trivial solution exists.";
    *_z = Eigen::VectorXd::Zero(n);
    return err;
  }

  *_z = Eigen::VectorXd::Zero(2 * n);
  int iter = 0;
  double theta = 0;
  double ratio = 0;
  int leaving  = 0;
  Eigen::VectorXd Be = Eigen::VectorXd::Constant(n, 1);
  Eigen::VectorXd x = _q;
  std::vector<int> bas;
  std::vector<int> nonbas;

  int t = 2 * n + 1;
  int entering = t;

  bas.clear();
  nonbas.clear();

  for (int i = 0; i < n; ++i) {
    bas.push_back(i);
  }

  Eigen::MatrixXd B = -Eigen::MatrixXd::Identity(n, n);

  for (int i = 0; i < bas.size(); ++i) {
    B.col(bas[i]) = _M.col(bas[i]);
  }

  x = -B.partialPivLu().solve(_q);

  Eigen::VectorXd minuxX = -x;
  int lvindex;
  double tval = minuxX.maxCoeff(&lvindex);
  leaving = bas[lvindex];
  bas[lvindex] = t;

  Eigen::VectorXd U = Eigen::VectorXd::Zero(n);
  for (int i = 0; i < n; ++i) {
    if (x[i] < 0)
      U[i] = 1;
  }
  Be = -(B * U);
  x += tval * U;
  x[lvindex] = tval;
  B.col(lvindex) = Be;

  for (iter = 0; iter < maxiter; ++iter) {
    if (leaving == t) {
      break;
    } else if (leaving < n) {
      entering = n + leaving;
      Be = Eigen::VectorXd::Zero(n);
      Be[leaving] = -1;
    } else {
      entering = leaving - n;
      Be = _M.col(entering);
    }

    Eigen::VectorXd d = B.partialPivLu().solve(Be);

    std::vector<int> j;
    for (int i = 0; i < n; ++i) {
      if (d[i] > piv_tol)
        j.push_back(i);
    }
    if (j.empty()) {
      // err = 2;
      break;
    }

    int jSize = static_cast<int>(j.size());
    Eigen::VectorXd minRatio(jSize);
    for (int i = 0; i < jSize; ++i) {
      minRatio[i] = (x[j[i]] + zer_tol) / d[j[i]];
    }
    double theta = minRatio.minCoeff();

    std::vector<int> tmpJ;
    std::vector<double> tmpMinRatio;
    for (int i = 0; i < jSize; ++i) {
      if (x[j[i]] / d[j[i]] <= theta) {
        tmpJ.push_back(j[i]);
        tmpMinRatio.push_back(minRatio[i]);
      }
    }
//    if (tmpJ.empty())
//    {
//.........这里部分代码省略.........

double SumExp(const Eigen::VectorXd& x) {
  const double x_max = x.maxCoeff();
  return (x.array() - x_max).exp().sum() * exp(x_max);
};

int main(int argc, char *argv[])
{

  // Load a mesh
  igl::readOBJ(TUTORIAL_SHARED_PATH "/inspired_mesh.obj", V, F);

  printf("--Initialization--\n");
  V_border = igl::is_border_vertex(V,F);
  igl::adjacency_list(F, VV);
  igl::vertex_triangle_adjacency(V,F,VF,VFi);
  igl::triangle_triangle_adjacency(F,TT,TTi);
  igl::edge_topology(V,F,E,F2E,E2F);

  // Generate "subdivided" mesh for visualization of curl terms
  igl::false_barycentric_subdivision(V, F, Vbs, Fbs);

  // Compute scale for visualizing fields
  global_scale =  .2*igl::avg_edge_length(V, F);

  //Compute scale for visualizing texture
  uv_scale = 0.6/igl::avg_edge_length(V, F);

  // Compute face barycenters
  igl::barycenter(V, F, B);

  // Compute local basis for faces
  igl::local_basis(V,F,B1,B2,B3);

  //Generate random vectors for constraints
  generate_constraints();

  // Interpolate a 2-PolyVector field to be used as the original field
  printf("--Initial solution--\n");
  igl::n_polyvector(V, F, b, bc, two_pv_ori);

  // Post process original field
  // Compute curl_minimizing matchings and curl
  Eigen::MatrixXi match_ab, match_ba;  // matchings across interior edges
  printf("--Matchings and curl--\n");
  double avgCurl = igl::polyvector_field_matchings(two_pv_ori, V, F, true, match_ab, match_ba, curl_ori);
  double maxCurl = curl_ori.maxCoeff();
  printf("original curl -- max: %.5g, avg: %.5g\n", maxCurl,  avgCurl);

  printf("--Singularities--\n");
  // Compute singularities
  igl::polyvector_field_singularities_from_matchings(V, F, V_border, VF, TT, E2F, F2E, match_ab, match_ba, singularities_ori);
  printf("original #singularities: %ld\n", singularities.rows());

  printf("--Cuts--\n");
 // Get mesh cuts based on singularities
  igl::polyvector_field_cut_mesh_with_singularities(V, F, VF, VV, TT, TTi, singularities_ori, cuts_ori);

  printf("--Combing--\n");
// Comb field
  Eigen::MatrixXd combed;
  igl::polyvector_field_comb_from_matchings_and_cuts(V, F, TT, E2F, F2E, two_pv_ori, match_ab, match_ba, cuts_ori, combed);

  printf("--Cut mesh--\n");
  // Reconstruct integrable vector fields from combed field
  igl::cut_mesh(V, F, VF, VFi, TT, TTi, V_border, cuts_ori, Vcut_ori, Fcut_ori);

  printf("--Poisson--\n");
  double avgPoisson = igl::polyvector_field_poisson_reconstruction(Vcut_ori, Fcut_ori, combed, scalars_ori, two_pv_poisson_ori, poisson_error_ori);
  double maxPoisson = poisson_error_ori.maxCoeff();
  printf("poisson error -- max: %.5g, avg: %.5g\n", maxPoisson, avgPoisson);


  // Set the curl-free 2-PolyVector to equal the original field
  two_pv = two_pv_ori;
  singularities = singularities_ori;
  curl = curl_ori;
  cuts = cuts_ori;
  two_pv_poisson = two_pv_poisson_ori;
  poisson_error = poisson_error_ori;
  Vcut = Vcut_ori;
  Fcut = Fcut_ori;
  scalars = scalars_ori;

  printf("--Integrable - Precomputation--\n");
  // Precompute stuff for solver
  igl::integrable_polyvector_fields_precompute(V, F, b, bc, blevel, two_pv_ori, ipfdata);

  cerr<<"Done. Press keys 1-0 for various visualizations, 'A' to improve integrability." <<endl;

  igl::viewer::Viewer viewer;
  viewer.callback_key_down = &key_down;
  viewer.core.show_lines = false;
  key_down(viewer,'2',0);

  // Replace the standard texture with an integer shift invariant texture
  line_texture(texture_R, texture_G, texture_B);

  viewer.launch();

  return 0;
}

bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
{

  if (key == '1')
  {
    display_mode = 1;
    update_display(viewer);
  }

  if (key == '2')
  {
    display_mode = 2;
    update_display(viewer);
  }

  if (key == '3')
  {
    display_mode = 3;
    update_display(viewer);
  }

  if (key == '4')
  {
    display_mode = 4;
    update_display(viewer);
  }

  if (key == '5')
  {
    display_mode = 5;
    update_display(viewer);
  }

  if (key == '6')
  {
    display_mode = 6;
    update_display(viewer);
  }

  if (key == '7')
  {
    display_mode = 7;
    update_display(viewer);
  }

  if (key == '8')
  {
    display_mode = 8;
    update_display(viewer);
  }

  if (key == '9')
  {
    display_mode = 9;
    update_display(viewer);
  }

  if (key == '0')
  {
    display_mode = 0;
    update_display(viewer);
  }

  if (key == 'A')
  {
    //do a batch of iterations
    printf("--Improving Curl--\n");
    for (int bi = 0; bi<5; ++bi)
    {
      printf("\n\n **** Batch %d ****\n", iter);
      igl::integrable_polyvector_fields_solve(ipfdata, params, two_pv, iter ==0);
      iter++;
      params.wSmooth *= params.redFactor_wsmooth;
    }
    // Post process current field
    // Compute curl_minimizing matchings and curl
    printf("--Matchings and curl--\n");
    Eigen::MatrixXi match_ab, match_ba;  // matchings across interior edges
    double avgCurl = igl::polyvector_field_matchings(two_pv, V, F, true, match_ab, match_ba, curl);
    double maxCurl = curl.maxCoeff();
    printf("curl -- max: %.5g, avg: %.5g\n", maxCurl,  avgCurl);
    // Compute singularities
    printf("--Singularities--\n");
    igl::polyvector_field_singularities_from_matchings(V, F, match_ab, match_ba, singularities);
    printf("#singularities: %ld\n", singularities.rows());
    // Get mesh cuts based on singularities
    printf("--Cuts--\n");
    igl::polyvector_field_cut_mesh_with_singularities(V, F, singularities, cuts);
    // Comb field
    printf("--Combing--\n");
    Eigen::MatrixXd combed;
    igl::polyvector_field_comb_from_matchings_and_cuts(V, F, two_pv, match_ab, match_ba, cuts, combed);
    // Reconstruct integrable vector fields from combed field
    printf("--Cut mesh--\n");
    igl::cut_mesh(V, F, cuts, Vcut, Fcut);
    printf("--Poisson--\n");
    double avgPoisson = igl::polyvector_field_poisson_reconstruction(Vcut, Fcut, combed, scalars, two_pv_poisson, poisson_error);
    double maxPoisson = poisson_error.maxCoeff();
    printf("poisson error -- max: %.5g, avg: %.5g\n", maxPoisson, avgPoisson);

//.........这里部分代码省略.........

void propa_2d()
{
    trace.beginBlock("2d propagation");

    typedef DiscreteExteriorCalculus<2, 2, EigenLinearAlgebraBackend> Calculus;
		typedef DiscreteExteriorCalculusFactory<EigenLinearAlgebraBackend> CalculusFactory;

    {
        trace.beginBlock("solving time dependent equation");

        const Calculus::Scalar cc = 8; // px/s
        trace.info() << "cc = " << cc << endl;

        const Z2i::Domain domain(Z2i::Point(0,0), Z2i::Point(29,29));
        const Calculus calculus = CalculusFactory::createFromDigitalSet(generateDiskSet(domain), false);

        //! [time_laplace]
        const Calculus::DualIdentity0 laplace = calculus.laplace<DUAL>() + 1e-8 * calculus.identity<0, DUAL>();
        //! [time_laplace]
        trace.info() << "laplace = " << laplace << endl;

        trace.beginBlock("finding eigen pairs");
        //! [time_eigen]
        typedef Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> EigenSolverMatrix;
        const EigenSolverMatrix eigen_solver(laplace.myContainer);

        const Eigen::VectorXd eigen_values = eigen_solver.eigenvalues();
        const Eigen::MatrixXd eigen_vectors = eigen_solver.eigenvectors();
        //! [time_eigen]
        trace.info() << "eigen_values_range = " << eigen_values.minCoeff() << "/" << eigen_values.maxCoeff() << endl;
        trace.endBlock();

        //! [time_omega]
        const Eigen::VectorXd angular_frequencies = cc * eigen_values.array().sqrt();
        //! [time_omega]

        Eigen::VectorXcd initial_wave = Eigen::VectorXcd::Zero(calculus.kFormLength(0, DUAL));
        //set_initial_phase_null(calculus, initial_wave);
        set_initial_phase_dir_yy(calculus, initial_wave);

        {
            Board2D board;
            board << domain;
            board << CustomStyle("KForm", new KFormStyle2D(-1, 1));
            board << Calculus::DualForm0(calculus, initial_wave.real());
            board.saveSVG("propagation_time_wave_initial_coarse.svg");
        }

        //! [time_init_proj]
        Eigen::VectorXcd initial_projections = eigen_vectors.transpose() * initial_wave;
        //! [time_init_proj]

        // low pass
        //! [time_low_pass]
        const Calculus::Scalar lambda_cutoff = 4.5;
        const Calculus::Scalar angular_frequency_cutoff = 2*M_PI * cc / lambda_cutoff;
        int cutted = 0;
        for (int kk=0; kk<initial_projections.rows(); kk++)
        {
            const Calculus::Scalar angular_frequency = angular_frequencies(kk);
            if (angular_frequency < angular_frequency_cutoff) continue;
            initial_projections(kk) = 0;
            cutted ++;
        }
        //! [time_low_pass]
        trace.info() << "cutted = " << cutted << "/" << initial_projections.rows() << endl;

        {
            const Eigen::VectorXcd wave = eigen_vectors * initial_projections;
            Board2D board;
            board << domain;
            board << CustomStyle("KForm", new KFormStyle2D(-1, 1));
            board << Calculus::DualForm0(calculus, wave.real());
            board.saveSVG("propagation_time_wave_initial_smoothed.svg");
        }

        const int kk_max = 60;
        trace.progressBar(0,kk_max);
        for (int kk=0; kk<kk_max; kk++)
        {
            //! [time_solve_time]
            const Calculus::Scalar time = kk/20.;
            const Eigen::VectorXcd current_projections = (angular_frequencies * std::complex<double>(0,time)).array().exp() * initial_projections.array();
            const Eigen::VectorXcd current_wave = eigen_vectors * current_projections;
            //! [time_solve_time]

            std::stringstream ss;
            ss << "propagation_time_wave_solution_" << std::setfill('0') << std::setw(3) << kk << ".svg";

            Board2D board;
            board << domain;
            board << CustomStyle("KForm", new KFormStyle2D(-1, 1));
            board << Calculus::DualForm0(calculus, current_wave.real());
            board.saveSVG(ss.str().c_str());

            trace.progressBar(kk+1,kk_max);
        }
        trace.info() << endl;

        trace.endBlock();
//.........这里部分代码省略.........

Eigen::RowVectorXd scaledValues(Eigen::VectorXd const& x)
{
    return x.unaryExpr(scaledValue(x.minCoeff(),x.maxCoeff())).transpose();
}