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

void checkDenseVectorAssignment(V1& v1, V2 const& v2){
	BOOST_REQUIRE_EQUAL(v1.size(),v2.size());
	//indexed access
	for(std::size_t i = 0; i != v2.size(); ++i){
		v1(i) = 0;
		BOOST_CHECK_EQUAL(v1(i),0);
		v1(i) = v2(i);
		BOOST_CHECK_EQUAL(v1(i),v2(i));
		v1(i) = 0;
		BOOST_CHECK_EQUAL(v1(i),0);
	}
	//iterator access rows
	typedef typename V1::iterator Iter;
	BOOST_REQUIRE_EQUAL(v1.end()-v1.begin(), v1.size());
	std::size_t k = 0;
	for(Iter it = v1.begin(); it != v1.end(); ++it,++k){
		BOOST_CHECK_EQUAL(k,it.index());
		*it = 0;
		BOOST_CHECK_EQUAL(v1(k),0);
		*it = v2(k);
		BOOST_CHECK_EQUAL(v1(k),v2(k));
		*it = 0;
		BOOST_CHECK_EQUAL(v1(k),0);
	}
	//test that the actual iterated length equals the number of elements
	BOOST_CHECK_EQUAL(k, v2.size());
}

void bi::cross(const M1 X, const M2 Y, const V1 muX, const V2 muY,
    M3 SigmaXY) {
  /* pre-conditions */
  BI_ASSERT(X.size2() == muX.size());
  BI_ASSERT(Y.size2() == muY.size());
  BI_ASSERT(X.size1() == Y.size1());
  BI_ASSERT(SigmaXY.size1() == muX.size() && SigmaXY.size2() == muY.size());

  const int N = X.size1();

  gemm(1.0/(N - 1.0), X, Y, 0.0, SigmaXY, 'T', 'N');
  ger(-N/(N - 1.0), muX, muY, SigmaXY);
}

void bi::MetropolisResamplerHost::ancestors(Random& rng, const V1 lws,
    V2 as, int B) {
  const int P1 = lws.size(); // number of particles
  const int P2 = as.size(); // number of ancestors to draw

  #pragma omp parallel
  {
    real alpha, lw1, lw2;
    int k, p1, p2, p;

    #pragma omp for
    for (p = 0; p < P2; ++p) {
      p1 = p;
      lw1 = lws(p);
      for (k = 0; k < B; ++k) {
        p2 = rng.uniformInt(0, P1 - 1);
        lw2 = lws(p2);
        alpha = rng.uniform<real>();

        if (bi::log(alpha) < lw2 - lw1) {
          /* accept */
          p1 = p2;
          lw1 = lw2;
        }
      }

      /* write result */
      as(p) = p1;
    }
  }
}

void bi::mean(const UniformPdf<V1>& q, V2 mu) {
  /* pre-condition */
  BI_ASSERT(q.size() == mu.size());

  axpy(0.5, q.lower(), mu, true);
  axpy(0.5, q.upper(), mu);
}

void bi::cov(const M1 X, const V1 w, const V2 mu, M2 Sigma) {
  /* pre-conditions */
  BI_ASSERT(X.size2() == mu.size());
  BI_ASSERT(X.size1() == w.size());
  BI_ASSERT(Sigma.size1() == mu.size() && Sigma.size2() == mu.size());

  typedef typename V1::value_type T;
  typename sim_temp_matrix<M2>::type Y(X.size1(), X.size2());
  typename sim_temp_matrix<M2>::type Z(X.size1(), X.size2());
  typename sim_temp_vector<V2>::type v(w.size());

  T Wt = sum_reduce(w);
  Y = X;
  sub_rows(Y, mu);
  sqrt_elements(w, v);
  gdmm(1.0, v, Y, 0.0, Z);
  syrk(1.0/Wt, Z, 0.0, Sigma, 'U', 'T');
  // alternative weight: 1.0/(Wt - W2t/Wt)
}

