C++ MatrixXd::sum方法代码示例

void ActivationFunctionsTestCase::linear()
{
  const int N = 1000;
  Eigen::MatrixXd a = Eigen::VectorXd::Random(N) * 10.0;
  Eigen::MatrixXd z = Eigen::VectorXd::Zero(N);
  OpenANN::linear(a, z);
  ASSERT_EQUALS(a.minCoeff(), z.minCoeff());
  ASSERT_EQUALS(a.maxCoeff(), z.maxCoeff());

  Eigen::MatrixXd gd = Eigen::VectorXd::Zero(N);
  Eigen::MatrixXd expected = Eigen::VectorXd::Ones(N);
  OpenANN::linearDerivative(gd);
  ASSERT_EQUALS(gd.sum(), expected.sum());
}

Eigen::MatrixXd caculateAlgebraicDistance(const Eigen::SparseMatrix<int,Eigen::RowMajor>& incidenceMatrix,
                                          const Eigen::SparseMatrix<double,Eigen::RowMajor>& updatedMatrix)
{
  Eigen::MatrixXd inc = Eigen::MatrixXd(incidenceMatrix);
  Eigen::MatrixXd update = Eigen::MatrixXd(updatedMatrix);
  int edgeNumber = inc.cols();
  int verticeNumber = inc.rows();
  Eigen::MatrixXd algebraicDistanceForEdge(edgeNumber,1);
  Eigen::VectorXd sumOfcolsVector(edgeNumber);
  sumOfcolsVector = inc.colwise().sum();
  for(int i = 0; i<edgeNumber;++i){
    int verticeNumberForEdge = sumOfcolsVector[i];
    Eigen::MatrixXd allEdgeMatrix(update.rows(),verticeNumberForEdge);
    int g=0;
    for(int j=0;j<verticeNumber;++j){
        if(inc(j,i)){
          ++g;
          for(int v = 0; v<update.rows();++v){
              allEdgeMatrix(v,g-1) = update(v,j);
          }
        }
    }

    if(g!=0){
      Eigen::MatrixXd maxValueOfEdge = getMaxDiffForeachEdge(allEdgeMatrix);
      maxValueOfEdge = maxValueOfEdge.transpose() * maxValueOfEdge;
      double edgeDistance = maxValueOfEdge.sum();
      edgeDistance = sqrt(edgeDistance);
      algebraicDistanceForEdge(i,0) = edgeDistance;

    }

  }
  return algebraicDistanceForEdge;
}

void ActivationFunctionsTestCase::rectifier()
{
  const int N = 1000;
  Eigen::MatrixXd a = Eigen::MatrixXd::Random(1, N) * 10.0;
  Eigen::MatrixXd z = Eigen::MatrixXd::Zero(1, N);
  OpenANN::rectifier(a, z);
  ASSERT_EQUALS(0.0, z.minCoeff());
  ASSERT_EQUALS(a.maxCoeff(), z.maxCoeff());

  Eigen::MatrixXd gd = Eigen::MatrixXd::Zero(1, N);
  Eigen::MatrixXd expected = Eigen::MatrixXd::Ones(1, N);
  for(int i = 0; i < N; i++)
    expected(i) *= (double)(z(i) > 0.0);
  OpenANN::rectifierDerivative(z, gd);
  ASSERT_EQUALS(gd.sum(), expected.sum());
}

void ActivationFunctionsTestCase::softmax()
{
  const int N = 1000;
  Eigen::MatrixXd a = Eigen::VectorXd::Random(N).transpose();
  OpenANN::softmax(a);
  ASSERT_EQUALS_DELTA(1.0, a.sum(), 1e-3);
  ASSERT_WITHIN(a.minCoeff(), 0.0, 1.0);
  ASSERT_WITHIN(a.maxCoeff(), 0.0, 1.0);
}

