C++ VectorXd::sum方法代碼示例

void normalize(Eigen::VectorXd& v)
{
	if (v.sum() == 0.0)
	{
		std::cout << "ERRO no portfolio:\n";
		std::cout << v;
		std::cout << "ERRO no portfolio.";
		system("pause");
	}
	v = (1.0/ v.sum())*v;
}

gaussian_process::gaussian_process( const Eigen::MatrixXd& domains
                                , const Eigen::VectorXd& targets
                                , const gaussian_process::covariance& covariance
                                , double self_covariance )
    : domains_( domains )
    , targets_( targets )
    , covariance_( covariance )
    , self_covariance_( self_covariance )
    , offset_( targets.sum() / targets.rows() )
    , K_( domains.rows(), domains.rows() )
{
    if( domains.rows() != targets.rows() ) { COMMA_THROW( comma::exception, "expected " << domains.rows() << " row(s) in targets, got " << targets.rows() << " row(s)" ); }
    targets_.array() -= offset_; // normalise
    //use m_K as Kxx + variance*I, then invert it
    //fill Kxx with values from covariance function
    //for elements r,c in upper triangle
    for( std::size_t r = 0; r < std::size_t( domains.rows() ); ++r )
    {
        K_( r, r ) = self_covariance_;
        const Eigen::VectorXd& row = domains.row( r );
        for( std::size_t c = r + 1; c < std::size_t( domains.rows() ); ++c )
        {
            K_( c, r ) = K_( r, c ) = covariance_( row, domains.row( c ) );
        }
    }
    L_.compute( K_ ); // invert Kxx + variance * I to become (by definition) B
    alpha_ = L_.solve( targets_ );
}

void PCATestCase::decorrelation()
{
  OpenANN::RandomNumberGenerator rng;
  const int N = 100;
  const int D = 2;
  Eigen::MatrixXd X(N, D);
  rng.fillNormalDistribution(X);

  // Linear transformation (correlation)
  Eigen::MatrixXd A(D, D);
  A << 1.0, 0.5, 0.5, 2.0;
  Eigen::MatrixXd Xt = X * A.transpose();
  // Decorrelation (without dimensionality reduction)
  OpenANN::PCA pca(D);
  pca.fit(Xt);
  Eigen::MatrixXd Y = pca.transform(Xt);

  // Covariance matrix should be identity matrix
  Eigen::MatrixXd cov = Y.transpose() * Y;
  ASSERT_EQUALS_DELTA(cov(0, 0), (double) N, 1e-5);
  ASSERT_EQUALS_DELTA(cov(1, 1), (double) N, 1e-5);
  ASSERT_EQUALS_DELTA(cov(0, 1), 0.0, 1e-5);
  ASSERT_EQUALS_DELTA(cov(1, 0), 0.0, 1e-5);

  Eigen::VectorXd evr = pca.explainedVarianceRatio();
  double evrSum = evr.sum();
  ASSERT_EQUALS_DELTA(evrSum, 1.0, 1e-5);
}

 Eigen::VectorXd grad(const Eigen::VectorXd &x1, const Eigen::VectorXd &x2) const {
   Eigen::VectorXd grad(_input_dim + 1);
   Eigen::VectorXd z = (x1 - x2).cwiseQuotient(_ell).array().square();
   double k = _sf2 * exp(-0.5 * z.sum());
   grad.head(_input_dim) = z * k;
   grad(_input_dim) = 2.0 * k;
   return grad;
 }

    Solver(const int N, const int m, double* alpha_in, double* beta_in, double* t_in, double* yerr) {


        Eigen::Map<Eigen::VectorXd> alpha(alpha_in, m);
        Eigen::Map<Eigen::VectorXd> beta(beta_in, m);
        Eigen::Map<Eigen::VectorXd> t(t_in, N);
        Eigen::Map<Eigen::VectorXd> d(yerr, N);
        for (int i=0; i<N; i++) {
            d(i) = alpha.sum() + yerr[i]*yerr[i];
        }
        this->N = N;
        this->m = m;
        grp_ = new GRP(N, m, alpha, beta, t, d);

    };

	void calcMeanAndCovarWeighedVectorized(const Eigen::MatrixXf &input, const Eigen::VectorXd &inputWeights, Eigen::MatrixXf &out_covMat, Eigen::VectorXf &out_mean,Eigen::MatrixXf &temp)
	{
		out_mean=input.col(0); //to resize
		out_mean.setZero();
		double wSumInv=1.0/inputWeights.sum();
		for (int k=0;k<inputWeights.size();k++){
			double w=inputWeights[k];
			out_mean+=input.col(k)*(float)(w*wSumInv);
		}
		out_mean = input.rowwise().mean();
		temp = (input.colwise() - out_mean);
		for (int k=0;k<inputWeights.size();k++){
			temp.col(k) *= (float)(sqrt(inputWeights(k)*wSumInv));	//using square roots, as we only want the normalized weights to be included once for each result element in the multiplication below
		}
		out_covMat = temp*temp.transpose();
	}

 SimuEuro(const option & o, long path, const std::vector<double> & RN){
     opt=o;
     N= path;
     asset_price.resize(N);
     asset_price.setZero();
     option_value= asset_price;
     
     for (long i=0; i< N; i++) {
         asset_price(i)=opt.S* exp((opt.r- opt.q)*opt.T-.5*opt.sigma*opt.sigma*opt.T+ opt.sigma* sqrt(opt.T)* RN[i]);
         
         if(opt.Call) option_value(i)= fmax(asset_price(i)- opt.K,0.0);
         else option_value(i)= fmax(-asset_price(i)+opt.K, 0.0);
     }
     
     mean= option_value.sum()/ option_value.size() * exp(-opt.T*opt.r);
     stdiv= option_value.squaredNorm()/ option_value.size()* exp(-opt.r*opt.T *2);
     stdiv= stdiv- pow(mean,2.0);
     stdiv= sqrt(stdiv/ N);
 };

