C++ AVEC::push_back方法代码示例

AVEC createAVEC(AVAR x1, AVAR x2, AVAR x3) {
  AVEC v;
  v.push_back(x1);
  v.push_back(x2);
  v.push_back(x3);
  return v;
}

TEST(AgradRevMatrix, trace_gen_quad_form_mat_grad_dvv) {
  using stan::math::trace_gen_quad_form;
  using stan::math::sum;
  using stan::agrad::matrix_v;
  using stan::math::matrix_d;
  
  matrix_v av(4,4);
  matrix_d ad(4,4);
  matrix_d bd(4,2);
  matrix_v bv(4,2);
  matrix_d cd(2,2);
  matrix_v cv(2,2);
  AVAR res;
  AVEC vars;
  VEC grad;
  size_t i,j,pos;
  
  
  bd << 100, 10,
  0,  1,
  -3, -3,
  5,  2;
  bv << 100, 10,
  0,  1,
  -3, -3,
  5,  2;
  ad << 2.0,  3.0, 4.0,   5.0, 
  6.0, 10.0, 2.0,   2.0,
  7.0,  2.0, 7.0,   1.0,
  8.0,  2.0, 1.0, 112.0;
  av << 2.0,  3.0, 4.0,   5.0, 
  6.0, 10.0, 2.0,   2.0,
  7.0,  2.0, 7.0,   1.0,
  8.0,  2.0, 1.0, 112.0;
  cd.setIdentity(2,2);
  cv.setIdentity(2,2);
  
  matrix_d dqda(bd*cd.transpose()*bd.transpose());
  matrix_d dqdb(ad*bd*cd.transpose() + ad.transpose()*bd*cd);
  
  res = trace_gen_quad_form(cd,av,bv);
  
  vars.clear();
  for (i = 0; i < 4; i++)
    for (j = 0; j < 2; j++)
      vars.push_back(bv(i,j));
  for (i = 0; i < 4; i++)
    for (j = 0; j < 4; j++)
      vars.push_back(av(i,j));
  grad = cgradvec(res,vars);
  pos = 0;
  for (i = 0; i < 4; i++)
    for (j = 0; j < 2; j++, pos++)
      EXPECT_FLOAT_EQ(grad[pos], dqdb(i,j));
  for (i = 0; i < 4; i++)
    for (j = 0; j < 4; j++, pos++)
      EXPECT_FLOAT_EQ(grad[pos], dqda(i,j));
}

TEST(AgradRevMatrix, trace_gen_inv_quad_form_ldlt_grad_vvd) {
  using stan::math::sum;
  using stan::agrad::matrix_v;
  using stan::math::matrix_d;
  
  matrix_v av(4,4);
  matrix_d ad(4,4);
  matrix_d bd(4,2);
  matrix_d cd(2,2);
  matrix_v cv(2,2);
  AVAR res;
  AVEC vars;
  VEC grad;
  size_t i,j,pos;
  
  
  bd << 100, 10,
  0,  1,
  -3, -3,
  5,  2;
  ad << 9.0,  3.0, 3.0,   3.0, 
        3.0, 10.0, 2.0,   2.0,
        3.0,  2.0, 7.0,   1.0,
        3.0,  2.0, 1.0, 112.0;
  av << 9.0,  3.0, 3.0,   3.0, 
        3.0, 10.0, 2.0,   2.0,
        3.0,  2.0, 7.0,   1.0,
        3.0,  2.0, 1.0, 112.0;
  cd.setIdentity(2,2);
  cv.setIdentity(2,2);

  stan::math::LDLT_factor<stan::agrad::var,-1,-1> ldlt_av;
  ldlt_av.compute(av);
  ASSERT_TRUE(ldlt_av.success());

  matrix_d ainv(ad.inverse());
  matrix_d dqda(-ainv*bd*cd.transpose()*bd.transpose()*ainv);
  matrix_d dqdc(bd.transpose()*ainv.transpose()*bd);
  
  // var-var
  res = trace_gen_inv_quad_form_ldlt(cv,ldlt_av,bd);
  
  vars.clear();
  for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
      vars.push_back(cv(i,j));
  for (i = 0; i < 4; i++)
    for (j = 0; j < 4; j++)
      vars.push_back(av(i,j));
  grad = cgradvec(res,vars);
  pos = 0;
  for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++, pos++)
      EXPECT_FLOAT_EQ(grad[pos], dqdc(i,j));
  for (i = 0; i < 4; i++)
    for (j = 0; j < 4; j++, pos++)
      EXPECT_FLOAT_EQ(grad[pos], dqda(i,j));
}

