C++ mat::cols方法代码示例

// m=m-m*v*v'
void sub_m_v_vT(mat &m, const vec &v)
{
  vec v2(m.rows());
  double tmp, *v2p;
  const double *vp;
  int i, j;

  it_assert(v.size() == m.cols(), "sub_m_v_vT()");

  v2p = v2._data();
  for (i = 0; i < m.rows(); i++) {
    tmp = 0.0;
    vp = v._data();
    for (j = 0; j < m.cols(); j++)
      tmp += *(vp++) * m._elem(i, j);
    *(v2p++) = tmp;
  }

  v2p = v2._data();
  for (i = 0; i < m.rows(); i++) {
    vp = v._data();
    for (j = 0; j < m.cols(); j++)
      m._elem(i, j) -= *v2p * *(vp++);
    v2p++;
  }
}

mat find_inliers(mat data, const parameters& params) {
  uint_fast32_t n_inliers = 0;
  double outlier_proportion;
  
  uvec best_indices;
  uvec indices = rand_indices(0, data.n_cols);
  
  uint_fast32_t iter = 0;
  uint_fast32_t iter_max = -1;
  while (iter < iter_max) {
    mat maybe_inliers = data.cols(indices.head(params.nfit_points));
    
    mat model = params.model_function(maybe_inliers);
    
    mat distances = params.distance_function(model, data);
    
    uvec inliers = find(distances < params.distance_threshold);
    
    if (inliers.n_elem > n_inliers) {
      n_inliers = inliers.n_elem;
      best_indices = inliers;
  
      outlier_proportion = 1.0 - (double) inliers.n_elem / (double) data.n_cols; 
      iter_max = lround(
        log(1.0 - 0.99) / log(1.0 - pow(1.0 - outlier_proportion, params.nfit_points))
      );
    }

    indices = shuffle(indices); // generate new random indexes
    ++iter;
  }

  return data.cols(best_indices);
}

mat ls_solve_od(const mat &A, const mat &B)
{
    int m=A.rows(), n=A.cols(), N=B.cols(), j;
    double beta;
    mat A2(A), B2(B), B3(n,N), submat, submat2;
    vec tmp(n), v;

//    it_assert1(m >= n, "The system is under-determined!");
    //it_assert1(m == B.rows(), "The number of rows in A must equal the number of rows in B!");

    // Perform a Householder QR factorization
    for (j=0; j<n; j++) {
	house(rvectorize(A2(j, m-1, j, j)), v, beta);
	v *= sqrt(beta);
	// 	submat.ref(A2, j,m-1, j,n-1);
 	submat = A2(j,m-1,j,n-1);
	sub_v_vT_m(submat, v);
	// 	submat.ref(B2, j,m-1, 0,N-1);
 	submat = B2(j,m-1,0,N-1);
	sub_v_vT_m(submat, v);
    }

    //    submat.ref(A2, 0,n-1,0,n-1);
    //    submat2.ref(B2, 0,n-1,0,N-1);
    submat = A2(0,n-1,0,n-1);
    submat2 = B2(0,n-1,0,N-1);
    for (j=0; j<N; j++) {
	backward_substitution(submat, submat2.get_col(j), tmp);
	B3.set_col(j, tmp);
    }
    
    return B3;
}

vec mix::process(bvec ce, mat x)
{
 vec y;	
 int N;
 bvec one = ("1");

	#if (DEBUG_LEVEL==3)
	cout << "***** mix::process *****" << endl;	
	cout << "ce=" << ce << endl;
	cout << "x=" << x << endl;
	sleep(1000);
	#endif

	N=ce.length();
	if (x.rows()!=N) { 
		throw sci_exception("mix::process - ce.size <> x.rows()", x.rows() );
	}
	if (x.cols()!=2) { 
		throw sci_exception("mix::process - x=[x1,x1] - x.cols()!=2", x.cols() );
	}

	y.set_length(N);
	for (int i=0; i<N; i++) {
		if ( bool(ce[i])) {
			y0 = G.process(one, to_vec(x(i,0)*x(i,1)))(0);
		}
		y[i]=y0;
	}
	#if (DEBUG_LEVEL==3)
	cout << "y=" << y << endl;
	cout << "+++++ mix::process +++++" << endl;	
	sleep(1000);
	#endif
	return (y);
}

    static mat qp_off_diag(const mat& Q, const mat& A) {
      mat res(Q.cols() + A.rows(), 
	      Q.cols() + A.rows());
      
      res << 
	off_diag(Q), -A.transpose(),
	A, mat::Zero(A.rows(), A.rows());
      
      return res;
    }

cmat operator*(const std::complex<double> &s, const mat &m)
{
  it_assert_debug(m.rows() > 0 && m.cols() > 0, "operator*(): Matrix of zero length");

  cmat temp(m.rows(), m.cols());

  for (int i = 0;i < m._datasize();i++) {
    temp._data()[i] = s * m._data()[i];
  }
  return temp;
}

static
void assert_mat(const mat &ref, const mat &act)
{
  ASSERT_EQ(ref.rows(), act.rows());
  ASSERT_EQ(ref.cols(), act.cols());
  for (int n = 0; n < ref.rows(); ++n) {
    for (int k = 0; k < ref.cols(); ++k) {
      ASSERT_NEAR(ref(n,k), act(n,k), tol);
    }
  }
}

    mat qp_matrix(const mat& Q, const mat& A) {
      assert(A.cols() == Q.cols());
      assert(Q.rows() == Q.cols());
      
      mat M;
      M.resize(Q.rows() + A.rows(), Q.cols() + A.rows() );

      M << Q, -A.transpose(),
	A, mat::Zero(A.rows(), A.rows() );

      return std::move(M);
    }

