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

 static void apply(Alpha a, const V1 &x, const V2 &y, Beta b, V3 &z)
 {
     if (!math::is_zero(b))
         z.array() = a * x.array() * y.array() + b * z.array();
     else
         z.array() = a * x.array() * y.array();
 }

inline void bi::exp_rows(M2 X, const V2& is) {
  BOOST_AUTO(iter, is.begin());
  BOOST_AUTO(end, is.end());
  for (; iter != end; ++iter) {
    BOOST_AUTO(row1, row(X, *iter));
    exp_elements(row1, row1);
  }
}

inline void bi::exp_columns(M2 X, const V2& is) {
  BOOST_AUTO(iter, is.begin());
  BOOST_AUTO(end, is.end());
  for (; iter != end; ++iter) {
    BOOST_AUTO(col, column(X, *iter));
    exp_elements(col, col);
  }
}

bool V2::operator == (V2 const& rhs) const{
	if (x() == rhs.x() and y() == rhs.y()) {
		return true;
	}
	else {
		return false;
	}
}

void bi::det_rows(const M2 X, const V2& is, V3 det) {
  BOOST_AUTO(iter, is.begin());
  BOOST_AUTO(end, is.end());

  set_elements(det, 1.0);
  for (; iter != end; ++iter) {
    mul_elements(det, column(X, *iter), det);
  }
}

    V5() {
	thunk_OK = this;
	if( bar() < 1000 ) {
	    V2 *p = this;
	    p->foo( bar() );
	    V3 *q = this;
	    q->foo( bar() );
	}
    }

void
MiscCheckTheParallelEnvironment(const V1& vec1, const V2& vec2)
{
  const BaseEnvironment& env = vec1.env();

  if (env.numSubEnvironments() == (unsigned int) env.fullComm().NumProc()) {
    UQ_FATAL_TEST_MACRO(env.subRank() != 0,
                        env.worldRank(),
                        "MiscCheckTheParallelEnvironment<V1,V2>()",
                        "there should exist only one processor per sub environment");
    UQ_FATAL_TEST_MACRO((vec1.numOfProcsForStorage() != 1) ||
                        (vec2.numOfProcsForStorage() != 1),
                        env.worldRank(),
                        "MiscCheckTheParallelEnvironment<V1,V2>()",
                        "only 1 processor (per sub environment) should be necessary for the storage of a parameter vector");
  }
  else if (env.numSubEnvironments() < (unsigned int) env.fullComm().NumProc()) {
    UQ_FATAL_TEST_MACRO(env.fullComm().NumProc()%env.numSubEnvironments() != 0,
                        env.worldRank(),
                        "MiscCheckTheParallelEnvironment<V1,V2>()",
                        "total number of processors should be a multiple of the number of sub environments");
    unsigned int numProcsPerSubEnvironment = env.fullComm().NumProc()/env.numSubEnvironments();
    UQ_FATAL_TEST_MACRO(env.subComm().NumProc() != (int) numProcsPerSubEnvironment,
                        env.worldRank(),
                        "MiscCheckTheParallelEnvironment<V1,V2>()",
                        "inconsistent number of processors per sub environment");
    if ((vec1.numOfProcsForStorage() == 1) &&
        (vec2.numOfProcsForStorage() == 1)) {
      // Ok
    }
    else if ((vec1.numOfProcsForStorage() == numProcsPerSubEnvironment) &&
             (vec2.numOfProcsForStorage() == numProcsPerSubEnvironment)) {
      UQ_FATAL_TEST_MACRO(true,
                          env.worldRank(),
                          "MiscCheckTheParallelEnvironment<V1,V2>()",
                          "parallel vectors are not supported yet");
    }
    else {
      UQ_FATAL_TEST_MACRO(true,
                          env.worldRank(),
                          "MiscCheckTheParallelEnvironment<V1,V2>()",
                          "number of processors required for a vector storage should be equal to either 1 or to the number of processors in the sub environment");
    }
  }
  else {
    UQ_FATAL_TEST_MACRO(true,
                        env.worldRank(),
                        "MiscCheckTheParallelEnvironment<V1,V2>()",
                        "number of processors per sub environment is less than 1!");
  }

  return;
}

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