Eigen::MatrixXd localWeighting( const Eigen::MatrixXd &W, bool isFull, bool isSymmetric)
{
	int n = W.rows();
	double Ls = (W.count()-n)/2; //count number of edges of the graph

	Eigen::MatrixXd C = Eigen::MatrixXd::Zero(n,n);
	double Ws = 0.5*W.sum();

	Eigen::VectorXd D = W.colwise().sum();

	if (isFull)
	{
		if (isSymmetric)
		{
			computeC_symmetric_full(const_cast<Eigen::MatrixXd &>(W),C,D,n);
		}
		else
		{
			// this is a trick to ensure vectorizatoin and no cache misses! some tranpositions have to be made though
			const_cast<Eigen::MatrixXd &>(W).transposeInPlace();
			computeC_symmetric_full(const_cast<Eigen::MatrixXd &>(W),C,D,n);
			const_cast<Eigen::MatrixXd &>(W).transposeInPlace();
			C.transposeInPlace();
			// the original code use vertical access to rows, but is slower
			//compute_C_asymmetric_full(const_cast<Eigen::MatrixXd &>(W),C,D,n);
		}
	}
	else
	{
		if (isSymmetric)
		{
			computeC_symmetric_sparse(const_cast<Eigen::MatrixXd &>(W),C,D,n);
		}
		else
		{
			compute_C_asymmetric_sparse(const_cast<Eigen::MatrixXd &>(W),C,D,n);
		}
	}

	Eigen::MatrixXd G = ((Ls/Ws)*W).cwiseProduct(C);

	Eigen::VectorXd DG = G.colwise().sum();

	for (int i=0; i<n; i++)
	{
		G.row(i)/=DG(i);
	}
	
	return G;
}

 inline Real stdDev(const Eigen::MatrixXd& input){
   return ::sqrt(((Eigen::MatrixXd)((input.array()-input.sum()/input.rows()).pow(2.0))).sum()/(input.rows()-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++)
    {
//.........这里部分代码省略.........

/** unittest for oapackage
 *
 * Returns UNITTEST_SUCCESS if all tests are ok.
 *
 */
int oaunittest (int verbose, int writetests = 0, int randval = 0) {
        double t0 = get_time_ms ();
        const char *bstr = "OA unittest";
        cprintf (verbose, "%s: start\n", bstr);

        srand (randval);

        int allgood = UNITTEST_SUCCESS;

        Combinations::initialize_number_combinations (20);

        /* constructors */
        {
                cprintf (verbose, "%s: interaction matrices\n", bstr);

                array_link al = exampleArray (2);
                Eigen::MatrixXd m1 = array2xfeigen (al);
                Eigen::MatrixXd m2 = arraylink2eigen (array2xf (al));

                Eigen::MatrixXd dm = m1 - m2;
                int sum = dm.sum ();

                myassert (sum == 0, "unittest error: construction of interaction matrices\n");
        }

		cprintf(verbose, "%s: reduceConferenceTransformation\n", bstr);
		myassert(unittest_reduceConferenceTransformation()==0, "unittest unittest_reduceConferenceTransformation failed");

        /* constructors */
        {
                cprintf (verbose, "%s: array manipulation operations\n", bstr);

                test_array_manipulation (verbose);
        }

        /* double conference matrices */
        {
                cprintf (verbose, "%s: double conference matrices\n", bstr);

                array_link al = exampleArray (36, verbose);
                myassert (al.is_conference (2), "check on double conference design type");

                myassert (testLMC0checkDC (al, verbose >= 2), "testLMC0checkDC");

        }

        /* conference matrices */
        {
                cprintf (verbose, "%s: conference matrices\n", bstr);

                int N = 4;
                conference_t ctype (N, N, 0);

                arraylist_t kk;
                array_link al = ctype.create_root ();
                kk.push_back (al);

                for (int extcol = 2; extcol < N; extcol++) {
                        kk = extend_conference (kk, ctype, 0);
                }
                myassert (kk.size () == 1, "unittest error: conference matrices for N=4\n");
        }

        {
                cprintf (verbose, "%s: generators for conference matrix extensions\n", bstr);
                test_conference_candidate_generators (verbose);
        }

        {
                cprintf (verbose, "%s: conference matrix Fvalues\n", bstr);
                array_link al = exampleArray (22, 0);
                if (verbose >= 2)
                        al.show ();
                if (0) {
                        std::vector< int > f3 = al.FvaluesConference (3);
                        if (verbose >= 2) {
                                printf ("F3: ");
                                display_vector (f3);
                                printf ("\n");
                        }
                }

                const int N = al.n_rows;
                jstructconference_t js (N, 4);
                std::vector< int > f4 = al.FvaluesConference (4);
                std::vector< int > j4 = js.Jvalues ();

                if (verbose >= 2) {
                        printf ("j4: ");
                        display_vector (j4);
                        printf ("\n");
                        printf ("F4: ");
                        display_vector (f4);
                        printf ("\n");
                }
//.........这里部分代码省略.........

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