C++ VectorXf::size方法代码示例

void Layouter::laplacianEigMap( const MatrixXd& laplacian, const VectorXf& radiusVec, const VectorXi& hashID, MatrixXf& finalPos2D, float& finalRadius, float sparseFactor /*= 1.f*/ )
{
	if (laplacian.rows() <= 0 || laplacian.cols() <= 0 || radiusVec.size() <= 0)
	{
		finalPos2D = MatrixXf::Zero(1,2);
		finalRadius = 0;
		return;
	}
	if (laplacian.rows() == 1 && laplacian.cols() == 1 && radiusVec.size() == 1)
	{
		finalPos2D = MatrixXf::Zero(1,2);
		finalRadius = radiusVec[0];
		return;
	}
	MatrixXd pos2D, finalPosd;
	LaplacianSolver::compute(laplacian, pos2D); 

	VectorXd minPos, maxPos;
	MDSPostProcesser m_postProcessor(500, 1.0f, 1.0, 1.02, radiusVec.minCoeff());
	m_postProcessor.setSparseFactor(sparseFactor);
	m_postProcessor.set2DPos(pos2D, radiusVec.cast<double>(), &hashID);
	m_postProcessor.compute();
	m_postProcessor.getFinalPos(finalPosd);
	m_postProcessor.getFinalBoundingRect(minPos, maxPos);
	finalPos2D = finalPosd.cast<float>();
	// 	finalPos2D = pos2D.cast<float>();
	// 	minPos = pos2D.colwise().minCoeff();
	// 	maxPos = pos2D.colwise().maxCoeff();
	//  VectorXd size = maxPos - minPos;
	finalRadius = m_postProcessor.getFinalRadius();// size.norm() * 0.5;//(size[0] > size[1] ? size[0] : size[1]) * 0.5f;
}

float LinearInterpolation (VectorXf X , VectorXf Y, float X_PointOfInterest){
//Produce Y_point_of_interest given X,Y and target X_point_of_interest
//X : vector containing the X variables of the interpolant 
//Y : vector containing the Y variables of the interpolant 
//PointOfInterest : Point of X to estimate the new point of Y
    
  float   xk, xkp1, yk, ykp1;  //Points adjecent to the point of interpolation
  if ( X.size() != Y.size() ){cout << "Problem with vector sizes" << endl; return(-1);}
//cout <<  " X(0): " <<  X(0) <<" X(Y.size()-1): " <<X(Y.size()-1)   <<   " Point of interest: " << X_PointOfInterest<< endl;
  if ( X_PointOfInterest < X(0) || X_PointOfInterest > X(Y.size()-1) ){cout << "You interpolate out of the curve boundaries" << endl; return(-1);}
//Find the points right before and right after the point of interest
  for (int i=1; i<X.size() ; i++){
    if (X(i)>= X_PointOfInterest){
      xkp1 = X(i);
      xk = X(i-1);
      ykp1 = Y(i);
      yk = Y(i-1);
      break;}
  }
//point-slope form for a line formula
  float t = (X_PointOfInterest -xk)/(xkp1 -xk);
  float yPOI = (1-t) * yk + t * ykp1;  // estimate point of interest
// cout << "(" << xk << ",  " << X_PointOfInterest << " , " << xkp1 << ") & (" << yk << ", " << yPOI  << ", " << ykp1 << ")"<< endl;
  return (yPOI);   
} 

double vector_stdv_mad( VectorXf residues)
{
    // Return the standard deviation of vector with MAD estimation
    int n_samples = residues.size();
    sort( residues.derived().data(),residues.derived().data()+residues.size());
    double median = residues( n_samples/2 );
    residues << ( residues - VectorXf::Constant(n_samples,median) ).cwiseAbs();
    sort(residues.derived().data(),residues.derived().data()+residues.size());
    double MAD = residues( n_samples/2 );
    return 1.4826 * MAD;
}

