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

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;
}

void KF_joseph_update(VectorXf &x, MatrixXf &P,float v,float R, MatrixXf H)
{
    VectorXf PHt = P*H.transpose();
    MatrixXf S = H*PHt;
    S(0,0) += R;
    MatrixXf Si = S.inverse();
    Si = make_symmetric(Si);
    MatrixXf PSD_check = Si.llt().matrixL(); //chol of scalar is sqrt
    PSD_check.transpose();
    PSD_check.conjugate();

    VectorXf W = PHt*Si;
    x = x+W*v;
    
    //Joseph-form covariance update
    MatrixXf eye(P.rows(), P.cols());
    eye.setIdentity();
    MatrixXf C = eye - W*H;
    P = C*P*C.transpose() + W*R*W.transpose();  

    float eps = 2.2204*pow(10.0,-16); //numerical safety 
    P = P+eye*eps;

    PSD_check = P.llt().matrixL();
    PSD_check.transpose();
    PSD_check.conjugate(); //for upper tri
}

double NonRigid::computeArap(VectorXf &p_vec, VectorXf &g_vec)
{
    VectorXf d_vec(VectorXf::Map(d_cur.data(), d_cur.cols()*d_cur.rows()));

    g_vec = L*p_vec - d_vec;

    return (0.5*p_vec.transpose()*L*p_vec - d_vec.transpose()*p_vec)(0, 0);
}

MatrixXf Sphere::make_initial_simplex( const VectorXf &pars, float size )
{
    /*
    * Make the initial tetrahedron
    */
    int npar = pars.size();

    MatrixXf simplex = MatrixXf::Zero(npar+1,npar);

    simplex.rowwise() += pars.transpose();

    for (int k = 1; k < npar+1; k++) {
        simplex(k,k-1) += size;
    }

    return simplex;
}

void drwnNNGraphLearner::nearestNeighbourUpdate(unsigned nCycle, unsigned maxIterations)
{
    DRWN_FCN_TIC;

#ifndef USE_THREADING
    // nearest neighbour search
    for (unsigned nIterations = 0; nIterations < maxIterations; nIterations++) {

        // perform an update
        drwnNNGraphMoves::update(_posGraph, drwnNNGraphLabelsEqualMetric());
        drwnNNGraphMoves::update(_negGraph, drwnNNGraphLabelsNotEqualMetric());

        // check energy
        pair<double, double> e_pos = _posGraph.energy();
        pair<double, double> e_neg = _negGraph.energy();
        DRWN_LOG_MESSAGE("...search iteration " << nCycle << "." << nIterations
            << "; energy (+) " << e_pos.first << " (-) " << e_neg.first);
    }
#else
    drwnNNGraphMoveUpdateThread<drwnNNGraphLabelsEqualMetric> posThread(&_posGraph, drwnNNGraphLabelsEqualMetric());
    drwnNNGraphMoveUpdateThread<drwnNNGraphLabelsNotEqualMetric> negThread(&_negGraph, drwnNNGraphLabelsNotEqualMetric());

    // threaded nearest neighbour search
    {
        drwnThreadPool threadPool(2);
        for (unsigned nIterations = 0; nIterations < maxIterations; nIterations++) {

            // perform an update
            threadPool.start();
            threadPool.addJob(&posThread);
            threadPool.addJob(&negThread);
            threadPool.finish();

            // report energy
            DRWN_LOG_MESSAGE("...search iteration " << nCycle << "." << nIterations
                << "; energy (+) " << posThread.energy.first << " (-) " << negThread.energy.first);
        }
    }
#endif

    // cache data matrix
    DRWN_LOG_VERBOSE("...caching data matrix of size " << _dim << "-by-" << _dim);
    _X = MatrixXd::Zero(_dim, _dim);
#ifndef USE_THREADING
    drwnNNGraphNodeIndex u(0, 0);
    for (u.imgIndx = 0; u.imgIndx < _posGraph.numImages(); u.imgIndx++) {
        MatrixXf localX = MatrixXf::Zero(_dim, _dim);
        for (u.segId = 0; u.segId < _posGraph[u.imgIndx].numNodes(); u.segId++) {
            const drwnNNGraphEdgeList& e = _posGraph[u].edges;
            for (drwnNNGraphEdgeList::const_iterator kt = e.begin(); kt != e.end(); ++kt) {
                const drwnNNGraphNodeIndex v(kt->targetNode);
                const VectorXf delta(_graph[u].features - _graph[v].features);
                localX += delta * delta.transpose();
            }
        }
        _X += localX.cast<double>();
    }
#else
    {
        // prepare thread data
        const unsigned nJobs = std::min((unsigned)_posGraph.numImages(),
            std::max((unsigned)1, drwnThreadPool::MAX_THREADS));
        vector<set<unsigned> > imgIndxes(nJobs);
        for (unsigned imgIndx = 0; imgIndx < _posGraph.numImages(); imgIndx++) {
            imgIndxes[imgIndx % nJobs].insert(imgIndx);
        }

        // start threads
        drwnThreadPool threadPool(nJobs);
        threadPool.start();
        vector<drwnNNGraphDataMatrixThread *> jobs(nJobs);
        for (unsigned i = 0; i < nJobs; i++) {
            jobs[i] = new drwnNNGraphDataMatrixThread(this, &_X, imgIndxes[i]);
            threadPool.addJob(jobs[i]);
        }
        threadPool.finish();

        for (unsigned i = 0; i < jobs.size(); i++) {
            delete jobs[i];
        }
        jobs.clear();
    }
#endif

    DRWN_FCN_TOC;
}