int cMathUtil::SampleDiscreteProb(const Eigen::VectorXd& probs)
{
	assert(std::abs(probs.sum() - 1) < 0.00001);
	double rand = RandDouble();

	int rand_idx = gInvalidIdx;
	int num_probs = static_cast<int>(probs.size());
	for (int i = 0; i < num_probs; ++i)
	{
		double curr_prob = probs[i];
		rand -= curr_prob;

		if (rand <= 0)
		{
			rand_idx = i;
			break;
		}
	}
	return rand_idx;
}

      IEFSolver aniso_solver(symm);
      aniso_solver.buildAnisotropicMatrix(cavity, gf_i, gf_o, op);

      IEFSolver iso_solver(symm);
      iso_solver.buildIsotropicMatrix(cavity, gf_i, gf_o, op);

      int size = cavity.size();
      Eigen::VectorXd fake_mep = computeMEP(cavity.elements(), charge);

      THEN("the total apparent surface charge is") {
        Eigen::VectorXd aniso_fake_asc = aniso_solver.computeCharge(fake_mep);
        Eigen::VectorXd iso_fake_asc = iso_solver.computeCharge(fake_mep);

        int nr_irrep = cavity.pointGroup().nrIrrep();
        double totalAnisoASC = aniso_fake_asc.sum() * nr_irrep;
        double totalIsoASC = iso_fake_asc.sum() * nr_irrep;

        CAPTURE(totalASC);
        CAPTURE(totalAnisoASC);
        CAPTURE(totalASC - totalAnisoASC);
        REQUIRE(totalIsoASC == Approx(totalAnisoASC));
        REQUIRE(totalASC == Approx(totalAnisoASC).epsilon(1.0e-03));
      }
    }

    /*! \class IEFSolver
     *  \test \b pointChargeGePolC2 tests IEFSolver using a point charge with a GePol
     * cavity in C2 symmetry
     *  We are forcing the usage of the buildAnisotropicMatrix method.
     *  The results are compared with Gauss' theorem and the results from the

SEXP cleanZ(SEXP X,SEXP Z,SEXP bic_vide_vect,SEXP methode_BIC,SEXP plot,SEXP bla,SEXP Rstar)
{
 
  BEGIN_RCPP
  //déclaration des varibles
  const Map<MatrixXd> matZ(as<Map<MatrixXd> >(Z));//Z
  const Map<MatrixXd> matX(as<Map<MatrixXd> >(X));//X
  const Map<VectorXd> Bic_vide_vect(as<Map<VectorXd> >(bic_vide_vect));//bic_vide_vect
  bool star = Rcpp::as<bool>(Rstar);     // BICstar

  Rcpp::NumericVector met_BIC(methode_BIC),Plot(plot),Bla(bla);
  typedef Rcpp::NumericVector::iterator vec_iterator;
  vec_iterator imet_BIC=met_BIC.begin(),iplot=Plot.begin(),ibla=Bla.begin();
  int p=matZ.cols();//nombre de colonne de la matrice Z
  Eigen::MatrixXd Zopt;//meilleur Z obtenu
  double Bicbest;//BIC associe au meilleur modèle
  Eigen::VectorXd bicvect;//vecteur BIC des matrices Z retenues
  Eigen::MatrixXd newZ;//matrice Z modifie a chaque etapes
  Eigen::MatrixXd list_cand;//matrice qui contient les coordonnées des candidats (liste candidats)
  int compte;//permet de cree la matrice liste en désignant le numéro du candidat (liste candidats)
  int nbcand;//nombre de candidats
  int numcand;//numero du candidat
  int i_loc;//premiere coordonnee du candidat (modification Z)
  int j_loc;//duexieme coordonnee du candidat (modification Z)
  Eigen::MatrixXd Zcand;//matrice Z du candidat (modification Z)
  Eigen::VectorXd BIC_cand;//vecteur qui contient les BIC de chaque colonne des la matrice Zcand (calcul du BIC)
  double Sum_BIC_cand;//somme des BIC de BIC_cand (calcul BIC)
  Eigen::VectorXd stock_bool;//vecteur qui permet de savoir si un candidat peut etre choisi (stock)
  Eigen::VectorXd stock_BIC;//vecteur qui contient le BIC de tout les candidats (stock)
  double sumbic;//BIC de chaque etapes
  double BIC_min;
  bool minimum;
  int i;
  int iret;
  nbcand=matZ.sum();
  Eigen::VectorXd bic_etape(nbcand);//vecteur qui stock le BIC de chaque etapes
  Eigen::VectorXd complexite_etape(nbcand);//vecteur qui stock la complexite de Z a chaque etapes
  Eigen::VectorXd BIC_step;//vecteur qui stock le BIC de chaque etapes
  Eigen::VectorXd complexite_step ;//vecteur qui stock la compléxité de Z a chaque etapes
//  Eigen::ArrayXXd SumCol(1,p);//SumCol est un vecteur qui contient la somme de chaque colonne de Zcand (modification Z)

    
  //initialisation
  bicvect=BicZ_cpp2(matX,matZ,Bic_vide_vect,imet_BIC[0]);
//somme a la main
  sumbic=bicvect.sum();
  if (star)
  {
      sumbic=sumbic-ProbaZ_cpp(matZ);
  }
  if (ibla[0]>0)
  {
    Rcout<<sumbic<<"\n";
  }
  Bicbest =sumbic;
  Zopt=matZ;
  newZ=matZ;
  compte=nbcand+1;
  int step =0;
  
  while(step<nbcand)  
  {
    list_cand.resize(compte-1,2);//le nombre de candidats est la complexite du modele
    compte=0;//initialisation du vecteur liste designe le numero du candidat
  //liste candidats (couples [i,j])
//    list_cand.resize(nbcand,2);//le nombre de candidats est la complexite du modele
    
    for(int m=0;m<p;m++)//parcours ligne
    {
      for(int n=0;n<p;n++)//parcours colonne
      {
        if(newZ(m,n)==1)//on cherche les candidats
        {
          list_cand(compte,0)=m;//stock la ligne du candidat
          list_cand(compte,1)=n;//stock la colonne du candidat
          compte=compte+1;//passage au candidat suivant
        }
      }
    }
    //pour chaque candidat
    for (numcand=0;numcand<compte;numcand++)
    {
    //modification (calcul du Z)    (Zcand avec methode (rejet ou relax)   Z,i,j,methode="relax", p2max=inf,rmax=inf)
      i_loc=list_cand(numcand,0);
      j_loc=list_cand(numcand,1);
      Zcand=newZ;
      Zcand(i_loc,j_loc)=1-Zcand(i_loc,j_loc);//modification de Z par rapport au candidat selectionne
      
    //calcul du bic (du nouveau Z genere)      (bicZ)
      BIC_cand=BicZ_cpp(matX,Zcand,Bic_vide_vect,bicvect,imet_BIC[0],newZ);//calcul du vecteur BIC du candidat
      Sum_BIC_cand=BIC_cand.sum();
      if (star)
      {
        Sum_BIC_cand=Sum_BIC_cand-ProbaZ_cpp(Zcand);
      }
    //stockage du BIC
      stock_BIC.resize(compte);
      stock_BIC(numcand)=Sum_BIC_cand;
    }
    //choix du candidat retenu (tirage au sort ou meilleur selon le cas)
//.........這裏部分代碼省略.........

SEXP rechercheZ_relax(SEXP X,SEXP Z,SEXP bic_vide_vect,SEXP methode_tirage,SEXP methode_BIC,SEXP Rmax,SEXP p2max,SEXP Maxiter,SEXP plot,SEXP best,SEXP better,SEXP random,SEXP bla,SEXP nb_opt_max,SEXP Rexact,SEXP Rstar)
{
 
  BEGIN_RCPP
  //déclaration des varibles
  const Map<MatrixXd> matZ(as<Map<MatrixXd> >(Z));//Z
  const Map<MatrixXd> matX(as<Map<MatrixXd> >(X));//X
  const Map<VectorXd> Bic_vide_vect(as<Map<VectorXd> >(bic_vide_vect));//bic_vide_vect
  bool Exact = Rcpp::as<bool>(Rexact);     // length vector 
  bool star = Rcpp::as<bool>(Rstar);     // BICstar

  Rcpp::NumericVector met_tirage(methode_tirage),met_BIC(methode_BIC),rmax(Rmax),P2max(p2max),maxiter(Maxiter),Plot(plot),Best(best),Better(better),Random(random),Bla(bla),Nb_opt_max(nb_opt_max);
  typedef Rcpp::NumericVector::iterator vec_iterator;
  vec_iterator imet_tirage = met_tirage.begin(),imet_BIC=met_BIC.begin(),irmax=rmax.begin(),ip2max=P2max.begin(),imaxiter=maxiter.begin(),iplot=Plot.begin(),ibest=Best.begin(),ibetter=Better.begin(),irandom=Random.begin(),ibla=Bla.begin(),inb_opt_max=Nb_opt_max.begin();
  int p=matZ.cols();//nombre de colonne de la matrice Z
  Eigen::MatrixXd Zopt;//meilleur Z obtenu
  double Bicbest;//BIC associé au meilleur modèle
  Eigen::VectorXd bicvect;//vecteur BIC des matrices Z retenues
  Eigen::MatrixXd newZ;//matrice Z modifié à chaque étapes
  Eigen::MatrixXd list_cand;//matrice qui contient les coordonnées des candidats
  int nbcand=0;//nombre de candidats
  int nb_opt;//nombre de fois ou on a retrouve bicbest
  int step_opt=0;//étape où l'on découvre BIC_opt
  int sumZ=0;
  int k;//nombre du tirage aleatoire (liste candidats)
  int compte;//permet de cree la matrice liste en désignant le numéro du candidat (liste candidats)
  int i;//coordonnée (x ou y) des candidats (liste candidats)
  int rand1;//nombre aleatoire pour la 3ème méthode de tirage des candidats (liste candidats)
  int rand2;//nombre aleatoire pour la 3ème méthode de tirage des candidats (liste candidats)
  int numcand;//numero du candidat
  int i_loc;//premiere coordonnée du candidat (modification Z)
  int j_loc;//duexième coordonnée du candidat (modification Z)
  int realisable;//booléen qui permet de savoir si on peut effectuer un changement dans Z (modification Z)
  Eigen::MatrixXd Zcand;//matrice Z du candidat (modification Z)
  Eigen::ArrayXXd SumCol(1,Zcand.cols());//SumCol est un vecteur qui contient la somme de chaque colonne de Zcand (modification Z)
  int val;//permet de mettre une colonne et une ligne a 0 (modification Z)
  Eigen::VectorXd BIC_cand;//vecteur qui contient les BIC de chaque colonne des la matrice Zcand (calcul du BIC)
  double Sum_BIC_cand;//somme des BIC de BIC_cand (calcul BIC)
  Eigen::VectorXd stock_bool;//vecteur qui permet de savoir si un candidat peut etre choisi (stock)
  Eigen::VectorXd stock_BIC;//vecteur qui contient le BIC de tout les candidats (stock)
  double sumbic;//BIC de chaque etapes
  
  Eigen::VectorXd bic_etape (imaxiter[0]);//vecteur qui stock le BIC de chaque etapes
  Eigen::VectorXd complexite_etape (imaxiter[0]);//vecteur qui stock la compléxité de Z a chaque etapes
  Eigen::VectorXd etape (imaxiter[0]);//vecteur qui stock le type de changement de chaque etapes
  bool station;//permet de savoir si on est stationnaire ou non
    
  //initialisation
  Zopt = MatrixXd::Zero(p,p);
//somme a la main
  Bicbest = Bic_vide_vect.sum();
  bicvect=BicZ_cpp(matX,matZ,Bic_vide_vect,Bic_vide_vect,imet_BIC[0],Zopt);
//somme a la main
  sumbic=bicvect.sum();
  if(star){
      sumbic=sumbic-ProbaZ_cpp(matZ);
   }
  if (ibla[0]>0)
  {
    Rcout<<sumbic<<"\n";
  }
  if (sumbic<Bicbest)
  {
    Bicbest =sumbic;
    Zopt=matZ;
  }
  newZ=matZ;
  
  int step =0;
  
  if (imet_tirage[0]==0)//methode de changement de la ligne et de la colonne
  {
    nbcand=2*p-2;
    list_cand.resize(nbcand,2);//2p-2 candidats
  }
  else if(imet_tirage[0]==-1)
  {
    nbcand=p-1;
    list_cand.resize(nbcand,2);//p-1 candidats (colonne uniquement)
  }
  else if(imet_tirage[0]==-2)
  {
    nbcand=p*(p-1);
    list_cand.resize(nbcand,2);//p-1 candidats (colonne uniquement)
  }
  else if(imet_tirage[0]>0)//methode de tirage aléatoire
  {
    nbcand=imet_tirage[0];
    list_cand.resize(nbcand,2);//le nombre de candidats est determine pas l'utilisateur
  }
  else if(imet_tirage[0]!=-3) //si methode a une valeur aberante
  {
    throw std::range_error("methode de tirage incorrecte");
  }
  
  if (irandom[0]==0)
  {
    ibetter[0]=1;
    ibest[0]=1;
  }
//.........這裏部分代碼省略.........

SEXP cleancolZ(SEXP X,SEXP Z,SEXP bic_vide_vect,SEXP methode_BIC,SEXP plot,SEXP bla,SEXP Rstar)
{
  BEGIN_RCPP
  //declaration des variables
  const Map<MatrixXd> matZ(as<Map<MatrixXd> >(Z));//Z
  const Map<MatrixXd> matX(as<Map<MatrixXd> >(X));//X
  const Map<VectorXd> Bic_vide_vect(as<Map<VectorXd> >(bic_vide_vect));//bic_vide_vect
  bool star = Rcpp::as<bool>(Rstar);     // BICstar

  Rcpp::NumericVector met_BIC(methode_BIC),Plot(plot),Bla(bla);
  typedef Rcpp::NumericVector::iterator vec_iterator;
  vec_iterator imet_BIC=met_BIC.begin(),iplot=Plot.begin(),ibla=Bla.begin();
  int p=matZ.cols();//nombre de colonne de la matrice Z
  Eigen::MatrixXd Zopt;//meilleur Z obtenu
  double Bicbest;//BIC associe au meilleur modele
  Eigen::VectorXd bicvect;//vecteur BIC des matrices Z retenues
  Eigen::MatrixXd newZ;//matrice Z modifie a chaque etapes
  Eigen::VectorXd list_cand;//vecteur qui contient les indices des colonnes candidates (liste candidats)
  int compte=0;//permet de cree la matrice liste en designant le numero du candidat (liste candidats)
  int nbcand;//nombre de candidats
  int numcand;//numero du candidat
  int i_loc;//premiere coordonnee du candidat (modification Z)
  Eigen::MatrixXd Zcand;//matrice Z du candidat (modification Z)
  Eigen::ArrayXXd SumCol(1,p);//SumCol est un vecteur qui contient la somme de chaque colonne de Zcand (modification Z)
  Eigen::VectorXd BIC_cand;//vecteur qui contient les BIC de chaque colonne des la matrice Zcand (calcul du BIC)
  double Sum_BIC_cand;//somme des BIC de BIC_cand (calcul BIC)
  Eigen::VectorXd stock_BIC;//vecteur qui contient le BIC de tout les candidats (stock)
  double sumbic;//BIC de chaque etapes
  double BIC_min;
  bool minimum;
  int i;
  int iret=1;
  Eigen::VectorXd bic_etape;//vecteur qui stock le BIC de chaque etapes
  Eigen::VectorXd complexite_etape;//vecteur qui stock la complexite de Z a chaque etapes
  Eigen::VectorXd BIC_step;//vecteur qui stock le BIC de chaque etapes
  Eigen::VectorXd complexite_step ;//vecteur qui stock la complexite de Z a chaque etapes
  //initialisation
  bicvect=BicZ_cpp2(matX,matZ,Bic_vide_vect,imet_BIC[0]);
//somme a la main
  sumbic=bicvect.sum();
  if(star)
  {
      sumbic=sumbic-ProbaZ_cpp(matZ);
   }
  if (ibla[0]>0)
  {
    Rcout<<sumbic<<"\n";
  }
  Bicbest =sumbic;
  Zopt=matZ;
  newZ=matZ;
  
  int step =0;
  SumCol=matZ.colwise().sum();//nombre d'elements dans chaque colonne
  nbcand=(SumCol>0).count();//nombre de colonnes candidates (quand il y'a au moins 1 element dans une colonne)
  list_cand.resize(nbcand);//le nombre de candidats est la complexite du modele
  bic_etape.resize(nbcand);
  complexite_etape.resize(nbcand);
  while(step<nbcand)  
  {
    compte=0;//initialisation du vecteur liste designe le numero du candidat
    
  //liste candidats (couples [i,j])
    for(int n=0;n<p;n++)//parcours les colonnes
    {
      if(newZ.col(n).sum()>0)//on cherche les candidats
      {
        list_cand(compte)=n;//stock la colonne du candidat
        compte=compte+1;//passage au candidat suivant
      }
    }
    
    stock_BIC.resize(compte);
    //pour chaque candidat
    for (numcand=0;numcand<compte;numcand++)
    {
    //modification (calcul du Z)    (Zcand avec methode (rejet ou relax)   Z,i,j,methode="relax", p2max=inf,rmax=inf)
      
      Zcand=newZ;
      for(i_loc=0;i_loc<p;i_loc++)//parcours des lignes
      {
        Zcand(i_loc,list_cand(numcand))=0;//modification de Z par rapport a la colonne candidate
      }
      
    //calcul du bic (du nouveau Z genere)
      BIC_cand=BicZ_cpp(matX,Zcand,Bic_vide_vect,bicvect,imet_BIC[0],newZ);//calcul du vecteur BIC du candidat
      Sum_BIC_cand=BIC_cand.sum();
      if (star)
      {
         Sum_BIC_cand=Sum_BIC_cand-ProbaZ_cpp(Zcand);
      }
    //stockage du BIC
      
      stock_BIC(numcand)=Sum_BIC_cand;
    }
    //on regarde la meilleur candidat
    BIC_min=Bicbest;//initialisation BIC_min
    minimum=false;
    for(i=0;i<compte;i++)
    {
//.........這裏部分代碼省略.........

int main(int argc, char *argv[])
{
  using namespace std;
  // Load a mesh in OFF format
  igl::readOBJ(TUTORIAL_SHARED_PATH "/camel_b.obj", V, F);

  Eigen::MatrixXd bnd_uv, uv_init;

  Eigen::VectorXd M;
  igl::doublearea(V, F, M);
  std::vector<std::vector<int>> all_bnds;
  igl::boundary_loop(F, all_bnds);

  // Heuristic primary boundary choice: longest
  auto primary_bnd = std::max_element(all_bnds.begin(), all_bnds.end(), [](const std::vector<int> &a, const std::vector<int> &b) { return a.size()<b.size(); });

  Eigen::VectorXi bnd = Eigen::Map<Eigen::VectorXi>(primary_bnd->data(), primary_bnd->size());

  igl::map_vertices_to_circle(V, bnd, bnd_uv);
  bnd_uv *= sqrt(M.sum() / (2 * igl::PI));
  if (all_bnds.size() == 1)
  {
    if (bnd.rows() == V.rows()) // case: all vertex on boundary
    {
      uv_init.resize(V.rows(), 2);
      for (int i = 0; i < bnd.rows(); i++)
        uv_init.row(bnd(i)) = bnd_uv.row(i);
    }
    else
    {
      igl::harmonic(V, F, bnd, bnd_uv, 1, uv_init);
      if (igl::flipped_triangles(uv_init, F).size() != 0)
        igl::harmonic(F, bnd, bnd_uv, 1, uv_init); // fallback uniform laplacian
    }
  }
  else
  {
    // if there is a hole, fill it and erase additional vertices.
    all_bnds.erase(primary_bnd);
    Eigen::MatrixXi F_filled;
    igl::topological_hole_fill(F, bnd, all_bnds, F_filled);
    igl::harmonic(F_filled, bnd, bnd_uv ,1, uv_init);
    uv_init = uv_init.topRows(V.rows());
  }

  Eigen::VectorXi b; Eigen::MatrixXd bc;
  igl::scaf_precompute(V, F, uv_init, scaf_data, igl::MappingEnergyType::SYMMETRIC_DIRICHLET, b, bc, 0);

  // Plot the mesh
  igl::opengl::glfw::Viewer viewer;
  viewer.data().set_mesh(V, F);
  const auto& V_uv = uv_scale * scaf_data.w_uv.topRows(V.rows());
  viewer.data().set_uv(V_uv);
  viewer.callback_key_down = &key_down;

  // Enable wireframe
  viewer.data().show_lines = true;

  // Draw checkerboard texture
  viewer.data().show_texture = true;


  std::cerr << "Press space for running an iteration." << std::endl;
  std::cerr << "Press 1 for Mesh 2 for UV" << std::endl;

  // Launch the viewer
  viewer.launch();
}

template <typename PointSource, typename PointTarget> int
TransformationEstimationJointOptimize<PointSource, PointTarget>::OptimizationFunctor::operator() (const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const
{
	const pcl::PointCloud<PointSource> & src_points = *estimator_->tmp_src_;
	const pcl::PointCloud<PointTarget> & tgt_points = *estimator_->tmp_tgt_;
	const std::vector<int> & src_indices = *estimator_->tmp_idx_src_;
	const std::vector<int> & tgt_indices = *estimator_->tmp_idx_tgt_;
	const std::vector<int> & src_indices_dfp = *estimator_->tmp_idx_src_dfp_;
	const std::vector<int> & tgt_indices_dfp = *estimator_->tmp_idx_tgt_dfp_;
	const std::vector<int> & src_indices_handles = *estimator_->tmp_idx_src_handles_;
	const std::vector<int> & tgt_indices_handles = *estimator_->tmp_idx_tgt_handles_;
	const float visualWeight = estimator_->visualFeatureWeight;
	const float denseWeight = estimator_->denseCloudWeight;
	const float handleWeight = estimator_->handleFeatureWeight;

	// Initialize the warp function with the given parameters
	Eigen::VectorXf params = x.cast<float> ();
	estimator_->warp_point_->setParam (params);

	Eigen::Matrix4f curr_transformation_matrix = estimator_->warp_point_->getTransform ();
	//std::cout << "[OptimizationFunctor::operator()] current transform = " << std::endl << curr_transformation_matrix << std::endl;

	//std::cerr << "[error function] before visual points: " << fvec.sum() << "\n";
	// Sum of squared distances of distinctive feature points
//	double diff_value_dfp = 0;
	const double dfp_factor = visualWeight/number_dfp;
	for (int i = 0; i < number_dfp; ++i)
	{
		std::cerr << "[error function] entered visual points loop: " << fvec.sum() << "\n";
		const PointSource & p_src = src_points.points[src_indices_dfp[i]];
		const PointTarget & p_tgt = tgt_points.points[tgt_indices_dfp[i]];

		// Transform the source point based on the current warp parameters
		PointSource p_src_warped;
		estimator_->warp_point_->warpPoint (p_src, p_src_warped);

		// Estimate the distance (cost function)
//		diff_value_dfp += estimator_->computeDistance (p_src_warped, p_tgt);
		fvec[i] = dfp_factor * estimator_->computeDistance (p_src_warped, p_tgt);
	}
	//std::cerr << "[error function] after visual points: " << fvec.sum() << "\n";

	const double p_factor = (denseWeight)/number_p;
	for (int i = 0; i < number_p; ++i)
	{
		const PointSource & p_src = src_points.points[src_indices[i]];
		const PointTarget & p_tgt = tgt_points.points[tgt_indices[i]];

		// Transform the source point based on the current warp parameters
		PointSource p_src_warped;
		estimator_->warp_point_->warpPoint (p_src, p_src_warped);

		// Estimate the distance (cost function)
//		diff_value_p += estimator_->computeDistance (p_src_warped, p_tgt);
		fvec[i+number_dfp] = p_factor * estimator_->computeDistancePointToPlane (p_src_warped, p_tgt);

		if(fvec[i+number_dfp] != fvec[i+number_dfp])
		{
			std::cerr << "i: " << i << " is nan! fvec[i]= " << fvec[i+number_dfp] << "\n";
			std::cerr << "source point       : " << p_src << "\n";
			std::cerr << "source point warped: " << p_src_warped << "\n";
			std::cerr << "target point       : " << p_tgt << "\n";
		}

		//	std::cerr << "fvec point " << (int)i << ":" << fvec[i+number_dfp] << "\n";
	}
	//std::cerr << "[error function] after points and features: " << fvec.sum() << "\n";

	const double handle_factor = (handleWeight)/number_handle_p;
	double totalError = 0.0;
	for (int i = 0; i < number_handle_p; ++i)
	{
		const PointSource & p_src = src_points.points[src_indices_handles[i]];
		const PointTarget & p_tgt = tgt_points.points[tgt_indices_handles[i]];

		// Transform the source point based on the current warp parameters
		PointSource p_src_warped;
		estimator_->warp_point_->warpPoint (p_src, p_src_warped);

		// Estimate the distance (cost function)
//		diff_value_p += estimator_->computeDistance (p_src_warped, p_tgt);
		fvec[i+number_dfp+number_p] = handle_factor * estimator_->computeDistance (p_src_warped, p_tgt);
		totalError += fvec[i+number_dfp+number_p];
		//std::cerr << "i: " << i << " src indice: " << src_indices_handles[i] << " src point: " << p_src_warped
		//                        << " tgt indice: " << tgt_indices_handles[i] << " tgt point: " << p_tgt
		//std::cerr << "i: " << i << " Distance error: " << fvec[i+number_dfp+number_p] << "\n";

		//std::cerr << "fvec handle point " << (int)i << ":" << fvec[i+number_dfp+number_p] << "\n";
	}
	//std::cerr << "[error function] after points and features and handles: " << fvec.sum() << "\n";
	//fvec[0] = totalError;
	//for (int i = 1; i < (number_handle_p + number_p + number_dfp); ++i)
	//{

	//	std::cerr << "i: " << i << "error: " << fvec[i] << "\n";
	//	fvec[i] = 0.0;
	//}
	//std::cerr << "[error functino] fvec total error : " << fvec.sum() << "\n";
	//std::cerr << "[error functino] my   total error : " << totalError << "\n";
	//std::cerr << "##############Other debug info: ##########\n"
//.........這裏部分代碼省略.........

	void state::compute(int want, FitContext *fc)
	{
		state *st = (state*) this;
		auto *oo = this;

		for (auto c1 : components) {
			if (c1->fitFunction) {
				omxFitFunctionCompute(c1->fitFunction, want, fc);
			} else {
				omxRecompute(c1, fc);
			}
		}
		if (!(want & FF_COMPUTE_FIT)) return;

		int nrow = components[0]->rows;
		for (auto c1 : components) {
			if (c1->rows != nrow) {
				mxThrow("%s: component '%s' has %d rows but component '%s' has %d rows",
					 oo->name(), components[0]->name(), nrow, c1->name(), c1->rows);
			}
		}

		Eigen::VectorXd expect;
		Eigen::VectorXd rowResult;
		int numC = components.size();
		Eigen::VectorXd tp(numC);
		double lp=0;
		for (int rx=0; rx < nrow; ++rx) {
			if (expectation->loadDefVars(rx) || rx == 0) {
				omxExpectationCompute(fc, expectation, NULL);
				if (!st->transition || rx == 0) {
					EigenVectorAdaptor Einitial(st->initial);
					expect = Einitial;
					if (expect.rows() != numC || expect.cols() != 1) {
						omxRaiseErrorf("%s: initial prob matrix must be %dx%d not %dx%d",
							       name(), numC, 1, expect.rows(), expect.cols());
						return;
					}
				}
				if (st->transition && (st->transition->rows != numC || st->transition->cols != numC)) {
					omxRaiseErrorf("%s: transition prob matrix must be %dx%d not %dx%d",
						       name(), numC, numC, st->transition->rows, st->transition->cols);
					return;
				}
			}
			for (int cx=0; cx < int(components.size()); ++cx) {
				EigenVectorAdaptor Ecomp(components[cx]);
				tp[cx] = Ecomp[rx];
			}
			if (st->verbose >= 4) {
				mxPrintMat("tp", tp);
			}
			if (st->transition) {
				EigenMatrixAdaptor Etransition(st->transition);
				expect = (Etransition * expect).eval();
			}
			rowResult = tp.array() * expect.array();
			double rowp = rowResult.sum();
			rowResult /= rowp;
			lp += log(rowp);
			if (st->transition) expect = rowResult;
		}
		oo->matrix->data[0] = Global->llScale * lp;
		if (st->verbose >= 2) mxLog("%s: fit=%f", oo->name(), lp);
	}