void bi::var(const M1 X, const V1 w, const V2 mu, V3 sigma) {
  /* pre-conditions */
  BI_ASSERT(X.size2() == mu.size());
  BI_ASSERT(X.size1() == w.size());
  BI_ASSERT(sigma.size() == mu.size());

  typedef typename V1::value_type T1;
  typename sim_temp_matrix<M1>::type Z(X.size1(), X.size2());
  typename sim_temp_matrix<M1>::type Y(X.size1(), X.size2());
  typename sim_temp_vector<V2>::type v(w.size());

  T1 Wt = sum_reduce(w);
  Z = X;
  sub_rows(Z, mu);
  sqrt_elements(w, v);
  gdmm(1.0, v, Z, 0.0, Y);
  dot_columns(Y, sigma);
  divscal_elements(sigma, Wt, sigma);
  // alternative weight: 1.0/(Wt - W2t/Wt)
}

 template <class V1, class V2> double dot_impl(
     const V1 &v1, const V2 &v2) {
   assert(v1.size() == v2.size());
   if(v1.stride() > 0 && v2.stride() > 0){
     return ddot(v1.size(),
                 v1.data(), v1.stride(),
                 v2.data(), v2.stride());
   }else{
     double ans = 0;
     for(int i = 0; i < v1.size(); ++i){
       ans += v1[i] * v2[i];
     }
     return ans;
   }
 }

void bi::MultinomialResamplerHost::ancestors(Random& rng, const V1 lws, V2 as,
    MultinomialPrecompute<ON_HOST>& pre)
    throw (ParticleFilterDegeneratedException) {
  typedef typename V1::value_type T1;

  const int P = as.size();
  const int lwsSize = lws.size();

  T1 lW;

  /* weights */
  if (pre.W > 0) {
    lW = bi::log(pre.W);

    #pragma omp parallel
    {
      int Q = P/bi_omp_max_threads;
      int start = bi_omp_tid*Q + bi::min(bi_omp_tid, P % bi_omp_max_threads); // min() handles leftovers
      if (bi_omp_tid < P % bi_omp_max_threads) {
        ++Q; // pick up a leftover
      }

      int i, j = lwsSize;
      T1 lMax = 0.0, lu;
      for (i = Q; i > 0; --i) {
        lMax += bi::log(rng.uniform<T1>())/i;
        lu = lW + lMax;

        while (j > 0 && lu < bi::log(pre.Ws(j - 1))) {
          --j;
        }
        if (pre.sort) {
          as(start + i - 1) = pre.ps(j);
        } else {
          as(start + i - 1) = j;
        }
      }
    }
  } else {
    throw ParticleFilterDegeneratedException();
  }

  /* post-condition */
  BI_ASSERT(max_reduce(as) < lws.size());
}

void bi::hist(const V1 x, const V2 w, V3 c, V4 h) {
  /* pre-condition */
  BI_ASSERT(x.size() == w.size());
  BI_ASSERT(c.size() == h.size());
  BI_ASSERT(!V3::on_device);
  BI_ASSERT(!V4::on_device);

  typedef typename V1::value_type T1;
  typedef typename V2::value_type T2;

  const int P = x.size();
  const int B = c.size();
  T1 mx, mn;
  int i, j;
  typename temp_host_vector<T1>::type xSorted(P);
  typename temp_host_vector<T2>::type wSorted(P);
  xSorted = x;
  wSorted = w;

  bi::sort_by_key(xSorted, wSorted);
  mn = xSorted[0];
  mx = xSorted[xSorted.size() - 1];

  /* compute bin right edges */
  for (j = 0; j < B; ++j) {
    c[j] = mn + (j + 1)*(mx - mn)/B;
  }

  /* compute bin heights */
  h.clear();
  for (i = 0, j = 0; i < P; ++i) {
    if (xSorted[i] >= c[j] && j < B - 1) {
      ++j;
    }
    h[j] += wSorted[i];
  }

  /* compute bin centres */
  for (j = B - 1; j > 0; --j) {
    c[j] = 0.5*(c[j - 1] + c[j]);
  }
  c[0] = 0.5*(mn + c[0]);
}