float Sphere::opt_rad(const VectorXf &r0,const fitUserNew user)
{
  MatrixXf diff = user->rr.rowwise() - r0.transpose();
  return diff.rowwise().norm().mean();
}

void Sphere::calculate_cm_ave_dist (const MatrixXf &rr, VectorXf &cm, float &avep)
{
    cm = rr.colwise().mean();
    MatrixXf diff = rr.rowwise() - cm.transpose();
    avep = diff.rowwise().norm().mean();
}

inline
void r_and_t(MatrixXf &rot_cw, VectorXf &pos_cw,MatrixXf start_points, MatrixXf end_points,
             MatrixXf P1w,MatrixXf P2w,MatrixXf initRot_cw,VectorXf initPos_cw,
             int maxIterNum,float TerminateTh,int nargin)
{
    if(nargin<6)
    {
        return;
    }

    if(nargin<8)
    {
        maxIterNum = 8;
        TerminateTh = 1e-5;
    }

    int n = start_points.cols();

    if(n != end_points.cols() || n!= P1w.cols() || n!= P2w.cols())
    {
        return;
    }

    if(n<4)
    {
        return;
    }

    //first compute the weight of each line and the normal of
    //the interpretation plane passing through to camera center and the line

    VectorXf w = VectorXf::Zero(n);
    MatrixXf nc = MatrixXf::Zero(3,n);

    for(int i = 0 ; i < n ; i++)
    {
        //the weight of a line is the inverse of its image length
        w[i] = 1/(start_points.col(i)-end_points.col(i)).norm();
        vfloat3 v1 = start_points.col(i);
        vfloat3 v2 = end_points.col(i);
        vfloat3 temp = v1.cross(v2);
        nc.col(i) = temp/temp.norm();
    }

    MatrixXf rot_wc = initPos_cw.transpose();
    MatrixXf pos_wc = - initRot_cw.transpose() * initPos_cw;


    for(int iter = 1 ; iter < maxIterNum ; iter++)
    {
        //construct the equation (31)
        MatrixXf A = MatrixXf::Zero(6,7);
        MatrixXf C = MatrixXf::Zero(3,3);
        MatrixXf D = MatrixXf::Zero(3,3);
        MatrixXf F = MatrixXf::Zero(3,3);
        vfloat3 c_bar = vfloat3(0,0,0);
        vfloat3 d_bar = vfloat3(0,0,0);
        for(int i = 0 ; i < n ; i++)
        {
            //for first point on line
            vfloat3 Pi = rot_wc * P1w.col(i);
            vfloat3 Ni = nc.col(i);
            float wi = w[i];
            vfloat3 bi = Pi.cross(Ni);
            C = C + wi*Ni*Ni.transpose();
            D = D + wi*bi*bi.transpose();
            F = F + wi*Ni*bi.transpose();
            vfloat3 tempi = Pi + pos_wc;
            float scale = Ni.transpose() * tempi;
            scale *= wi;
            c_bar = c_bar + scale * Ni;
            d_bar = d_bar + scale*bi;
            //for second point on line
            Pi = rot_wc * P2w.col(i);
            Ni = nc.col(i);
            wi = w[i];
            bi = Pi.cross(Ni);
            C  = C + wi*Ni*Ni.transpose();
            D  = D + wi*bi*bi.transpose();
            F  = F + wi*Ni*bi.transpose();
            scale = (Ni.transpose() * (Pi + pos_wc));
            scale *= wi;
            c_bar = c_bar + scale * Ni;
            d_bar = d_bar + scale * bi;
        }
        A.block<3,3>(0,0) = C;
        A.block<3,3>(0,3) = F;
        (A.col(6)).segment(0,2) = c_bar;
        A.block<3,3>(3,0) = F.transpose();
        A.block<3,3>(2,2) = D;
        (A.col(6)).segment(3,5) = d_bar;
        //sovle the system by using SVD;
        JacobiSVD<MatrixXf> svd(A, ComputeThinU | ComputeThinV);
        VectorXf vec(7);
        //the last column of Vmat;
        vec = (svd.matrixV()).col(6);
        //the condition that the last element of vec should be 1.
        vec = vec/vec[6];

        //update the rotation and translation parameters;
//.........这里部分代码省略.........

int main(void)
{
  cout << "Eigen v" << EIGEN_WORLD_VERSION << "." << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION << endl;
  static const int R = 288;
  static const int N = R*(R+1)/2;
  static const int M = 63;
  static const int HALF_M = M/2;
  static const float nsigma = 2.5f;

  MatrixXf data = MatrixXf::Random(M, N);
  MatrixXf mask = MatrixXf::Zero(M, N);
  MatrixXf result = MatrixXf::Zero(1, N);
  VectorXf std = VectorXf::Zero(N);
  VectorXf centroid = VectorXf::Zero(N);
  VectorXf mean = VectorXf::Zero(N);
  VectorXf minval = VectorXf::Zero(N);
  VectorXf maxval = VectorXf::Zero(N);

  cout << "computing..." << flush;
  double t = GetRealTime();

  // computes the exact median
  if (M&1)
  {
#pragma omp parallel for
    for (int i = 0; i < N; i++)
    {
      vector<float> row(data.data()+i*M, data.data()+(i+1)*M);
      nth_element(row.begin(), row.begin()+HALF_M, row.end());
      centroid(i) = row[HALF_M];
    }
  }
  // nth_element guarantees x_0,...,x_{n-1} < x_n
  else
  {
#pragma omp parallel for
    for (int i = 0; i < N; i++)
    {
      vector<float> row(data.data()+i*M, data.data()+(i+1)*M);
      nth_element(row.begin(), row.begin()+HALF_M, row.end());
      centroid(i) = row[HALF_M];
      centroid(i) += *max_element(row.begin(), row.begin()+HALF_M);
      centroid(i) *= 0.5f;
    }
  }

  // compute the mean
  mean = data.colwise().mean();

  // compute std (x) = sqrt ( 1/N SUM_i (x(i) - mean(x))^2 )
  std = (((data.rowwise() - mean.transpose()).array().square()).colwise().sum() *
         (1.0f / M))
            .array()
            .sqrt();

  // compute n sigmas from centroid
  minval = centroid - std * nsigma;
  maxval = centroid + std * nsigma;
  
  // compute clip mask
  for (int i = 0; i < N; i++)
  {
    mask.col(i) = (data.col(i).array() > minval(i)).select(VectorXf::Ones(M), 0.0f);
    mask.col(i) = (data.col(i).array() < maxval(i)).select(VectorXf::Ones(M), 0.0f);
  }

  // apply clip mask to data
  data.array() *= mask.array();

  // compute mean such that we ignore clipped data, this is our final result
  result = data.colwise().sum().array() / mask.colwise().sum().array();

  t = GetRealTime() - t;
  cout << "[done]" << endl << endl;

  size_t bytes = data.size()*sizeof(float);
  cout << "data: " << M << "x" << N << endl;
  cout << "size: " << bytes*1e-6f << " MB" << endl;
  cout << "rate: " << bytes/(1e6f*t) << " MB/s" << endl;
  cout << "time: " << t << " s" << endl;

  return 0;
}

int Pca::Calculate(vector<float> &x, 
              const unsigned int &nrows, 
              const unsigned int &ncols, 
              const bool is_corr, 
              const bool is_center, 
              const bool is_scale) {
  _ncols = ncols;
  _nrows = nrows;
  _is_corr = is_corr;
  _is_center = is_center;
  _is_scale = is_scale;
  if (x.size()!= _nrows*_ncols) {
    return -1;
  }
  if ((1 == _ncols) || (1 == nrows)) {
    return -1;
  }
  // Convert vector to Eigen 2-dimensional matrix
  //Map<MatrixXf> _xXf(x.data(), _nrows, _ncols);
  _xXf.resize(_nrows, _ncols);
  for (unsigned int i = 0; i < _nrows; ++i) {
    for (unsigned int j = 0; j < _ncols; ++j) {
      _xXf(i, j) = x[j + i*_ncols];
    }
  }
  // Mean and standard deviation for each column
  VectorXf mean_vector(_ncols);
  mean_vector = _xXf.colwise().mean();
  VectorXf sd_vector(_ncols);
  unsigned int zero_sd_num = 0;
  float denom = static_cast<float>((_nrows > 1)? _nrows - 1: 1);
  for (unsigned int i = 0; i < _ncols; ++i) {
     VectorXf curr_col  = VectorXf::Constant(_nrows, mean_vector(i)); // mean(x) for column x
     curr_col = _xXf.col(i) - curr_col; // x - mean(x)
     curr_col = curr_col.array().square(); // (x-mean(x))^2  
     sd_vector(i) = sqrt((curr_col.sum())/denom);
     if (0 == sd_vector(i)) {
      zero_sd_num++;
     }
  }
  // If colums with sd == 0 are too many,
  // don't continue calculations
  if (1 > _ncols-zero_sd_num) {
    return -1;
  }
  // Delete columns where sd == 0
  MatrixXf tmp(_nrows, _ncols-zero_sd_num);
  VectorXf tmp_mean_vector(_ncols-zero_sd_num);
  unsigned int curr_col_num = 0;
  for (unsigned int i = 0; i < _ncols; ++i) {
    if (0 != sd_vector(i)) {
      tmp.col(curr_col_num) = _xXf.col(i);
      tmp_mean_vector(curr_col_num) = mean_vector(i);
      curr_col_num++;
    } else {
      _eliminated_columns.push_back(i);
    }
  }
  _ncols -= zero_sd_num;
  _xXf = tmp;
  mean_vector = tmp_mean_vector;
  tmp.resize(0, 0); tmp_mean_vector.resize(0);
  // Shift to zero
  if (true == _is_center) {
    for (unsigned int i = 0; i < _ncols; ++i) {
      _xXf.col(i) -= VectorXf::Constant(_nrows, mean_vector(i));
    }
  }
  // Scale to unit variance
  if ( (false == _is_corr) || (true == _is_scale)) {
    for (unsigned int i = 0; i < _ncols; ++i) {
     _xXf.col(i) /= sqrt(_xXf.col(i).array().square().sum()/denom);
    }
  }
  #ifdef DEBUG
    cout << "\nScaled matrix:\n";
    cout << _xXf << endl;
    cout << "\nMean before scaling:\n" << mean_vector.transpose();
    cout << "\nStandard deviation before scaling:\n" << sd_vector.transpose();
  #endif
  // When _nrows < _ncols then svd will be used.
  // If corr is true and _nrows > _ncols then will be used correlation matrix
  // (TODO): What about covariance?
  if ( (_nrows < _ncols) || (false == _is_corr)) { // Singular Value Decomposition is on
    _method = "svd";
    JacobiSVD<MatrixXf> svd(_xXf, ComputeThinV);
    VectorXf eigen_singular_values = svd.singularValues();
    VectorXf tmp_vec = eigen_singular_values.array().square();
    float tmp_sum = tmp_vec.sum();
    tmp_vec /= tmp_sum;
    // PC's standard deviation and
    // PC's proportion of variance
    _kaiser = 0;
    unsigned int lim = (_nrows < _ncols)? _nrows : _ncols;
    for (unsigned int i = 0; i < lim; ++i) {
      _sd.push_back(eigen_singular_values(i)/sqrt(denom));
      if (_sd[i] >= 1) {
        _kaiser = i + 1;
      }
      _prop_of_var.push_back(tmp_vec(i));
//.........这里部分代码省略.........