C++ ublas::vector类代码示例

void ribi::QtUblasVectorDoubleModel::SetRawData(const boost::numeric::ublas::vector<double>& data)
{
  if (!Matrix::VectorsDoubleAreEqual(m_data,data))
  {
    const int new_size = boost::numeric_cast<int>(data.size());
    const int cur_size = this->rowCount();
    if (cur_size < new_size)
    {
      this->insertRows(cur_size,new_size - cur_size,QModelIndex());
    }
    else if (cur_size > new_size)
    {
      this->removeRows(cur_size,cur_size - new_size,QModelIndex());
    }
    //Set the data before emitting signals, as the response to that signal
    //will be dependent on that data
    m_data = data;

    assert(this->rowCount() == boost::numeric_cast<int>(data.size())
      && "So emit layoutChange can work on the newest layout");

    //If you forget this line, the view displays a different number of rows than m_data has
    emit layoutChanged();

    //emit dataChanged(QModelIndex(),QModelIndex());
    const QModelIndex top_left = this->index(0,0);
    const QModelIndex bottom_right = this->index(m_data.size() - 1, 0);
    emit dataChanged(top_left,bottom_right);
  }

  assert(this->rowCount() == boost::numeric_cast<int>(this->m_data.size()));
  assert(this->rowCount() == boost::numeric_cast<int>(m_header_vertical_text.size()));
  assert(this->columnCount() == (this->rowCount() == 0 ? 0 : 1));
}

boost::python::tuple
stft_glue(const boost::numeric::ublas::vector<double> &x,
        const int length,
        const int shift)
{
    const int n_frame = ((int)x.size() - length) / shift + 1;
    const int n_freq = length / 2 + 1;
//    std::cerr << n_frame << " " << n_freq << " " << length << " " << shift << std::endl;

    std::vector<double> y(x.size() + length * 10);
    std::cerr << x.size() << " " << y.size() << std::endl;
    for(size_t i = 0; i != y.size(); ++i)
        y[i] = i < x.size() ? x[i] : 0.0;
//    std::cerr << y.size() << std::endl;
    STFT stft(length, shift, n_frame);
    //allocation
    twoDimArray<double> abs_spec(n_frame, n_freq);
    twoDimArray<std::complex<double> > phase_spec(n_frame, n_freq);

    // exec stft
    stft.exec(&(y[0]), &abs_spec, &phase_spec);

    boost::numeric::ublas::matrix<double> abs_spec_boost(n_frame, n_freq);
    boost::numeric::ublas::matrix<double> phase_spec_boost(n_frame, n_freq);

    for(int i = 0; i != n_frame; ++i){
        for(int j = 0; j != n_freq; ++j){
            abs_spec_boost(i,j) = abs_spec[i][j];
            phase_spec_boost(i,j) = std::arg(phase_spec[i][j]);
        }
    }

    auto ret = boost::python::make_tuple(abs_spec_boost, phase_spec_boost);
    return ret;
}

  int operator()(V& v) const {
     using namespace boost::numeric::bindings::blas ;

     typedef typename V::value_type                                                        value_type ;
     typedef typename bindings::remove_imaginary<value_type>::type real_type ;

     // Copy vector from reference
     for (typename V::size_type i=0; i<v_ref_.size(); ++i)
        v[i] = v_ref_(i);

     // Test blas routines and compare with reference
     real_type nrm = nrm2( v );
     if ( std::abs(nrm - norm_2(v_ref_)) > std::numeric_limits< real_type >::epsilon() * norm_2(v_ref_)) {
       std::cout << "nrm2 : " << std::abs(nrm - norm_2(v_ref_)) << " > " << std::numeric_limits< real_type >::epsilon() * norm_2(v_ref_) << std::endl ;
       return 255 ;
     }

     nrm = asum( v );
     if ( std::abs(nrm - abs_sum(v_ref_)) > std::numeric_limits< real_type >::epsilon() * abs_sum(v_ref_)) {
       std::cout << "asum : " << std::abs(nrm - abs_sum(v_ref_)) << " > " << std::numeric_limits< real_type >::epsilon() * abs_sum(v_ref_) << std::endl ;
       return 255 ;
     }

     scal( value_type(2.0), v );
     for (typename V::size_type i=0; i<v_ref_.size(); ++i)
        if (std::abs( v[i] - real_type(2.0)*v_ref_(i) ) > real_type(2.0)*std::abs(v_ref_(i))) return 255 ;

     return 0;
  }

void
UZRectMollerup::output (Log& log) const
{
  output_lazy (std::vector<double> (Theta_error.begin (), Theta_error.end ()),
               "Theta_error", log);
  output_lazy (std::vector<double> (Kedge.begin (), Kedge.end ()),
               "Kedge", log);
}

      void pprint(const boost::numeric::ublas::vector<double>& v) {
	cout << "[";
	for( uint i=0; i < v.size(); i++) {
	  cout << v(i);
	  if(i != v.size()-1)
	    cout << " , ";
	}
	cout << "]\n";
      }

double compute_distance(const ublas::vector<double>& pt1, const ublas::vector<double>& pt2)
{
	assert(pt1.size()==pt2.size());
	double sum = 0;
	for(UINT i = 0; i<pt1.size(); i++)
	{
		sum += pow((double)(pt1(i) - pt2(i)), 2.0);
	}
	return sqrt(sum);
}

double norm(const boost::numeric::ublas::vector<double>& v1, const boost::numeric::ublas::vector<double>& v2) {
    assert (v1.size() == v2.size());
    double sum = 0;
    for (size_t i=0; i<v1.size(); ++i)
    {
        double minus = v1(i) - v2(i);
        double r = minus * minus;
        sum += r;
    }

    return sqrt(sum);
}