MatrixX3f Surface::compute_normals(const MatrixX3f& rr, const MatrixX3i& tris)
{
    printf("\tcomputing normals\n");
    // first, compute triangle normals
    MatrixX3f r1(tris.rows(),3); MatrixX3f r2(tris.rows(),3); MatrixX3f r3(tris.rows(),3);

    for(qint32 i = 0; i < tris.rows(); ++i)
    {
        r1.row(i) = rr.row(tris(i, 0));
        r2.row(i) = rr.row(tris(i, 1));
        r3.row(i) = rr.row(tris(i, 2));
    }

    MatrixX3f x = r2 - r1;
    MatrixX3f y = r3 - r1;
    MatrixX3f tri_nn(x.rows(),y.cols());
    tri_nn.col(0) = x.col(1).cwiseProduct(y.col(2)) - x.col(2).cwiseProduct(y.col(1));
    tri_nn.col(1) = x.col(2).cwiseProduct(y.col(0)) - x.col(0).cwiseProduct(y.col(2));
    tri_nn.col(2) = x.col(0).cwiseProduct(y.col(1)) - x.col(1).cwiseProduct(y.col(0));

    //   Triangle normals and areas
    MatrixX3f tmp = tri_nn.cwiseProduct(tri_nn);
    VectorXf normSize = tmp.rowwise().sum();
    normSize = normSize.cwiseSqrt();

    for(qint32 i = 0; i < normSize.size(); ++i)
        if(normSize(i) != 0)
            tri_nn.row(i) /= normSize(i);

    MatrixX3f nn = MatrixX3f::Zero(rr.rows(), 3);

    for(qint32 p = 0; p < tris.rows(); ++p)
    {
        Vector3i verts = tris.row(p);
        for(qint32 j = 0; j < verts.size(); ++j)
            nn.row(verts(j)) = tri_nn.row(p);
    }

    tmp = nn.cwiseProduct(nn);
    normSize = tmp.rowwise().sum();
    normSize = normSize.cwiseSqrt();

    for(qint32 i = 0; i < normSize.size(); ++i)
        if(normSize(i) != 0)
            nn.row(i) /= normSize(i);

    return nn;
}

VectorXf AutoEncoder::Calculate(int index)
{
    /*
     *  Description:
     *  Calculate the i-th samples by feedforward and store hiddenvector
     *  in the row i of hidden matrix, output in row i of output matrix
     *
     *  @return outputVector: The output of FF
     */
    
    VectorXf HiddenVector = Weight_encode * OriginalData->row(index).transpose() + Bias_encode;
    for (int i = 0; i < HiddenVector.size(); i++)
    {
        HiddenVector(i) = sigmoid(HiddenVector(i));
    }
    HiddenMatrix.row(index) = HiddenVector.transpose();

    VectorXf output_vector = VectorXf(IO_dim);
    output_vector = Weight_decode * HiddenVector + Bias_decode;
    for (int i = 0; i < output_vector.size(); i++)
    {
        output_vector(i) = sigmoid(output_vector(i));
    }
    OutputMatrix.row(index) = output_vector.transpose();
    
    return output_vector;
}

vector<int> Util::descending_order(VectorXf& values) {
    vector<int> indices(values.size());
    for (int i = 0; i != indices.size(); ++i) indices[i] = i;
    DescentCompareIndicesByAnotherVectorValues comp(values);
    sort(indices.begin(), indices.end(), comp);
    return indices;
}

float sparsify_absolute (
    const VectorXf & dense,
    VectorSf & sparse,
    const Vector<float> & thresh)
{
  ASSERT1_LE(0, min(thresh));

  const int I = dense.size();

  sparse.resize(I);

  float loss = 0;

  for (int i = 0; i < I; ++i) {

    const float dense_i = dense(i);
    const float abs_i = fabsf(dense_i);

    if (abs_i > thresh[i]) sparse.insert(i) = dense_i;
    else loss += abs_i;
  }

  sparse.finalize();

  return loss;
}