TEST(AgradRevMatrix, quad_form_sym_mat_grad_vv) {
  using stan::math::quad_form_sym;
  using stan::math::sum;
  using stan::agrad::matrix_v;
  using stan::math::matrix_d;
  
  matrix_v av(4,4);
  matrix_d ad(4,4);
  matrix_d bd(4,2);
  matrix_v bv(4,2);
  AVAR res;
  AVEC vars;
  VEC grad;
  size_t i,j,pos;
  
  
  bd << 100, 10,
  0,  1,
  -3, -3,
  5,  2;
  bv << 100, 10,
  0,  1,
  -3, -3,
  5,  2;
  ad << 2.0,  3.0, 4.0,   5.0, 
  3.0, 10.0, 2.0,   2.0,
  4.0,  2.0, 7.0,   1.0,
  5.0,  2.0, 1.0, 112.0;
  av << 2.0,  3.0, 4.0,   5.0, 
  3.0, 10.0, 2.0,   2.0,
  4.0,  2.0, 7.0,   1.0,
  5.0,  2.0, 1.0, 112.0;
  
  matrix_d dqda(bd*matrix_d::Ones(2,2)*bd.transpose());
  matrix_d dqdb((ad*bd + ad.transpose()*bd)*matrix_d::Ones(2,2));
  
  // var-var
  res = sum(quad_form_sym(av,bv));
  
  vars.clear();
  for (i = 0; i < 4; i++)
    for (j = 0; j < 2; j++)
      vars.push_back(bv(i,j));
  for (i = 0; i < 4; i++)
    for (j = 0; j < 4; j++)
      vars.push_back(av(i,j));
  grad = cgradvec(res,vars);
  pos = 0;
  for (i = 0; i < 4; i++)
    for (j = 0; j < 2; j++, pos++)
      EXPECT_FLOAT_EQ(grad[pos], dqdb(i,j));
  for (i = 0; i < 4; i++)
    for (j = 0; j < 4; j++, pos++)
      EXPECT_FLOAT_EQ(grad[pos], dqda(i,j));
}

TEST(AgradRevMatrix, append_col_matrix) {
  matrix_v a(2,2);
  matrix_v a_square(2,2);
  MatrixXd b(2,2);

  a << 2.0, 3.0,
       9.0, -1.0;

  b << 4.0, 3.0,
       0.0, 1.0;

  AVEC x;
  for (int i = 0; i < 2; ++i) {
    for (int j = 0; j < 2; ++j) {
      x.push_back(a(i,j));
      a_square(i, j) = square(a(i, j));
    }
  }

  AVAR append_col_ab = sum(append_col(a_square, b));

  VEC g = cgradvec(append_col_ab, x);

  size_t idx = 0;
  for (int i = 0; i < 2; i++)
    for (int j = 0; j < 2; j++)
      EXPECT_FLOAT_EQ(a(i, j).val()*2.0, g[idx++]);
  stan::math::recover_memory();
}

TEST(AgradRevMatrix,inverse_spd_inverse_spd_sum) {
  using stan::math::sum;
  using stan::math::inverse_spd;
  using stan::agrad::matrix_v;
  
  matrix_v a(4,4);
  a << 1.0, 0.0, 0.0, 0.0, 
  0.0, 1.0, 0.0, 0.0,
  0.0, 0.0, 1.0, 0.0,
  0.0, 0.0, 0.0, 1.0;
  
  AVEC x;
  for (int i = 0; i < 4; ++i)
    for (int j = 0; j < 4; ++j)
      x.push_back(a(i,j));
  
  AVAR a_inv_inv_sum = sum(inverse_spd(inverse_spd(a)));
  
  VEC g;
  a_inv_inv_sum.grad(x,g);

  size_t k = 0;
  for (int i = 0; i < 4; ++i) {
    for (int j = 0; j < 4; ++j) {
      EXPECT_FLOAT_EQ(1.0,g[k]);
      k++;
    }
  }
}

