C конечно-разностное дифференцирование — проектирование

#c #templates #design-patterns #traits #expression-templates

#c #шаблоны #проектирование-шаблоны #Трейты #выражение-шаблоны

Вопрос:

пусть A будет:

 class A {
  std::vector<double> values_;
public:
  A(const std::vector<double> amp;values) : values_(values){};

  void bumpAt(std::size_t i, const double amp;df) {
    values_[i]  = df;

  virtual method1();
  virtual method2();
  ...
}

class B : public A {
  overrides methods
  ...
}
  

для простоты рассмотрим функцию:

 double foo(input1, input2, ..., const A amp;a, const B amp;b, inputK, ...) {
 /* do complex stuff */
 return ...;
}
  

мы хотим дифференцироваться foo() в отношении его аргументов. Следовательно, чувствительность первого порядка d foo/d a является std::vector<double> с размером, равным a.size() . То же самое рассуждение относится к d foo/d b .

Наивная реализация будет выглядеть следующим образом:

 
// compute d foo/d a 
std::vector<double> computeDfDa(input1, input2, ..., const A amp;a, const B amp;b, inputK, ..., double da = 1.0){
  std::vector<double> dfda = {};

  auto aUp = a.copy();
  auto aDown = a.copy();
  for (auto i = 0; i < a.size();   i) {

      // bump up
      aUp.bumpAt(i, da);
      // bump down
      aDown.bumpAt(i, -da);

      auto up = foo(input1, input2, ..., aUp, b, inputK, ...);
      auto down = foo(input1, input2, ..., aDown, b, inputK, ...);

      auto derivative = (up - down) / 2.0 / da;
      dfda.pushback(derivative);

      // revert bumps
      aUp.bumpAt(i, -da);
      aDown.bumpAt(i, da);
  }
   return dfda;
}

// compute d foo/d b 
std::vector<double> computeDfDb(input1, input2, ..., const A amp;a, const B amp;b, inputK, ..., double db = 0.01){
  std::vector<double> dfdb = {};

  auto bUp = b.copy();
  auto bDown = b.copy();
  for (auto i = 0; i < a.size();   i) {

      // bump up
      bUp.bumpAt(i, db);
      // bump down
      bDown.bumpAt(i, -db);

      auto up = foo(input1, input2, ..., a, bUp, inputK, ...);
      auto down = foo(input1, input2, ..., a, bDown, inputK, ...);

      auto derivative = (up - down) / 2.0 / db;
      dfdb.pushback(derivative);

      // revert bumps
      bUp.bumpAt(i, -db);
      bDown.bumpAt(i, db);
  }
   return dfdb;
}
  

Это работает хорошо, однако у нас в основном один и тот же код для computeDfDa() и для computeDfDb() .

Существует ли какой-либо шаблон проектирования, который позволил бы иметь уникальную (возможно, шаблонную) функцию, которая автоматически понимала бы, какой ввод следует использовать?

Пожалуйста, обратите внимание, что положение a и b во входных данных не является коммутативным.

Если сложность и количество входных данных foo() намного больше, наивное решение сгенерирует много бесполезного кода, поскольку нам придется писать computeDfDx() функцию для каждого ввода x foo() .

Комментарии:

1. computeDfDa( a, b, 1.0) Имеет тот же результат, что и computeDfDa( b, a, 1.0) ? Вы не показали compute , является ли аргумент compute взаимозаменяемым?

2. @LouisGo нет, реальный случай заключается в том, что a и b являются сложными объектами, а вычисление не является коммутативным

3. С точки зрения дизайна ваша дифференцирующая логика не должна быть привязана к вашим классам. Существует функция с N аргументами, и вам нужна пара векторов с плавающей запятой длиной N, которые содержат локальную производную этой функции в некоторой точке. Если вы кодируете это независимо от ваших классов, скорее всего, сложность вашего кода снижается, а также дублирование кода.

4. @BitTickler вычисления внутри foo() сильно зависят от типов классов

5. @mastro Речь идет о простоте кода. Одним из инструментов для этого является разделение задач. Разность некоторой функции в точке — это непересекающаяся проблема, которая не должна переплетаться с тем, что еще должны делать ваши классы. Затем классы вызывают функцию дифференцирования и предоставляют аргументы по своему усмотрению. Также проще тестировать (модульные тесты). И для повторного использования. Вы упомянули, что это критически важно для безопасности. Поэтому важно иметь хорошие тесты. И последнее, что не менее важно. Не все должно быть ООП. Это отличный пример для этого.

Ответ №1:

Поскольку compute порядок тот же, но цикл итерации проходит через разные контейнеры, вы можете реорганизовать эту функцию.

 std::vector<double> computeLoop( std::vector<double> amp;v, std::vector<double> const amp;computeArg1, std::vector<double> const amp;computeArg2, double d = 1.0 )
{
  std::vector<double> dfd = {};
  for (auto i = 0; i < v.size();   i) {
      // bump up
      v[i]  = d;
      auto up = compute(computeArg1, computeArg2);
      v[i] -= d;

      // bump down
      v[i] -= d;
      auto down = compute(computeArg1, computeArg2);
      v[i]  = d;

      auto derivative = (up - down) / 2.0 / d;
      dfd.pushback(derivative);
  }
}
  

Фактический вызов.

 auto dfda = computeLoop( a, a, b );
auto dfdb = computeLoop( b, a, b );
  

Пусть v передается по ссылке, но это может вызвать проблемы с обслуживанием. Потому что v может быть та же ссылка, что и computeArg1 или computeArg2 , однако в computeLoop это не очевидно. Кто-то может неосознанно нарушить код в будущем.

Комментарии:

1. приведенное выше вычисление приведет к тому, что производная всегда будет равна 0 — вы сталкиваетесь v , но нигде не используете ее

2. @mastro, я редактирую свой ответ, чтобы он v передавался по ссылке, поэтому его также можно было изменить.

3. спасибо, однако это небезопасно, приведенный выше код используется в критически важном для безопасности приложении

Ответ №2:

ИМХО, проблема в том, что происходит сдвиг уровня абстракции.

Эти классы A, B являются либо …

  1. просто замаскированные векторы
  2. не векторы, а что-то совсем другое.

В случае (1) должна быть возможность переписать в форму, подобную этой:

 #include <... the usual suspects ... >

using VecF64 = std::vector<double>;
using VecVecF64 = std::vector<VecF64>;

// foo, dropping allusions of encapsulations. It KNOWS those A, B are vectors.
// Hence we can write it without knowledge of A, B types.
double foo(const VecVecF64amp; args) {
  return <result-of-complicated-computations>;
}

// With that, it is also easier to write the differentiation function
VecVecF64 fooGradients( const VecVecF64amp; args, double delta = 0.01) {
  VecVecF64 resu<
  result.reserve(args.size());
  // for each arg vector in args, do your thing.
  // .. 
  return resu<
}
  

В случае (2), если вы все инкапсулированы и скрываете природу A, B,
откуда вы знаете количество удвоений, которые foo может использовать для своего вычисления? Что приводит к вопросу о размере вашего вектора градиента для A и делает невозможным реализацию этой идеи bump.

Я предполагаю, что у вас проблема типа case 1.

Комментарии:

1. спасибо, это случай (2) — внутри foo другие методы из A для извлечения соответствующих данных с помощью интерполяции, т. Е. Что-то вроде auto v = a.interpolate(input1). Foo не обязательно знать размер A, но это все равно можно сделать с помощью метода A.size() (как это делается в цикле). Все это в настоящее время работает более чем нормально, поэтому я бы сказал, что это невозможно, может быть, я просто недостаточно четко выразился