310 std::vector<double> & values,
317 Assert(values.size() == this->n() || values.size() == 0,
319 Assert(grads.size() == this->n() || grads.size() == 0,
321 Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
323 Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
325 Assert(fourth_derivatives.size() == this->n() ||
326 fourth_derivatives.size() == 0,
330 update_grads = (grads.size() == this->n()),
331 update_grad_grads = (grad_grads.size() == this->n()),
333 update_4th_derivatives = (fourth_derivatives.size() == this->n());
336 unsigned int n_values_and_derivatives = 0;
338 n_values_and_derivatives = 1;
340 n_values_and_derivatives = 2;
341 if (update_grad_grads)
342 n_values_and_derivatives = 3;
344 n_values_and_derivatives = 4;
345 if (update_4th_derivatives)
346 n_values_and_derivatives = 5;
353 const unsigned int n_polynomials = polynomials.size();
354 boost::container::small_vector<ndarray<double, dim, 5>, 20> values_1d(
356 if (n_values_and_derivatives == 1)
357 for (
unsigned int i = 0; i < n_polynomials; ++i)
358 for (
unsigned int d = 0; d < dim; ++d)
359 values_1d[i][d][0] = polynomials[i].value(p(d));
361 for (
unsigned int i = 0; i < n_polynomials; ++i)
362 for (
unsigned d = 0; d < dim; ++d)
363 polynomials[i].value(p(d),
364 n_values_and_derivatives,
365 values_1d[i][d].data());
367 unsigned int indices[3];
368 unsigned int ind = 0;
369 for (indices[2] = 0; indices[2] < (dim > 2 ? n_polynomials : 1); ++indices[2])
370 for (indices[1] = 0; indices[1] < (dim > 1 ? n_polynomials : 1);
372 if (n_values_and_derivatives == 1)
373 for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
375 double value = values_1d[indices[0]][0][0];
376 for (
unsigned int d = 1; d < dim; ++d)
377 value *= values_1d[indices[d]][d][0];
378 values[index_map_inverse[ind]] = value;
381 for (indices[0] = 0; indices[0] < n_polynomials; ++indices[0], ++ind)
383 const unsigned int i = index_map_inverse[ind];
387 double value = values_1d[indices[0]][0][0];
388 for (
unsigned int x = 1; x < dim; ++x)
389 value *= values_1d[indices[x]][x][0];
394 for (
unsigned int d = 0; d < dim; ++d)
397 for (
unsigned int x = 0; x < dim; ++x)
398 grad *= values_1d[indices[x]][x][(d == x) ? 1 : 0];
402 if (update_grad_grads)
403 for (
unsigned int d1 = 0; d1 < dim; ++d1)
404 for (
unsigned int d2 = 0; d2 < dim; ++d2)
407 for (
unsigned int x = 0; x < dim; ++x)
409 unsigned int derivative = 0;
415 der2 *= values_1d[indices[x]][x][derivative];
417 grad_grads[i][d1][d2] = der2;
421 for (
unsigned int d1 = 0; d1 < dim; ++d1)
422 for (
unsigned int d2 = 0; d2 < dim; ++d2)
423 for (
unsigned int d3 = 0; d3 < dim; ++d3)
426 for (
unsigned int x = 0; x < dim; ++x)
428 unsigned int derivative = 0;
436 der3 *= values_1d[indices[x]][x][derivative];
438 third_derivatives[i][d1][d2][d3] = der3;
441 if (update_4th_derivatives)
442 for (
unsigned int d1 = 0; d1 < dim; ++d1)
443 for (
unsigned int d2 = 0; d2 < dim; ++d2)
444 for (
unsigned int d3 = 0; d3 < dim; ++d3)
445 for (
unsigned int d4 = 0; d4 < dim; ++d4)
448 for (
unsigned int x = 0; x < dim; ++x)
450 unsigned int derivative = 0;
460 der4 *= values_1d[indices[x]][x][derivative];
462 fourth_derivatives[i][d1][d2][d3][d4] = der4;
701 std::vector<double> & values,
707 Assert(values.size() == this->n() || values.size() == 0,
709 Assert(grads.size() == this->n() || grads.size() == 0,
711 Assert(grad_grads.size() == this->n() || grad_grads.size() == 0,
713 Assert(third_derivatives.size() == this->n() || third_derivatives.size() == 0,
715 Assert(fourth_derivatives.size() == this->n() ||
716 fourth_derivatives.size() == 0,
720 update_grads = (grads.size() == this->n()),
721 update_grad_grads = (grad_grads.size() == this->n()),
723 update_4th_derivatives = (fourth_derivatives.size() == this->n());
728 unsigned int n_values_and_derivatives = 0;
730 n_values_and_derivatives = 1;
732 n_values_and_derivatives = 2;
733 if (update_grad_grads)
734 n_values_and_derivatives = 3;
736 n_values_and_derivatives = 4;
737 if (update_4th_derivatives)
738 n_values_and_derivatives = 5;
744 std::size_t max_n_polynomials = 0;
745 for (
unsigned int d = 0; d < dim; ++d)
746 max_n_polynomials =
std::max(max_n_polynomials, polynomials[d].size());
750 for (
unsigned int d = 0; d < dim; ++d)
751 for (
unsigned int i = 0; i < polynomials[d].size(); ++i)
752 polynomials[d][i].value(p(d),
753 n_values_and_derivatives - 1,
756 for (
unsigned int i = 0; i < this->n(); ++i)
762 std::array<unsigned int, dim> indices;
763 compute_index(i, indices);
768 for (
unsigned int x = 0; x < dim; ++x)
769 values[i] *= v(x, indices[x])[0];
773 for (
unsigned int d = 0; d < dim; ++d)
776 for (
unsigned int x = 0; x < dim; ++x)
777 grads[i][d] *= v(x, indices[x])[d == x ? 1 : 0];
780 if (update_grad_grads)
781 for (
unsigned int d1 = 0; d1 < dim; ++d1)
782 for (
unsigned int d2 = 0; d2 < dim; ++d2)
784 grad_grads[i][d1][d2] = 1.;
785 for (
unsigned int x = 0; x < dim; ++x)
787 unsigned int derivative = 0;
793 grad_grads[i][d1][d2] *= v(x, indices[x])[derivative];
798 for (
unsigned int d1 = 0; d1 < dim; ++d1)
799 for (
unsigned int d2 = 0; d2 < dim; ++d2)
800 for (
unsigned int d3 = 0; d3 < dim; ++d3)
802 third_derivatives[i][d1][d2][d3] = 1.;
803 for (
unsigned int x = 0; x < dim; ++x)
805 unsigned int derivative = 0;
813 third_derivatives[i][d1][d2][d3] *=
814 v(x, indices[x])[derivative];
818 if (update_4th_derivatives)
819 for (
unsigned int d1 = 0; d1 < dim; ++d1)
820 for (
unsigned int d2 = 0; d2 < dim; ++d2)
821 for (
unsigned int d3 = 0; d3 < dim; ++d3)
822 for (
unsigned int d4 = 0; d4 < dim; ++d4)
824 fourth_derivatives[i][d1][d2][d3][d4] = 1.;
825 for (
unsigned int x = 0; x < dim; ++x)
827 unsigned int derivative = 0;
837 fourth_derivatives[i][d1][d2][d3][d4] *=
838 v(x, indices[x])[derivative];