TEST(AgradRevMatrix, distance_vv) {
  using stan::math::vector_v;

  vector_v a(3), b(3);
  AVAR c;
  for (int i = -1; i < 2; i++) { // a = (-1, 0, 1), b = (1, 2, 3)
    a(i+1) = i;
    b(i+1) = i + 2;
  }
  c = stan::math::distance(a, b);
  EXPECT_FLOAT_EQ(3.464102, c.val());
  AVEC ab;
  VEC grad;
  for (size_t i = 0; i < 3; i++) {
    ab.push_back(a[i]);
    ab.push_back(b[i]);
  }
  c.grad(ab, grad);
  EXPECT_FLOAT_EQ((a(0).val() - b(0).val()) / 3.464102, grad[0]);
  EXPECT_FLOAT_EQ(-(a(0).val() - b(0).val()) / 3.464102, grad[1]);
  EXPECT_FLOAT_EQ((a(1).val() - b(1).val()) / 3.464102, grad[2]);
  EXPECT_FLOAT_EQ(-(a(1).val() - b(1).val()) / 3.464102, grad[3]);
  EXPECT_FLOAT_EQ((a(2).val() - b(2).val()) / 3.464102, grad[4]);
  EXPECT_FLOAT_EQ(-(a(2).val() - b(2).val()) / 3.464102, grad[5]);
}

TEST(AgradRevMatrix,eigenval_sum) {
  using stan::math::sum;
  using stan::agrad::matrix_v;
  using stan::agrad::vector_v;
  using stan::math::eigenvalues_sym;

  matrix_v a(3,3);
  a << 
    1.0, 2.0, 3.0,
    2.0, 5.0, 7.9,
    3.0, 7.9, 1.08;
  AVEC x = createAVEC(a(0,0), a(1,1), a(2,2), a(1,2));
  x.push_back(a(0,1));
  x.push_back(a(2,0));

  // grad sum eig = I
  vector_v a_eigenvalues = eigenvalues_sym(a);
  
  AVAR sum_a_eigenvalues = sum(a_eigenvalues);
  
  VEC g = cgradvec(sum_a_eigenvalues,x);

  EXPECT_NEAR(1.0,g[0],1.0E-11);
  EXPECT_NEAR(1.0,g[1],1.0E-11);
  EXPECT_NEAR(1.0,g[2],1.0E-11);

  EXPECT_NEAR(0.0,g[3],1.0E-10);
  EXPECT_NEAR(0.0,g[4],1.0E-10);
  EXPECT_NEAR(0.0,g[5],1.0E-10);
}

TEST(AgradRevMatrix, dot_product_vd) {
  AVEC a;
  VEC b;
  AVAR c;
  for (int i = -1; i < 2; i++) { // a = (-1, 0, 1), b = (1, 2, 3)
    a.push_back(i);
    b.push_back(i + 2);
  }
  c = dot_product(&a[0], &b[0], 3);
  EXPECT_EQ(2, c);
  VEC grad;
  c.grad(a, grad);
  EXPECT_EQ(grad[0], 1);
  EXPECT_EQ(grad[1], 2);
  EXPECT_EQ(grad[2], 3);
}

TEST(AgradRevMatrix, trace_quad_form_ldlt_mat_grad_dv) {
  using stan::math::sum;
  using stan::math::matrix_v;
  using stan::math::matrix_d;
  
  matrix_d ad(4,4);
  matrix_d bd(4,2);
  matrix_v bv(4,2);
  AVAR res;
  AVEC vars;
  VEC grad;
  size_t i,j,pos;
  
  
  bd << 100, 10,
  0,  1,
  -3, -3,
  5,  2;
  bv << 100, 10,
  0,  1,
  -3, -3,
  5,  2;
  ad << 9.0,  3.0, 3.0,   3.0, 
        3.0, 10.0, 2.0,   2.0,
        3.0,  2.0, 7.0,   1.0,
        3.0,  2.0, 1.0, 112.0;
  
  stan::math::LDLT_factor<double,-1,-1> ldlt_ad;
  ldlt_ad.compute(ad);
  ASSERT_TRUE(ldlt_ad.success());
  
  matrix_d ainv(ad.inverse());
  matrix_d dqdb(ainv*bd + ainv.transpose()*bd);
  
  // var-var
  res = trace_inv_quad_form_ldlt(ldlt_ad,bv);
  
  vars.clear();
  for (i = 0; i < 4; i++)
    for (j = 0; j < 2; j++)
      vars.push_back(bv(i,j));
  grad = cgradvec(res,vars);
  pos = 0;
  for (i = 0; i < 4; i++)
    for (j = 0; j < 2; j++, pos++)
      EXPECT_FLOAT_EQ(grad[pos], dqdb(i,j));
}