/**
 * Draw model.
 * Care about model shift during move.
 */
    void
View::drawModel(Cube *model)
{
    if (!model->isLost()) {
        V2 screenPos = getScreenPos(model);

        Anim::eSide side = Anim::SIDE_LEFT;
        if (!model->isLeft()) {
            side = Anim::SIDE_RIGHT;
        }
        model->anim()->drawAt(m_screen,
                screenPos.getX(), screenPos.getY(), side);
    }
}

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

inline void frenet(const V2& d1, V2& t, V2& n){
	t = d1;
	t.normalize();
	// normal according to right-hand rule
	n[0] =-t[1];
	n[1] = t[0];
}

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 rot (const T1 &t1, V1 &v1, const T2 &t2, V2 &v2) 
 {
         typedef typename promote_traits<typename V1::value_type, typename V2::value_type>::promote_type promote_type;
         vector<promote_type> vt (t1 * v1 + t2 * v2);
         v2.assign (- t2 * v1 + t1 * v2);
         v1.assign (vt);
     }

real bi::det_vector(const V2 x, const V3& is) {
  BOOST_AUTO(iter, is.begin());
  BOOST_AUTO(end, is.end());
  real det = 1.0;
  for (; iter != end; ++iter) {
    det *= *(x.begin() + *iter);
  }
  return det;
}

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

	V2 ControllerInput::GetMoveVector()
    {
        // direction vector
        ::Ogre::Vector3 move = ::Ogre::Vector3::ZERO;

        // get key states
        bool isForward		= OISKeyboard->isKeyDown(ActiveKeyBinding->MoveForward);
        bool isBackwards	= OISKeyboard->isKeyDown(ActiveKeyBinding->MoveBackward);
        bool isLeft			= OISKeyboard->isKeyDown(ActiveKeyBinding->MoveLeft);
        bool isRight		= OISKeyboard->isKeyDown(ActiveKeyBinding->MoveRight);
		bool uiReadsKeys	= ControllerUI::ProcessingInput;

        // orientation axes of controlled node
		Quaternion orient = Avatar->SceneNode->getOrientation();

        ::Ogre::Vector3 xAxis = orient.xAxis().normalisedCopy();
        ::Ogre::Vector3 zAxis = orient.zAxis().normalisedCopy();

        // apply keystates on direction vector
        if ((isForward && !uiReadsKeys) || IsBothMouseDown || isAutoMove)            
            move += -zAxis;
            
        if (isBackwards && !uiReadsKeys)
            move += zAxis;

        if (isLeft && !uiReadsKeys)
            move += -xAxis;

        if (isRight && !uiReadsKeys)
            move += xAxis;

        // get a V2 variant on ground
        V2 direction;
		direction.X = move.x;
		direction.Y = move.z;
		
		// normalize
		direction.Normalize();
      
        return direction;
    };

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

  const int N = X.size1();
  typename sim_temp_matrix<M1>::type Z(X.size2(), X.size1());
  Z = X;
  sub_rows(Z, mu);
  dot_columns(Z, sigma);
  scal(1.0/(N - 1.0), sigma);
}

/**
 * Move direct to the destination.
 * @param field destination field
 * @return whether a move was made
 */
bool
MouseControl::moveHardTo(const V2 &field) const
{
    bool moved = false;
    const Unit *unit = m_controls->getActive();
    if (unit) {
        V2 loc = unit->getLoc();
        if (field.getX() < loc.getX()) {
            moved = m_controls->makeMove(unit->myOrder(Dir::DIR_LEFT));
        }
        else if (loc.getX() + unit->getW() <= field.getX()) {
            moved = m_controls->makeMove(unit->myOrder(Dir::DIR_RIGHT));
        }
        else if (field.getY() < loc.getY()) {
            moved = m_controls->makeMove(unit->myOrder(Dir::DIR_UP));
        }
        else if (loc.getY() + unit->getH() <= field.getY()) {
            moved = m_controls->makeMove(unit->myOrder(Dir::DIR_DOWN));
        }
    }
    return moved;
}