void NonRigid::computeInflateDir(VectorXf &p_vec, VectorXf &g_vec)
{
    VectorXf arap_g = g_vec;
    VectorXf inflate_g = g_vec;
    VectorXf userCrsp_g = g_vec;
    double f_e = 0;
    computeArap(p_vec, arap_g);
    computeUserCrsp(p_vec, userCrsp_g);

    Vector3f cur_v(0,0,0);
    Vector3f cur_g(0,0,0);
    size_t P_Num = p_vec.size()/3;
    for (size_t i = 0; i < P_Num; ++i)
    {
        cur_v << p_vec(i + 0*P_Num), p_vec(i + 1*P_Num), p_vec(i + 2*P_Num);
        if (computeVertexGrad(cur_v, cur_g) > 0.0)
        {
            inflate_g(i + 0*P_Num) = cur_g(0);
            inflate_g(i + 1*P_Num) = cur_g(1);
            inflate_g(i + 2*P_Num) = cur_g(2);
        }
        else
        {
            inflate_g(i + 0*P_Num) = -P_Prime_N(0, i);
            inflate_g(i + 1*P_Num) = -P_Prime_N(1, i);
            inflate_g(i + 2*P_Num) = -P_Prime_N(2, i);
        }

    }

    g_vec = lamd_inflate*inflate_g + lamd_arap*arap_g + lamd_userCrsp*userCrsp_g;
    g_vec.normalize();
}

int MonotonicityCheck(VectorXf X){
   // evaluate vector monotonicity of vector X
   int N =  X.size(); 
   for (int i=0; i<(N-1) ; i++){  
     if ( X(i) >X (1+i) ) return (-1);   }
   return (1);
}

VectorXf fnvalspapi(VectorXf X, VectorXf Y, VectorXf X_target){
  //evaluate Y_target for X_target given X and Y
    
	int N =  X_target.size();
	VectorXf rr(N); 
	for (int i=0; i<N ; i++){  
      rr(i) =  LinearInterpolation(X, Y, X_target(i)) ;  }
	return(rr);
}

VectorXf NewSolution( VectorXf x0 , float Step, int point, VectorXf t_reg){ 
    //generate new solution on the simplex defined in [t_reg(0)-x0-t_reg(N-1)] using a displacement of size Step
    // x0   : initial solution
    // Step : displacement size
    // point: knot to perturb
    // t_reg: time_grid
    
    VectorXf InitConf (2+ x0.size());
    InitConf(0) = t_reg(0); InitConf( x0.size()+1) = t_reg.maxCoeff()  ; InitConf.segment(1, x0.size()) = x0; 
    float LowBou =  InitConf(point-1);
    float UppBou =  InitConf(point+1);
    float New_State = (UppBou - LowBou) * Step + LowBou; 
    InitConf(point) = New_State; 
    //if (  MonotonicityCheck( InitConf.segment(1, x0.size()) ) == -1) { 
    //    cout << "We generated a unacceptable solutiion" << endl << "Initial seed was : " << x0.transpose() 
    //         << endl << " and we produced :"<<InitConf.segment(1, x0.size()).transpose() << endl;}
    return (InitConf.segment(1, x0.size()) ) ;
  } 

void write_to_python (const VectorXf & x, ostream & os)
{
  const int entries_per_line = 8;

  os << "[";
  for (int i = 0; i < x.size(); ++i) {
    os << x[i] << ", ";
    if ((i + 1) % entries_per_line == 0) os << "\n  ";
  }
  os << "]";
}

	virtual float bestThreshold( const VectorXf & f, float * gain ) const {
		const int N = lbl_.maxCoeff()+1;
		const float EPS=1e-6;
		// Create the feature/label pairs
		std::vector< std::pair<float,int> > elements( f.size() );
		for( int i=0; i<f.size(); i++ )
			elements[i] = std::make_pair( f[i], i );
		std::sort(elements.begin(), elements.end() );
		
		// And compute the probabilities and cirterion
		ArrayXf wl = 1e-20*ArrayXf::Ones(N), wr = 1e-20*ArrayXf::Ones(N);
		for( int i=0; i<lbl_.size(); i++ )
			wr[ lbl_[i] ] += weight_[i];
		
		// Initialize the thresholds
		float best_gain = 0, tbest = (elements.front().first+elements.back().first)/2, last_t = elements.front().first;
		const float tot_s = score( wr/wr.sum() );
		// Find the best threshold
		for( auto i: elements ) {
			const float t = i.first;
			const int j = i.second;
			// If there is a threshold
			if( t - last_t > EPS ) {
				// Compute the score
				const float l = wl.sum(), r = wr.sum();
				const float sl = score(wl/l), sr = score(wr/r);
				const float g = tot_s - ( sl*l/(l+r) + sr*r/(l+r) );
				if( g > best_gain ) {
					best_gain = g;
					tbest = (last_t+t)/2.;
				}
			}
			// Update the probabilities
			wl[ lbl_[j] ] += weight_[j];
			wr[ lbl_[j] ] -= weight_[j];
			last_t = t;
		}
		if( gain )
			*gain = best_gain;
		return tbest;
	}