void test_sort_indices_desc(VEC val) {
  
  AVEC x;
  for(size_t i=0U; i<val.size(); i++)
    x.push_back(AVAR(val[i]));
  
  std::vector<int> val_sorted = sort_indices_desc(val);
  std::vector<int> x_sorted = sort_indices_desc(x);
  
  for(size_t i=0U; i<val.size(); i++)
    EXPECT_EQ(val_sorted[i],x_sorted[i]);

  for(size_t i=0U; i<val.size(); i++)
    for(size_t j=0U; j<val.size(); j++)
      if(val_sorted[i] == val[j])
        EXPECT_EQ(x_sorted[i],x[j]);
      else
        EXPECT_FALSE(x_sorted[i]==x[j]);
}

void test_sort_asc(VEC val) {
  using stan::math::sort_asc;
  using stan::agrad::sort_asc;
  
  AVEC x;
  for(size_t i=0U; i<val.size(); i++)
    x.push_back(AVAR(val[i]));
  
  VEC val_sorted = sort_asc(val);
  AVEC x_sorted = sort_asc(x);
  
  for(size_t i=0U; i<val.size(); i++)
    EXPECT_EQ(val_sorted[i],x_sorted[i].val());

  for(size_t i=0U; i<val.size(); i++)
    for(size_t j=0U; j<val.size(); j++)
      if(val_sorted[i] == val[j])
        EXPECT_EQ(x_sorted[i],x[j]);
      else
        EXPECT_FALSE(x_sorted[i]==x[j]);
}

TEST(AgradRevMatrix, trace_gen_inv_quad_form_ldlt_grad_vvv_basic) {
  using stan::math::sum;
  using stan::agrad::matrix_v;
  using stan::math::matrix_d;
  using stan::math::inverse;
  using stan::math::multiply;
  using stan::math::trace;
  using stan::math::transpose;
  
  matrix_v cv(2,2);
  matrix_v av(4,4);
  matrix_v bv(4,2);
  AVAR result, result_basic;
  double result_val, result_basic_val;
  AVEC vars;
  VEC grad, grad_basic;
  size_t i,j;
  stan::math::LDLT_factor<stan::agrad::var,-1,-1> ldlt_av;
  
  // calculate gradient using trace_gen_inv_quad_form_ldlt  
  bv << 100, 10,
          0,  1,
         -3, -3,
          5,  2;
  av << 9.0,  3.0, 3.0,   3.0, 
        3.0, 10.0, 2.0,   2.0,
        3.0,  2.0, 7.0,   1.0,
        3.0,  2.0, 1.0, 112.0;
  cv << 1, 2, 3, 4;

  ldlt_av.compute(av);
  ASSERT_TRUE(ldlt_av.success());
  result = trace_gen_inv_quad_form_ldlt(cv,ldlt_av,bv);
  
  vars.clear();
  for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
      vars.push_back(cv(i,j));
  for (i = 0; i < 4; i++)
    for (j = 0; j < 2; j++)
      vars.push_back(bv(i,j));
  for (i = 0; i < 4; i++)
    for (j = 0; j < 4; j++)
      vars.push_back(av(i,j));
  grad = cgradvec(result,vars);
  result_val = result.val();

  // calculate gradient using basic math
  bv << 100, 10,
          0,  1,
         -3, -3,
          5,  2;
  av << 9.0,  3.0, 3.0,   3.0, 
        3.0, 10.0, 2.0,   2.0,
        3.0,  2.0, 7.0,   1.0,
        3.0,  2.0, 1.0, 112.0;
  cv << 1, 2, 3, 4;

  matrix_v tmp = bv.transpose() * inverse(av) * bv;
  matrix_v gen_inv_quad_form = multiply(cv, tmp);
  result_basic = trace(gen_inv_quad_form);


  vars.clear();
  for (i = 0; i < 2; i++)
    for (j = 0; j < 2; j++)
      vars.push_back(cv(i,j));
  for (i = 0; i < 4; i++)
    for (j = 0; j < 2; j++)
      vars.push_back(bv(i,j));
  for (i = 0; i < 4; i++)
    for (j = 0; j < 4; j++)
      vars.push_back(av(i,j));
  grad_basic = cgradvec(result_basic,vars);
  result_basic_val = result_basic.val();
  
  // check values;
  EXPECT_FLOAT_EQ(result_basic_val, result_val);
  ASSERT_EQ(grad_basic.size(), grad.size());
  for (size_t n = 0; n < grad_basic.size(); n++) {
    EXPECT_FLOAT_EQ(grad_basic[n], grad[n]);
  }
}