T hotelling_t2_2test(
  boost::numeric::ublas::vector<T> mean1,
  boost::numeric::ublas::vector<T> mean2,
  boost::numeric::ublas::matrix<T> cov1,
  boost::numeric::ublas::matrix<T> cov2, 
  unsigned n1,
  unsigned n2){
  unsigned k=mean1.size();  
//int n = n1 + n2 -1;
  boost::numeric::ublas::vector<T > mean_diff=mean1-mean2;
  boost::numeric::ublas::matrix<T>  
  pooled_cov=(cov1*(n1-1)+cov2*(n2-1))*1.0/(n1+n2-2) ;
  pooled_cov *= (1.0/n1+1.0/n2 );
  boost::numeric::ublas::matrix<T>  pooled_cov_inv(pooled_cov);
  invert_matrix(   pooled_cov  ,pooled_cov_inv);
//std::cout << mean_diff << std::endl;  
  T t2_score=
  boost::numeric::ublas::inner_prod(
			      boost::numeric::ublas::prod( pooled_cov_inv,mean_diff)
			      ,mean_diff);
  T f_score= (t2_score*(n1 + n2 - 1- k))/(k*(  n1 + n2 - 2));

  std::cout << t2_score   << std::endl;
  std::cout << f_score   << std::endl;

  boost::math::fisher_f dist(k, n1+n2-1-k);
  T p_score =boost::math::cdf(dist, f_score);

  std::cout << p_score << std::endl;
  T p_comp =boost::math::cdf(complement(dist, f_score));
  std::cout << p_comp << std::endl;
  //return p_comp;
  return p_score;
}

int computeNoInliersFundamental(mat &data, mat &f, ublas::vector<bool> &inlier_mask, double &dist_std)
{
  double T_DIST = 0.5;//dist>0.1 pixels = outlier might want to make it sqrt(5.99)*sigma
  int num_inlier = 0;
  double m_dist=0;
  dist_std=0;
  vec dist(data.size1());
  inlier_mask.resize(data.size1());
  for (int i=0;i<data.size1();i++)
    inlier_mask(i)=false;
  for (int i=0;i<data.size1();i++)
    {
      vec x1(3),x2(3);
      x1(0)=data(i,0);x1(1)=data(i,1);x1(2)=1;
      x2(0)=data(i,2);x2(1)=data(i,3);x2(2)=1;
      dist(i) = fabs(inner_prod(prod (x2, f),x1));//x2 F x1
      if (dist(i)<T_DIST)
	{
	  num_inlier++;
	  inlier_mask(i)=true;
	  m_dist+=dist(i);
	}
    }
  m_dist = m_dist/ num_inlier;
  for (int i=0;i<data.size1();i++)
      if (inlier_mask(i))
	dist_std += pow(dist(i)-m_dist,2);
  dist_std = dist_std/(num_inlier-1);
  return num_inlier;
}

void generateCircularSpreading(boost::numeric::ublas::vector<double> &spr,
  const std::vector<double> &pars) {
  spr(0) = initiateCircularSpreading(pars);
  for (unsigned int i=1; i<spr.size(); ++i) {
    spr(i) = circularMap(spr(i-1), pars);
  }
}

T hotelling_t2_1test(
  boost::numeric::ublas::vector<T> mean1,
  boost::numeric::ublas::vector<T> mean2,
  boost::numeric::ublas::matrix<T> cov1,
  unsigned n1
){
  unsigned k=mean1.size();  

  boost::numeric::ublas::vector<T > mean_diff=mean1-mean2;
  boost::numeric::ublas::matrix<T>  cov_inv(cov1);
  invert_matrix(cov1  ,cov_inv);
  T t2_score=
  boost::numeric::ublas::inner_prod(
			      boost::numeric::ublas::prod( cov_inv,mean_diff)
			      ,mean_diff)*n1;
  T f_score= (t2_score*(n1 - k))/(k*(n1-1));

  std::cout << t2_score   << std::endl;
  std::cout << f_score   << std::endl;

  boost::math::fisher_f dist(k, n1-k);
  T p_score =boost::math::cdf(dist, f_score);

  std::cout << p_score << std::endl;
  T p_comp =boost::math::cdf(complement(dist, f_score));
  std::cout << p_comp << std::endl;
  //return p_comp;
  return p_score;
}

void clear(boost::numeric::ublas::vector<T> &x) {
    const ptrdiff_t n = x.size();

#pragma omp parallel for schedule(dynamic, 1024)
    for(ptrdiff_t i = 0; i < n; ++i)
        x[i] = 0;
}

void vector_print(ublas::vector<ScalarType>& v )
{

    std::cout << "vector_print!\n";
    for (unsigned int i = 0; i < v.size(); i++)
      std::cout << v(i) << "\t";
    std::cout << "\n";
}

void generateValleySpreading(boost::numeric::ublas::vector<double> &spr,
  const std::vector<double> &pars) {
  spr(0) =
    boost::random::uniform_real_distribution<>(pars[0], pars[2])(rng);
  for (unsigned int i=1; i<spr.size(); ++i) {
    spr(i) = valleyMap(spr(i-1), pars);
  }
}

bool
UZRectMollerup::converges (const ublas::vector<double>& previous,
			   const ublas::vector<double>& current) const
{
  size_t size = previous.size ();
  daisy_assert (current.size () == size);

  for (unsigned int i = 0; i < size; i++)
    {
      if (   fabs (current[i] - previous[i]) > max_absolute_difference
	  && (   iszero (previous[i])
              || iszero (current[i])
	      || (  fabs ((current[i] - previous[i]) / previous[i])
		  > max_relative_difference)))
	return false;
    }
  return true;
}