float rttemp1( VectorXf t_reg,  VectorXf curvei, VectorXf curvek,int nknots, float lambda, VectorXf initial){
  //Calculate the cost of the given warping
  // t_reg  : time grid of y_reg
  // curvei : query curve
  // curvek : reference curve
  // nknots : number of knots
  // lambda : time distortion penalty parameter
  // initial: position of knots on the simplex 
    
  int N = t_reg.size();
  VectorXf struct_ = VectorXf::LinSpaced( nknots+2, t_reg(0)  , t_reg.maxCoeff() ); 
  VectorXf hik(N);
  VectorXf Q(2+ initial.size()) ;           //Solution with the placement of the knots on the simplex
  Q(0) = t_reg(0) ; Q(1+ initial.size()) = t_reg.maxCoeff()  ; Q.segment(1, initial.size()) = initial;
  hik =  fnvalspapi( struct_, Q , t_reg);   // compute the new internal time-scale  
  
  //cout << "hik: " << hik.transpose() << endl;
  //cout << "Monotonicity Checked on Hik: " << MonotonicityCheck(hik) << endl;
  //if(  MonotonicityCheck(hik) == -1) { cout <<" Q.transpose()  is :"<< Q.transpose()  << endl;}
  
  return ( (fnvalspapi(t_reg,curvei,hik)-fnvalspapi(t_reg,curvek ,t_reg)).array().pow(2).sum() + lambda * (hik - t_reg).array().pow(2).sum() );
}

rthink_Output rthik_SA(VectorXf t_reg,  VectorXf curvei, VectorXf curvek,int nknots, float lambda){
  // Random Search solver for the optimization problem of pairwise warping
  // t_reg : time_grid
  // curvei: query curve
  // curvek: reference curve 
  // nknots : number of knots
  // lambda : time distortion penalty parameter
    
  int k=0; float OldSol, bk; 
  VectorXf xk(nknots);   
  bk = (t_reg.maxCoeff() - t_reg(0))/(1+ float(nknots ));   //Distance between adjacent knots and edges-knots
  xk  = VectorXf::LinSpaced(nknots , t_reg(0) + bk  ,  - bk + t_reg.maxCoeff()  ); //Initial candidate solution with equispaced knots 
  VectorXf xn(nknots);   VectorXf help(2+nknots); 
 
  float NewSol; 
  OldSol =  rttemp1(t_reg, curvei, curvek, nknots , lambda, xk);  //Cost of initial solution
  int z= 99*nknots;                                               //Number of random search to do (proportional to the # of knots)
  //srand(1);                                                     //Fix the seed to have reproducable behaviour
  VectorXf Steps(z); Steps = (ArrayXf::Random(z)+1.)/2.;          //Generate possible random pertubations magnitude
  VectorXi Posit(z); for (int u=0; u <z; u++ ) Posit(u) =1+ rand()%(nknots+0); //Generate list of positions to purturb
 
  k=0;
  while((OldSol > .0001) && (k<z)) {
    xn = NewSolution(xk,Steps(k), Posit(k), t_reg );              //Get a new solution
    NewSol = rttemp1(t_reg, curvei, curvek, nknots, lambda, xn);  //Cost of new solution
    if (  (NewSol < OldSol)  ) {                                  //If it's better than the old one, use it.
      OldSol= NewSol; xk= xn; 
    }
  k++;
}

  VectorXf x3 =  VectorXf::LinSpaced(2+nknots, t_reg(0), t_reg( t_reg.size()-1)); 
  help(0) =  t_reg(0) ; help(nknots+1) = t_reg( t_reg.size()-1) ; help.segment(1,nknots) = xk;  

  rthink_Output G;
  G.Val = OldSol;
  G.Mapping = fnvalspapi(  x3 , help  , t_reg);
  return G ;
}