mat cheb(int n, const mat &x)
{
  it_assert_debug((x.rows() > 0) && (x.cols() > 0), "cheb(): empty matrix");

  mat out(x.rows(), x.cols());
  for (int i = 0; i < x.rows(); ++i) {
    for (int j = 0; j < x.cols(); ++j) {
      out(i, j) = cheb(n, x(i, j));
    }
  }
  return out;
}

void network::datasetOnLineTrain(mat dataset,int TRAINNING_TYPE,int MAX_TRAINNING_TIME) {
    int time=0;
    //上一轮训练下来的误差
    double round_error;
    mat serror;
    //赋值训练最小误差为较大的数
    min_error=DBL_MAX;
    //样本的输入和输出
    rowvec sample,out;
    while(time<MAX_TRAINNING_TIME)
    {
        round_error=0;
        for(int i=0;i<dataset.n_rows;i++)
        {
            sample=dataset.row(i).cols(0,in_vec-1);
            out=dataset.row(i).cols(in_vec,dataset.n_cols-1);
            if(TRAINNING_TYPE==BP_NONE)
            {
                simplyBPTrain(sample,out);
            }
            else if(TRAINNING_TYPE==BP_WITH_SPARSE)
            {
                withSparseTrain(sample,out);
            }
            else{
                return;
            }
        }
        updateOut(dataset.cols(0,in_vec-1));
        serror=dataset.cols(in_vec,dataset.n_cols-1)-output;
        round_error=norm(sum(serror%serror/dataset.n_rows),2);
        time++;
        if(round_error<min_error)
        {
            min_error=round_error;
            for(int i=layer_num;i>=1;i--)
            {
                mw[i]=w[i];
                mo[i]=o[i];
            }
        }
        if(time%SHOW_TIME==0)
            cout<<"第"<<time<<"次训练的误差变化为"<<setprecision(50)<<round_error<<endl;
        if(round_error<tor_error)
        {
            error=round_error;
            break;
        }
    }
    cout<<"本次训练了"<<time<<"次，最小训练误差为"<<min_error<<endl<<"最终训练误差为"<<round_error<<endl;
}

void forward_substitution(const mat &L, int p, const vec &b, vec &x)
{
    assert( L.rows() == L.cols() && L.cols() == b.size() && b.size() == x.size() && p <= L.rows()/2 );
    int n = L.rows(), i, j;

    x=b;
  
    for (j=0;j<n;j++) {
	x(j)/=L(j,j);
	for (i=j+1;i<MIN(j+p+1,n);i++) {
	    x(i)-=L(i,j)*x(j);
	}
    }
}

inline void update_WtA(mat & WtA, const mat & W, const mat & W1, const mat & H2, const mat & A)
{
	// compute WtA = (W[:, 0:k-1], W1)^T (A - W[, k:end] H2^T)
	if (H2.empty())
		update_WtA(WtA, W, W1, A);
	else
	{
		int k = W.n_cols - H2.n_cols;
		//std::cout << "1.3" << std::endl;
		//A.print("A = ");
		//(W.cols(k, W.n_cols-1) * H2.t()).print("W[, k:] = ");
		update_WtA(WtA, W.cols(0, k-1), W1, A - W.cols(k, W.n_cols-1) * H2.t());
	}
}

void backward_substitution(const mat &U, int q, const vec &b, vec &x)
{
    assert( U.rows() == U.cols() && U.cols() == b.size() && b.size() == x.size() && q <= U.rows()/2);
    int n = U.rows(), i, j;

    x=b;
  
    for (j=n-1; j>=0; j--) {
	x(j) /= U(j,j);
	for (i=MAX(0,j-q); i<j; i++) {
	    x(i)-=U(i,j)*x(j);
	}
    }
}

bool schur(const mat &A, mat &U, mat &T)
{
  it_assert_debug(A.rows() == A.cols(), "schur(): Matrix is not square");

  char jobvs = 'V';
  char sort = 'N';
  int info;
  int n = A.rows();
  int lda = n;
  int ldvs = n;
  int lwork = 3 * n; // This may be choosen better!
  int sdim = 0;
  vec wr(n);
  vec wi(n);
  vec work(lwork);

  T.set_size(lda, n, false);
  U.set_size(ldvs, n, false);

  T = A; // The routine overwrites input matrix with eigenvectors

  dgees_(&jobvs, &sort, 0, &n, T._data(), &lda, &sdim, wr._data(), wi._data(),
         U._data(), &ldvs, work._data(), &lwork, 0, &info);

  return (info == 0);
}

/**
 * Compute the covariance matrix of a set of inputs
 * @param C The covariance matrix of X
 * @param X A matrix of inputs (one input per row)
 */
void CovarianceFunction::covariance(mat& C, const mat& X) const
{
	// ensure that data dimensions match supplied covariance matrix
	assert(C.rows() == X.rows());
	assert(C.cols() == X.rows());

	if (X.rows() == 1)
	{
	    C.set(0, 0, computeDiagonalElement(X.get_row(0)));
	    return;
	}
	
	// calculate the lower and upper triangles
	double d;
	
	for(int i=0; i<X.rows() ; i++)
	{
		for(int j=0; j<i; j++)
		{
		    d = computeElement(X.get_row(i), X.get_row(j));
		    C.set(i, j, d);
			C.set(j, i, d);
		}
	}

	// calculate the diagonal part
	for(int i=0; i<X.rows() ; i++)
	{
		C.set(i, i, computeDiagonalElement(X.get_row(i)));
	}
}