Comments 61
На шаблонах C++ этот метод выглядит конечно не так красиво как в оригинале в книге Structure and Interpretation of Computer Programs.
Я только не понял каким образом это помогло решить исходную задачу про функцию формулировка которой была разбросана по нескольким функциям и была достаточно сложной для записи.
Я только не понял каким образом это помогло решить исходную задачу про функцию формулировка которой была разбросана по нескольким функциям и была достаточно сложной для записи.
Я только не понял каким образом это помогло решить исходную задачу про функцию формулировка которой была разбросана по нескольким функциям и была достаточно сложной для записи.В подавляющем большинстве случаев достаточно изменить типы аргументов и возвращаемых значений для всех функций цепочки. Реже придется делать перегруженную версию (как, например, для sin и cos).
Если вы про раздел 2.3.2, то там, насколько я понял, речь про символьное дифференцирование — запишем выражение строкой (ну или списком), и будем работать с этой строкой.
Автоматическое дифференцирование — несколько другой подход, и, пожалуй, ответ на ваш второй вопрос довольно хорошо иллюстрирует эту разницу. Как решить исходную задачу? Как и написано в посте,
Пример (очень условный и выдуманный мною только что; компилируемость и работоспособность не проверял, но идея, думаю, понятна): пусть был код
и нам надо получить df/dx.
Заменяем везде, где речь идет про вычисляемые функции, double на Derivable:
Ну еще придется или перегрузить
Автоматическое дифференцирование — несколько другой подход, и, пожалуй, ответ на ваш второй вопрос довольно хорошо иллюстрирует эту разницу. Как решить исходную задачу? Как и написано в посте,
главное везде все промежуточные переменные (и в том числе результаты промежуточных функций) сделать Derivable, остальной код даже менять не потребуется! Пример (очень условный и выдуманный мною только что; компилируемость и работоспособность не проверял, но идея, думаю, понятна): пусть был код
double foo(double x, double y) {
if (x < y) return x * x - y;
else return y * y - x;
}
double bar(double x, double c) {
return sin(x / c) * c / x;
}
double f(double x, double Q, double t) {
double u = foo(x, Q * x / 2);
double v = bar(x / t, constant::c); // скорость света
return u / v;
}
и нам надо получить df/dx.
Заменяем везде, где речь идет про вычисляемые функции, double на Derivable:
Derivable foo(Derivable x, Derivable y) {
if (x < y) return x * x - y;
else return y * y - x;
}
Derivable bar(Derivable x, double c) { // считаем, что тут c никогда не зависит от x
return sin(x / c) * c / x;
}
Derivable f(double x, double Q, double t) {
Derivable u = foo(x, Q * x / 2);
Derivable v = bar(x / t, constants::c); // скорость света
return u / v;
}
Ну еще придется или перегрузить
operator<, или сделать приведение типов от Derivable обратно к double.Ну еще придется или перегрузить operator<, или сделать приведение типов от Derivable обратно к double.Вот второе делать нельзя ни в коем случае, иначе компилятор «проглотит» без всяких ошибок вызов немодифицированной функции:
Derivable foo(Derivable x, Derivable y) {
if (x < y) return x * x - y;
else return y * y - x;
}
double bar(double x, double c) { // Мы забыли преобразовать bar...
return sin(x / c) * c / x;
}
Derivable f(double x, double Q, double t) {
Derivable u = foo(x, Q * x / 2);
Derivable v = bar(x / t, constants::c); // но компилятор "догадался" сделать два пользовательских приведения типа
return u / v;
}
Кстати, комментарии вида «считаем, что тут c никогда не зависит от x» — тоже потенциальный источник ошибки. Если делать автоматическое дифференцирование, то надо менять типы всех аргументов, а там уже программа разберется. В конце концов, может, нам когда-нибудь захочется найти производную по другой переменной?
Так что правильный код будет таким:
Derivable foo(Derivable x, Derivable y) {
if (x < y) return x * x - y;
else return y * y - x;
}
Derivable bar(Derivable x, Derivable c) {
return sin(x / c) * c / x;
}
Derivable f(Derivable x, Derivable Q, Derivable t) {
Derivable u = foo(x, Q * x / 2);
Derivable v = bar(x / t, constants::c); // скорость света
return u / v;
}
// Вызывать так:
// Derivable f0 = f(Derivable::IndependendVariable(x0), Q, t)
Спасибо за развернутый ответ. Невнимательно просмотрев статью я подумал что дифференцирование происходит на стадии компиляции. Здесь применен другой любопытный приём.
Поправка: в коде таки баг, надо как минимум
А вообще, mayorovp выше написал еще правильнее.
xd = Derivable::IndependendVariable(x);
Derivable u = foo(xd, Q * xd / 2);
Derivable v = bar(xd / t, constants::c); // скорость света
А вообще, mayorovp выше написал еще правильнее.
На русском языке тоже есть литература, например «Символьный C++. Введение в компьютерную алгебру с использованием объектно-ориентированного программирования». А вообще дифференцирование — простая алгоритмическая процедура, так что тут можно и без книжек справиться.
Спасибо, автор. Очень хороший пример красивого использования возможностей C++. Один только вопрос. Многие функции, которые заданы путем «склейки» других функций на разных участках своего определения, недифференцируемы в точках склейки. Отрабатывает ли ваш пакет такие ситуации корректно? Что будет, если задать, например, функцию модуля Derivable abs(Derivable x)?
Пакет обработает такие ситуации корректно во всех точках, кроме точек излома.
Пример с abs будет выглядеть так (при решении «в лоб»):
Если пытаться исправить поведение в точках излома, то только так:
Пример с abs будет выглядеть так (при решении «в лоб»):
Derivable abs(Derivable x) {
if (x < 0)
return -x;
else
return x;
}
Если пытаться исправить поведение в точках излома, то только так:
Derivable abs(Derivable x) {
if (fabs(x.getValue() < eps)
return Derivable(0.0, NaN);
if (x < 0)
return -x;
else
return x;
}
Если я правильно понял, этот метод не заработает даже для композиции элементарных функций (например, sin(x^2)). В смысле, возможно работать будет, но выдавать будет лажу. Это если не пытаться рассматривать интегралы от элементарных функций и, как частный случай, спец. функции. Вы бы хоть от математической постановки задачи стартовали, а то фиг знает какие у вас тут грабли спрятаны.
И неясно, зачем было изобретать велосипед, если есть пакеты для символьных манипуляций, к тому же с трансляцией в байт-код (аналог JIT). Отказ от рефакторинга может быть слабым аргументом против принциальных ограничений выбранного подхода, но задачу, опять же, вы не сформулировали.
И неясно, зачем было изобретать велосипед, если есть пакеты для символьных манипуляций, к тому же с трансляцией в байт-код (аналог JIT). Отказ от рефакторинга может быть слабым аргументом против принциальных ограничений выбранного подхода, но задачу, опять же, вы не сформулировали.
Неправильно вы все поняли. Это метод прекрасно работает для композиции элементарных функций, и sin(x^2) подсчитает абсолютно правильно.
Конечно, если учесть правило дифференцирования сложной функции, будет работать. А теперь посмотрите в код топикстартера.
Собственно, в посте теория с практикой небрежно друг с другом связаны, в следующем смысле: если было достаточно решить конкретное уравнение (в котором композиции нет), то вопросов нет. Хотя получается, что вся ответственность за вывод формул на бумажке и все возможные при этом глюки переносятся на программиста.
Собственно, в посте теория с практикой небрежно друг с другом связаны, в следующем смысле: если было достаточно решить конкретное уравнение (в котором композиции нет), то вопросов нет. Хотя получается, что вся ответственность за вывод формул на бумажке и все возможные при этом глюки переносятся на программиста.
Пример с синусом работать не будет только потому, что код для синуса не реализован. Но если вы захотите посчитать производную cos(x^2), то всё замечательно сработает:
Derivable xd = Derivable::IndependendVariable(x);
Derivable xsq = xd*xd;
Derivable res = cos(xsq);
Получим:
xd=(x,1)
xsq=(x*x,2*x)
res=(cos(x*x),-sin(x*x)*2*x)
И что здесь не так?
Если вам нужен код для синуса (экспоненты, логарифма...) то он пишется по аналогии с косинусом в три строчки.
Derivable xd = Derivable::IndependendVariable(x);
Derivable xsq = xd*xd;
Derivable res = cos(xsq);
Получим:
xd=(x,1)
xsq=(x*x,2*x)
res=(cos(x*x),-sin(x*x)*2*x)
И что здесь не так?
Если вам нужен код для синуса (экспоненты, логарифма...) то он пишется по аналогии с косинусом в три строчки.
Спасибо, теперь понятнее. Помогло то, что вы сначала явно определяете xsq, а потом вручную вычисляете косинус от xsq.
Главная претензия у меня следующая: код топикстартера синтаксически поддерживает аддитивность и правило Лейбница (произведение производных), но нет синтаксической поддержки для композиции функций. Вы в своем примере композицию реализовывали руками. Как вы понимаете, это открывает широкое поле для всевозможных багов, особенно когда формулы становятся хоть сколько-нибудь сложными.
На всякий случай уточню — у меня есть некий опыт работы с системами компьютерной алгебры (CAS), включая такие экзотические как FORM. В большинстве CAS проблем из предыдущего абзаца нет в принципе, но в FORM все специфично. У него синтаксической поддержки для взятия производных нет, но возможно написать цикл подстановок, эквивалентный символьному дифференцированию. Фактически, приходится вручную прописывать аддитивность производной, правило Лейбница и производную композиции. В общем, это поучительный был опыт. А тут топикстартер наступает на грабли и главное не пытается вслух и громко сказать о глючевых ограничениях его дизайна.
Возвращаясь к миру быстрых вычислений (CAS в общем случае медленные), я встречал код для символьных манипуляций, специально заточенный под задачу «один раз символьно проманипулировать, потом очень много раз быстро вычислять». По сравнению с символьным движком предлагаемый стиль программирования выглядит ну… велосипедом с квадратными колесами. По крайней мере до того, как сформулированы критерии дизайна, из-за которых готовые решения не годятся. Да и начиналась статья в духе «поставили задачу, решил и вдруг потом открыл для себя новый мир аналогичных решений».
Главная претензия у меня следующая: код топикстартера синтаксически поддерживает аддитивность и правило Лейбница (произведение производных), но нет синтаксической поддержки для композиции функций. Вы в своем примере композицию реализовывали руками. Как вы понимаете, это открывает широкое поле для всевозможных багов, особенно когда формулы становятся хоть сколько-нибудь сложными.
На всякий случай уточню — у меня есть некий опыт работы с системами компьютерной алгебры (CAS), включая такие экзотические как FORM. В большинстве CAS проблем из предыдущего абзаца нет в принципе, но в FORM все специфично. У него синтаксической поддержки для взятия производных нет, но возможно написать цикл подстановок, эквивалентный символьному дифференцированию. Фактически, приходится вручную прописывать аддитивность производной, правило Лейбница и производную композиции. В общем, это поучительный был опыт. А тут топикстартер наступает на грабли и главное не пытается вслух и громко сказать о глючевых ограничениях его дизайна.
Возвращаясь к миру быстрых вычислений (CAS в общем случае медленные), я встречал код для символьных манипуляций, специально заточенный под задачу «один раз символьно проманипулировать, потом очень много раз быстро вычислять». По сравнению с символьным движком предлагаемый стиль программирования выглядит ну… велосипедом с квадратными колесами. По крайней мере до того, как сформулированы критерии дизайна, из-за которых готовые решения не годятся. Да и начиналась статья в духе «поставили задачу, решил и вдруг потом открыл для себя новый мир аналогичных решений».
Если бы я написал
res=cos(Derivable::IndependendVariable(x)*Derivable::IndependendVariable(x));
сработало бы всё точно так же. Просто воспринимайте Derivable как начало ряда Тейлора (f0+f1*dx), считая dx*dx=0 — и это будет обычное кольцо, на котором действует вся арифметика. А заодно и более сложные функции. Композицию, как таковую, никто не реализовывал: реализуется не f(g(x)), а элементарная функция от Derivable.
Хотя я могу написать
И замечательно посчитается производная сложной функции. Безо всяких усилий с нашей стороны. Потому что производная — это не магия, а арифметика. И производная сложной функции — это всего лишь подстановка одной линейной функции в другую.
Можно написать какой-нибудь Derivable2, в котором держать не только первую, но и вторую производную. И тоже не нужно будет прописывать производную суперпозиции в общем виде, достаточно написать её для 6 основных функций (sin, cos, exp, ln, a*b, a/b) — остальное сделает математика.
res=cos(Derivable::IndependendVariable(x)*Derivable::IndependendVariable(x));
сработало бы всё точно так же. Просто воспринимайте Derivable как начало ряда Тейлора (f0+f1*dx), считая dx*dx=0 — и это будет обычное кольцо, на котором действует вся арифметика. А заодно и более сложные функции. Композицию, как таковую, никто не реализовывал: реализуется не f(g(x)), а элементарная функция от Derivable.
Хотя я могу написать
Derivable f(Derivable x,double a,double b){
return a*cos(x)+b*sin(x);
}
Derivable res=f(f(Derivable.IndependendVariable(x),2,3),-3,-2);
И замечательно посчитается производная сложной функции. Безо всяких усилий с нашей стороны. Потому что производная — это не магия, а арифметика. И производная сложной функции — это всего лишь подстановка одной линейной функции в другую.
Можно написать какой-нибудь Derivable2, в котором держать не только первую, но и вторую производную. И тоже не нужно будет прописывать производную суперпозиции в общем виде, достаточно написать её для 6 основных функций (sin, cos, exp, ln, a*b, a/b) — остальное сделает математика.
Это забавно. Я вам говорю, что в предлагаемом решении приходится композицию функций реализовывать руками и это минус, а вы в ответ уверяете, что все замечательно работает. Да работать может даже если на ассемблере будете писать, но дык мой комментарий же не об этом.
Где вы видите «реализацию композиции руками»? То, что написано f(f(x))? Так в любом языке чтобы получить композицию, надо её написать. И внутри f() у меня нигде не используется формула производной сложной функции. Потому что в этом подходе она не нужна.
Извините, а вы зачем эту беседу ведете? Чего добиваетесь? Не хотите послушать рассказ о том, что бывает иначе чем вы привыкли — ну и не слушайте. Зачем мое время тратите?
Если всмотритесь в код вот в этом каменте (см ниже, habrahabr.ru/company/intel/blog/170729/#comment_5937833), увидите что композицию функций и правда можно не реализовывать руками, а дать специальной библиотеке символьных вычислений сделать свое дело. Причем сделать его гарантированно правильно, для любых выражений. В том числе тех, которые вам прямо сейчас не нужны, а через день потребуют рефакторинга вашего писаного самостоятельно кода — когда окажется, что вам и вторая производная нужна, и спец. функции бы неплохо добавить, и не ОДУ, а уравнения в частных производных и т.п.
Если всмотритесь в код вот в этом каменте (см ниже, habrahabr.ru/company/intel/blog/170729/#comment_5937833), увидите что композицию функций и правда можно не реализовывать руками, а дать специальной библиотеке символьных вычислений сделать свое дело. Причем сделать его гарантированно правильно, для любых выражений. В том числе тех, которые вам прямо сейчас не нужны, а через день потребуют рефакторинга вашего писаного самостоятельно кода — когда окажется, что вам и вторая производная нужна, и спец. функции бы неплохо добавить, и не ОДУ, а уравнения в частных производных и т.п.
Извините, а вы зачем эту беседу ведете?
Чтобы узнать что-то новое. В данном случае это удалось — про библиотеку СКА, встраиваемую в C++ я если и слышал, то давно забыл. Правда, сейчас не могу представить себе ситуацию, когда её было бы оправдано использовать — чтобы задача была достаточно сложной, чтобы не писать её на коленке, а проще было подключить чужой код, разобравшись в схеме подключения, его требованиях и т.п. — но при этом чтобы быть уверенным, что система точно справится с задачей и не скажет «ну, не смогла». Применение этого механизма вместо класса на 50 строчек (в случае, что больше в программе символьные вычисления не применяются) выглядит несколько громоздким. Но, наверное, когда понадобится работа с какими-нибудь дистрибутивными полиномами с алгебраическими коэффициентами, то этот генератор и интерпретатор байт-кода будет взять дешевле.
Эх. Если бы вы этот текст сказали с самого начала, беседа пошла совсем по иному пути. У меня не складывается впечатления, что вы хотите узнать что-то новое. Вы достаточно агрессивно отстаиваете точку зрения что предлагаемое и так работает. Вы и так уверены в том, что оно работает, зачем защищиться? А вот если бы вы хотели узнать новое, то просто переспросили бы пару раз вместо нападения.
Заметим, что этот метод (с подсчётом производной через Derivable) будет работать, даже если считать значение функции каким-нибудь итеративным алгоритмом, или другими численными методами. Пакеты символьных вычислений быстро скажут «не могу» из-за экспоненциального роста памяти при операциях, более сложных, чем сложение.
Кстати, что это за пакеты с трансляцией в байт-код? Сколько они стоят? Полученный код легко встраивать в другие программы и переносить отдельно от самого пакета?
Кстати, что это за пакеты с трансляцией в байт-код? Сколько они стоят? Полученный код легко встраивать в другие программы и переносить отдельно от самого пакета?
эммм, чего? Пакеты символьных вычислений будут экспоненциально быстро накапливать данные, если эти данные сохранять и не перезаписывать. В каждой уважающей себя системе компьютерной алгебры есть средство «забывания» ненужных переменных. Ну или можно просто присваивать новые значения на место старых.
По ссылкам будет скорее всего не то, чего вы ожидаете, но все же:
code.google.com/p/numexpr/
beltoforion.de/muparsersse/math_expression_compiler_en.html
muparser.beltoforion.de/mup_features.html#idDef2
www.codeproject.com/Articles/45797/Fast-Polymorphic-Math-Parser
Это из первых хитов поиска на «fast math formula evaluation JIT», хотя в сторону muparsersse я в какой-то момент смотрел.
Если интересны более продвинутые C++ движки символьных вычислений, но без JIT и фокуса на «много раз вычислить», то можно посмотреть на
www.ginac.de/
и
Ev3:
www.lix.polytechnique.fr/~liberti/Ev3.pdf
www.lix.polytechnique.fr/~liberti/Ev3-1.0.tar.gz
А кроме того, можно же ж и комбинировать: сначала символьно поманипулировать внутри ginac/Ev3, потом отдать быстровычислятелю в заданных точках.
По ссылкам будет скорее всего не то, чего вы ожидаете, но все же:
code.google.com/p/numexpr/
beltoforion.de/muparsersse/math_expression_compiler_en.html
muparser.beltoforion.de/mup_features.html#idDef2
www.codeproject.com/Articles/45797/Fast-Polymorphic-Math-Parser
Это из первых хитов поиска на «fast math formula evaluation JIT», хотя в сторону muparsersse я в какой-то момент смотрел.
Если интересны более продвинутые C++ движки символьных вычислений, но без JIT и фокуса на «много раз вычислить», то можно посмотреть на
www.ginac.de/
и
Ev3:
www.lix.polytechnique.fr/~liberti/Ev3.pdf
www.lix.polytechnique.fr/~liberti/Ev3-1.0.tar.gz
А кроме того, можно же ж и комбинировать: сначала символьно поманипулировать внутри ginac/Ev3, потом отдать быстровычислятелю в заданных точках.
Ну, допустим, решаете вы численно какой-нибудь диффур, скажем, f'=x-f*f (взял с потолка — мог написать что угодно, лишь бы не линейное). Простейшим способом:
И вздумалось вам посчитать производную результата по a. Заменив double на Derivable, вы сразу получите и значение, и производную. А если захочется найти эту производную символьно (оставив тот же цикл для решения уравнения)? Какой степени получится многочлен на выходе из цикла?
Ссылки посмотрю, спасибо.
double f=a;
int N=1000000;
double dx=1.0/N;
for(int k=0;k<N;k++){
f+=(k*dx-f*f)*dx;
}
И вздумалось вам посчитать производную результата по a. Заменив double на Derivable, вы сразу получите и значение, и производную. А если захочется найти эту производную символьно (оставив тот же цикл для решения уравнения)? Какой степени получится многочлен на выходе из цикла?
Ссылки посмотрю, спасибо.
Вы видимо себе представляете процесс символьного вычисления производной единственно возможным способом как построение скобочной структуры символьных итераций первых производных, которые затем символьно упрощаются (получаем многочлен), и уже потом вычисляется значение многочлена в точке. В такой схеме, разумеется, выражение распухнет.
Если уж предлагать все варианты, которые предлагает символьный счет, то их я вижу как минимум два:
— построение точного решения, то есть не численно. Для данного случайного диффура есть точное решение. В точке его вычислить — легко. При известных оговорках, в определенных случаях это может быть быстрее численного счета. Но я его ввиду не имел, а имел следующий:
— однажды совершить символьную манипуляцию, эквивалентную вашей итерации, перегнать результат в байт-код и на каждом шаге цикла быстро вычислять. Получится тот же численный процесс, сходящийся к производной решения ОДУ. Только в этом случае вы будете использовать инструментарий, который вам точно не позволит выстрелить в ногу.
А о том, как растут символьные данные я в курсе, у меня были символьные вычисления с полтерабайта данных для переопределенных систем.
Если уж предлагать все варианты, которые предлагает символьный счет, то их я вижу как минимум два:
— построение точного решения, то есть не численно. Для данного случайного диффура есть точное решение. В точке его вычислить — легко. При известных оговорках, в определенных случаях это может быть быстрее численного счета. Но я его ввиду не имел, а имел следующий:
— однажды совершить символьную манипуляцию, эквивалентную вашей итерации, перегнать результат в байт-код и на каждом шаге цикла быстро вычислять. Получится тот же численный процесс, сходящийся к производной решения ОДУ. Только в этом случае вы будете использовать инструментарий, который вам точно не позволит выстрелить в ногу.
А о том, как растут символьные данные я в курсе, у меня были символьные вычисления с полтерабайта данных для переопределенных систем.
— однажды совершить символьную манипуляцию, эквивалентную вашей итерации, перегнать результат в байт-код и на каждом шаге цикла быстро вычислять. Получится тот же численный процесс, сходящийся к производной решения ОДУ. Только в этом случае вы будете использовать инструментарий, который вам точно не позволит выстрелить в ногу.
Боюсь, что это «однажды» придётся повторять при любом изменении кода программы, вычисляющего функции, производные которых надо считать. И, поскольку в операции используется несколько разных инструментов, и код содержится в нескольких местах — придётся ведь, например, как-нибудь привязывать к байт-коду значения дополнительных параметров функции — а значит, и система символьных вычислений, и программа на C должны знать про данный фрагмент, — то разработка и развитие такой программы выглядит чрезмерно сложным делом. И шансов ошибиться в этой процедуре — при каждом новом изменении байт-кода! — намного больше, чем в «велосипеде», написанном один раз, причём ровно в той среде, где он впоследствии будет работать.
Мне кажется, вы столкнулись с непривычным вам решением и совершенно не прилагаете усилий понять что вам предлагают. Сходили бы по ссылкам-то, для начала, почитали что там пишут. Ну и про JIT тоже почитали бы.
Грубо говоря, я предлагаю сделать следующее:
Тут я псевдокод написал, но должно быть понятно, о чем я:
— перед генерацией байт-кода происходит автоматическое вычисление производной от символьного выражения, эквивалентного одному шагу численного алгоритма,
— для один раз созданного выражения генерация байт-кода происходит автоматически,
— внутри горячего цикла происходит выполнение бинарного кода, вычисляющего один шаг численного алгоритма, сходящегося к производной решения ОДУ.
При этом interpreter у вас уже есть и поддерживается другой командой разработчиков. Никакого «однажды при любом изменении кода программы» у вас нет — однажды добавили сторонюю библиотеку в проект и поехали.
Грубо говоря, я предлагаю сделать следующее:
<code>
BytecodeObject symbolic_iteration = interpreter.generate_bytecode(interpreter.evaluate("diff(k*dx-f*f)*dx,x)") );
double f=a;
double df = 0;
int N=1000000;
double dx=1.0/N;
for(int k=0;k<N;k++)
{
f+=(k*dx-f*f)*dx;
df += symbolic_iteration(k,dx,f);
}
</code>
Тут я псевдокод написал, но должно быть понятно, о чем я:
— перед генерацией байт-кода происходит автоматическое вычисление производной от символьного выражения, эквивалентного одному шагу численного алгоритма,
— для один раз созданного выражения генерация байт-кода происходит автоматически,
— внутри горячего цикла происходит выполнение бинарного кода, вычисляющего один шаг численного алгоритма, сходящегося к производной решения ОДУ.
При этом interpreter у вас уже есть и поддерживается другой командой разработчиков. Никакого «однажды при любом изменении кода программы» у вас нет — однажды добавили сторонюю библиотеку в проект и поехали.
У меня выше со скобками неверно, правильно читать так:
BytecodeObject symbolic_iteration = interpreter.generate_bytecode(interpreter.evaluate("diff( (k*dx-f*f)*dx, x )") );
Один вопрос. Чем ваш подход лучше следующего кода:
Он что, работает быстрее? В написании/поддержке проще? Дает результат правильнее?
Я, кстати, несколько сомневаюсь насчет правильности результата, полученного через Derivable, для функций, вычисленных не аналитически, а численно. Вот это сравнение увидеть и хотелось бы, а не то, что вы написали.
Derivable f=a;
double df = 0;
int N=1000000;
double dx=1.0/N;
for(int k=0;k<N;k++)
f += (k*dx-f*f)*dx;
Он что, работает быстрее? В написании/поддержке проще? Дает результат правильнее?
Я, кстати, несколько сомневаюсь насчет правильности результата, полученного через Derivable, для функций, вычисленных не аналитически, а численно. Вот это сравнение увидеть и хотелось бы, а не то, что вы написали.
В разы проще и быстрее в поддержке:
— нет класса Derivable
— нет рефакторинга, когда нужно добавить поддержку, например, второй производной.
— композиция функций и другие естественные конструкции инкапсулированы в строке, которую обрабатывает внешняя специализированная библиотека. Оно не размазано по коду в форме
трех последовательных присваиваний, за порядком которых приходится следить и сверяться с формулой на бумаге. Цитируемый код заменяется одной строкой в духе
— в итоге дает меньше шансов выстрелить себе в ногу
— и вообще, всякое явное представление концепции лучше, чем более низкоуровневое описание ее следствий (формула целиком против алгоритма, который ее вычисляет)
— работает, наверное, со сравнимой скоростью внутри горячего цикла, медленнее на старте (накладные расходы на символьные манипуляции) — это зависит от библиотеки, которую вы возьмете. Со сравнимой скоростью должен работать JIT подход, чисто символьное вычисление, конечно, будет медленным.
— работает одинаково хорошо для простых примеров, так и для сложных. Вторая производная, производная спец. функций, попытка аналитически решить уравнение до решения численно — это возможно, и инкапсулировано в специализированном пакете, который для того и писан. Посмотрите на возможности GiNaC — вы вряд ли хотите переписывать его функционал.
Если вы сомневаетесь в точности численного ответа, то у вас к вариантам сменить численный метод прибавляется еще вариант попробовать решить точно, символьно. Этот ответ лежит на поверхности. И я его тут даже указывал, в других комментариях.
Вообще-то в этом треде я отвечал юзеру Mrrl, который хотел именно численного метода для производной. Он хотел именно этого — я показал, что численно сойтись к производной возможно.
— нет класса Derivable
— нет рефакторинга, когда нужно добавить поддержку, например, второй производной.
— композиция функций и другие естественные конструкции инкапсулированы в строке, которую обрабатывает внешняя специализированная библиотека. Оно не размазано по коду в форме
Derivable xd = Derivable::IndependendVariable(x);
Derivable xsq = xd*xd;
Derivable res = cos(xsq);трех последовательных присваиваний, за порядком которых приходится следить и сверяться с формулой на бумаге. Цитируемый код заменяется одной строкой в духе
interpreter.evaluate("cos(x^2)", "x = 1");
— в итоге дает меньше шансов выстрелить себе в ногу
— и вообще, всякое явное представление концепции лучше, чем более низкоуровневое описание ее следствий (формула целиком против алгоритма, который ее вычисляет)
— работает, наверное, со сравнимой скоростью внутри горячего цикла, медленнее на старте (накладные расходы на символьные манипуляции) — это зависит от библиотеки, которую вы возьмете. Со сравнимой скоростью должен работать JIT подход, чисто символьное вычисление, конечно, будет медленным.
— работает одинаково хорошо для простых примеров, так и для сложных. Вторая производная, производная спец. функций, попытка аналитически решить уравнение до решения численно — это возможно, и инкапсулировано в специализированном пакете, который для того и писан. Посмотрите на возможности GiNaC — вы вряд ли хотите переписывать его функционал.
Если вы сомневаетесь в точности численного ответа, то у вас к вариантам сменить численный метод прибавляется еще вариант попробовать решить точно, символьно. Этот ответ лежит на поверхности. И я его тут даже указывал, в других комментариях.
Вот это сравнение увидеть и хотелось бы, а не то, что вы написали.
Вообще-то в этом треде я отвечал юзеру Mrrl, который хотел именно численного метода для производной. Он хотел именно этого — я показал, что численно сойтись к производной возможно.
— композиция функций и другие естественные конструкции инкапсулированы в строке, которую обрабатывает внешняя специализированная библиотека. Оно не размазано по коду в форме
трех последовательных присваиваний, за порядком которых приходится следить и сверяться с формулой на бумаге. Цитируемый код заменяется одной строкой в духеDerivable xd = Derivable::IndependendVariable(x); Derivable xsq = xd*xd; Derivable res = cos(xsq);
interpreter.evaluate("cos(x^2)", "x = 1");
Вы до сих пор не согласны, что следующий код тоже работает:
xd = Derivable::IndependendVariable(x);
return cos(xd*xd);
?
разумеется я с этим согласен. А остальные мои аргументы вам не важны?
Ок.
И вторая производная там же есть, и все, что надо. И, как верно Mrrl написал ниже, добавить код для вычисления второй производной элементарных функций достаточно просто. Ну или воспользоваться
А если и функция задается в силу своего физического смысла не одной строкой, а пятью (см. пример в комментах выше), то прострелить ногу из-за необходимости писать больше кода там еще проще.
С этой точки зрения не вижу разницы между
— нет класса DerivableКласс
— нет рефакторинга, когда нужно добавить поддержку, например, второй производной.
Derivable или его аналог есть, только не в вашем коде, а в том, который вы взяли. Вопрос о том, что лучше — готовые решения или самописные — в каждом случае должен решаться отдельно, но (и я вам это уже отвечал в другом комменте) и по автоматическому дифференцированию есть готовые библиотеки.И вторая производная там же есть, и все, что надо. И, как верно Mrrl написал ниже, добавить код для вычисления второй производной элементарных функций достаточно просто. Ну или воспользоваться
Derivable<Derivable<double> >, что еще проще.— композиция функций и другие естественные конструкцииС этим разобрались.
— в итоге дает меньше шансов выстрелить себе в ногуПо-моему, максимальное упрощение пользовательского кода, т.е. того, где вы вычисляете функцию, более существенно уменьшает шансы выстрелить в ногу. В вашем варианте выстрел в ногу получается элементарно, если из двух копий формулы забыть исправить одну.
А если и функция задается в силу своего физического смысла не одной строкой, а пятью (см. пример в комментах выше), то прострелить ногу из-за необходимости писать больше кода там еще проще.
— и вообще, всякое явное представление концепции лучше, чем более низкоуровневое описание ее следствий (формула целиком против алгоритма, который ее вычисляет)
С этой точки зрения не вижу разницы между
interpreter.evaluate("cos(x^2)+x^3-sqrt(x)", "x = 1");
xd = Derivable::IndependendVariable(x);
return cos(xd * xd) + xd*xd*xd - sqrt(xd);
— работает, наверное, со сравнимой скоростью внутри горячего цикла, медленнее на старте (накладные расходы на символьные манипуляции)не в те ворота
— это зависит от библиотеки, которую вы возьмете. Со сравнимой скоростью должен работать JIT подход, чисто символьное вычисление, конечно, будет медленным.
— работает одинаково хорошо для простых примеров, так и для сложных. Вторая производная, производная спец. функций, попытка аналитически решить уравнение до решения численно — это возможно, и инкапсулировано в специализированном пакете, который для того и писан. Посмотрите на возможности GiNaC — вы вряд ли хотите переписывать его функционал.Все, кроме аналитического решения, вполне реализуется и доступно в имеющихся библиотеках автоматического дифференцирования. Аналитическое решение — конечно нет, но это уже совсем другая задача.
Так вы бы с задачи начали — и этого разговора могло б не произойти.
И вторая производная там же есть, и все, что надо. И, как верно Mrrl написал ниже, добавить код для вычисления второй производной элементарных функций достаточно просто. Ну или воспользоваться Derivable<Derivable >, что еще проще
Итак, вы согласились с необходимостью рефакторинга.
— композиция функций и другие естественныеконструкции
С этим разобрались.
Нет, не разобрались. Как насчет естественной конструкции посчитать интеграл от функции? Если такая задача встанет, вам придется добавлять код численного интегрирования. В предлагаемом подходе — изменить формулу один раз.
По-моему, максимальное упрощение пользовательского кода, т.е. того, где вы вычисляете функцию, более существенно уменьшает шансы выстрелить в ногу. В вашем варианте выстрел в ногу получается элементарно, если из двух копий формулы забыть исправить одну.
Да, это возможный сценарий — исправить формулу в одном месте. С другой стороны, объект символьного движка копируем — зачем в двух местах держать одну и ту же формулу?
А если и функция задается в силу своего физического смысла не одной строкой, а пятью (см. пример в комментах выше), то прострелить ногу из-за необходимости писать больше кода там еще проще.
Если вы об этом примере
Derivable foo(Derivable x, Derivable y) { if (x < y) return x * x - y; else return y * y - x; } Derivable bar(Derivable x, Derivable c) { return sin(x / c) * c / x; } Derivable f(Derivable x, Derivable Q, Derivable t) { Derivable u = foo(x, Q * x / 2); Derivable v = bar(x / t, constants::c); // скорость света return u / v; } // Вызывать так: // Derivable f0 = f(Derivable::IndependendVariable(x0), Q, t)
то ответом будет: символьный движок позволяет записать несколько символьных выражений за одну «сессию» + поддерживать кусочно определенные функции. Так вы еще и от С++ синтаксического сахара избавитесь в виде функций foo и bar.
Наконец, вы не видите разницы между
interpreter.evaluate("cos(x^2)+x^3-sqrt(x)", "x = 1");
и
xd = Derivable::IndependendVariable(x); return cos(xd * xd) + xd*xd*xd - sqrt(xd);
А как насчет
interpreter.evaluate("int(cos(x^2)+x^3-sqrt(x),0,y", "y = 1");
? Имеется ввиду интеграл с переменным верхним пределом.
— нет рефакторинга, когда нужно добавить поддержку, например, второй производной.
Интересное замечание. Сколько усилий вам потребуется, чтобы добавить вычисление второй производной для функции, вычисление которой занимает несколько методов и циклов? Через Derivable это делается элементарно — меняется несколько строчек в самом классе, а вся остальная программа не трогается вообще. А сможете ли вы посчитать вторую производную по a для того же решения дифференциального уравнения — не меняя тела цикла, в котором оно решается?
если функция исходно задается аналитически, то добавить вычисление второй производной — одна строчка. Если аналитической формы нет — тогда увы. Собственно, ради этого я с самого начала спрашивал математическую постановку задачи. Без нее неясно, есть у вас неаналитически заданные функции или нет.
На последний вопрос — ответ да. Байт-код, вычисляющий вторую производную будет инкапсулирован в объекте symbolic_iteration ровно так же, как и байт-код, вычисляющий первую производную. Вам только предстоит поменять строчное представление формулы, передаваемой символьному движку.
На последний вопрос — ответ да. Байт-код, вычисляющий вторую производную будет инкапсулирован в объекте symbolic_iteration ровно так же, как и байт-код, вычисляющий первую производную. Вам только предстоит поменять строчное представление формулы, передаваемой символьному движку.
Вы пишите:
В исходном посте:
В этом случае (вот условный пример: habrahabr.ru/company/intel/blog/170729/#comment_5934705 ) тоже будет достаточно одной строчки? Или надо будет каждую из функций переписывать, еще и дублируя запись каждой функции?
если функция исходно задается аналитически, то добавить вычисление второй производной — одна строчка.
В исходном посте:
Эта левая часть у меня была довольно сложная — даже просто вычисление ее значений в программе было разбросано по нескольким функциям,
В этом случае (вот условный пример: habrahabr.ru/company/intel/blog/170729/#comment_5934705 ) тоже будет достаточно одной строчки? Или надо будет каждую из функций переписывать, еще и дублируя запись каждой функции?
Смотрю я на этот код:
interpreter.generate_bytecode(interpreter.evaluate(«diff( (k*dx-f*f)*dx, x )») );
и совершенно не понимаю, как он работает. Дифференцируется некоторое выражение (k*dx-f*f)*dx по переменной x (кстати, почему не по a)? Откуда интерпретатор знает, что в этом выражении зависит от a, а что нет? Даже если вы напишете
diff((k*dx-f(a)*f(a))*dx,a) — то для вычисления этого выражения интерпретатору придётся знать уже вычисленное выражение df/da — но вы его в код не передаёте (а если будете передавать — как объясните, какое из переданных чисел попадёт в этот аргумент?) Если будете вычислять вторую производную — то для вычисления интерпретатору понадобится знать k,dx,f(a),df/da и ddf/da^2 — причём про df/da заранее неизвестно (в общем случае, не анализируя формулу), потребуется оно или нет. В любом случае, придётся все три выражения считать отдельно и аккуратно передавать в байт-код. И такое изменение вы называете «нет рефакторинга»? Или для второй производной вы воспользуетесь каким-то более простым и понятным методом?
interpreter.generate_bytecode(interpreter.evaluate(«diff( (k*dx-f*f)*dx, x )») );
и совершенно не понимаю, как он работает. Дифференцируется некоторое выражение (k*dx-f*f)*dx по переменной x (кстати, почему не по a)? Откуда интерпретатор знает, что в этом выражении зависит от a, а что нет? Даже если вы напишете
diff((k*dx-f(a)*f(a))*dx,a) — то для вычисления этого выражения интерпретатору придётся знать уже вычисленное выражение df/da — но вы его в код не передаёте (а если будете передавать — как объясните, какое из переданных чисел попадёт в этот аргумент?) Если будете вычислять вторую производную — то для вычисления интерпретатору понадобится знать k,dx,f(a),df/da и ddf/da^2 — причём про df/da заранее неизвестно (в общем случае, не анализируя формулу), потребуется оно или нет. В любом случае, придётся все три выражения считать отдельно и аккуратно передавать в байт-код. И такое изменение вы называете «нет рефакторинга»? Или для второй производной вы воспользуетесь каким-то более простым и понятным методом?
Еще раз подчеркиваю, что весь мой код с interpreter условный, потому что нужно сперва определиться с инструментарием, затем четко ему следовать. Я вам накидал ссылок — вы почитали, определились?
Вы правы в том, что нужно как-то явно объявить зависимость от a. Способ такого объявления зависит от движка.
Вторую производную вы тут посчитаете так:
(сейчас использую аналог синтаксиса Maple, в общем случае зависит от движка)
Вы правы в том, что нужно как-то явно объявить зависимость от a. Способ такого объявления зависит от движка.
Вторую производную вы тут посчитаете так:
... diff((k*dx-f(a)*f(a))*dx,a,a)...
(сейчас использую аналог синтаксиса Maple, в общем случае зависит от движка)
Я посмотрел ссылки. Первые 4 дают всего лишь парсеры математических выражений — это не особо интересно, если у меня нет пользовательского ввода и разбора введенных строк.
GiNaC — работа с символьными объектами? Да, он может дать решение задачи, но совсем не там, где мы его ищем. Надо будет заменить double на ряды Тейлора/Лорана и аккуратно с ними работать, упрощая после каждой операции. Но это решение совершенно ортогонально тому, что вы описываете с символьным взятием второй производной — оно ближе к решению автора поста. Так что продолжим искать решение через diff.
Последний пакет — да, умеет работать с формулами, умеет дифференцировать, умеет подставлять значения переменных в формулы. Если надо будет сосчитать функцию, а потом найти её производную, это поможет. Но поможет ли в последней задаче, где надо дифференцировать выражение с символом функции, считая, что конкретные значения производных этой функции заданы где-то в переменных C (а сивольного выражения f мы вообще не считаем)? Беглый просмотр не дал такой возможности.
Наконец, «движок Maple». Я пока ещё не посмотрел, сколько будут стоить построенные на нём пакеты, но, учитывая, что даже свежая версия самого Maple мне не по карману (приходится пользоваться древней 6.0), боюсь что это будет немало. А лицензия на распространение программ, использующих этот движок — еще больше. Но предположим даже, что мы решились использовать его в коде.
Итак, записали мы выражение (пока речь не о С++, а о том, что могло бы происходить в движке)
q:=diff(f(a)+(k*dx-f(a)*f(a))*dx,a,a);
Maple раскрыл его в
q := diff(f(a),a,a)-2*diff(f(a),a)^2*dx-2*f(a)*diff(f(a),a,a)*dx
Самой этой формулы мы не видим!
Предположим, что мы смогли выполнить аналог операции
subs(f(a)=u,diff(f(a),a)=u1,diff(f(a),a,a)=u2,q)
(корректного способа сделать это в Maple я не нашел, а изучать его с самого начала, чтобы найти эту возможность, сейчас некогда) — и получить
u2-2*u1^2*dx-2*u*u2*dx.
Далее, построили байт-код вызовом вроде
— и наконец-то получили фрагмент байт-кода, пригодный для вызова и передачи ему нужных данных.
Аналогично, нам придётся написать
— ведь нам нужна первая производная, чтобы передать её в функцию для второй производной?
И только сейчас мы можем, наконец, написать наш цикл:
Разумеется, интерфейс всех вызовов условный.
По-моему, я ничего не забыл. И не сделал ничего лишнего, всё, что написано, необходимо, чтобы передать в движок нужную информацию — по крайней мере, на уровне моей фантазии. И мне такое решение не кажется простым и безопасным — шансов «отстрелить ногу» слишком много.
Какой интерфейс символьного движка может дать более удобное и надёжное решение?
GiNaC — работа с символьными объектами? Да, он может дать решение задачи, но совсем не там, где мы его ищем. Надо будет заменить double на ряды Тейлора/Лорана и аккуратно с ними работать, упрощая после каждой операции. Но это решение совершенно ортогонально тому, что вы описываете с символьным взятием второй производной — оно ближе к решению автора поста. Так что продолжим искать решение через diff.
Последний пакет — да, умеет работать с формулами, умеет дифференцировать, умеет подставлять значения переменных в формулы. Если надо будет сосчитать функцию, а потом найти её производную, это поможет. Но поможет ли в последней задаче, где надо дифференцировать выражение с символом функции, считая, что конкретные значения производных этой функции заданы где-то в переменных C (а сивольного выражения f мы вообще не считаем)? Беглый просмотр не дал такой возможности.
Наконец, «движок Maple». Я пока ещё не посмотрел, сколько будут стоить построенные на нём пакеты, но, учитывая, что даже свежая версия самого Maple мне не по карману (приходится пользоваться древней 6.0), боюсь что это будет немало. А лицензия на распространение программ, использующих этот движок — еще больше. Но предположим даже, что мы решились использовать его в коде.
Итак, записали мы выражение (пока речь не о С++, а о том, что могло бы происходить в движке)
q:=diff(f(a)+(k*dx-f(a)*f(a))*dx,a,a);
Maple раскрыл его в
q := diff(f(a),a,a)-2*diff(f(a),a)^2*dx-2*f(a)*diff(f(a),a,a)*dx
Самой этой формулы мы не видим!
Предположим, что мы смогли выполнить аналог операции
subs(f(a)=u,diff(f(a),a)=u1,diff(f(a),a,a)=u2,q)
(корректного способа сделать это в Maple я не нашел, а изучать его с самого начала, чтобы найти эту возможность, сейчас некогда) — и получить
u2-2*u1^2*dx-2*u*u2*dx.
Далее, построили байт-код вызовом вроде
diff2=interpreter.generate_bytecode(
"proc(k,dx,u,u1,u2) \
return subs(f(a)=u,diff(f(a),a)=u1,diff(f(a),a,a)=u2,\
diff(f(a)+(k*dx-f(a)*f(a))*dx,a,a);");
— и наконец-то получили фрагмент байт-кода, пригодный для вызова и передачи ему нужных данных.
Аналогично, нам придётся написать
diff1=interpreter.generate_bytecode(
"proc(k,dx,u,u1) \
return subs(f(a)=u,diff(f(a),a)=u1,\
diff(f(a)+(k*dx-f(a)*f(a))*dx,a);");
— ведь нам нужна первая производная, чтобы передать её в функцию для второй производной?
И только сейчас мы можем, наконец, написать наш цикл:
double f=a,f1=0,f2=0;
int N=1000000;
double dx=1.0/N;
for(int k=0;k<N;k++){
double f2_new=interpreter.evaluate(diff2,k,dx,f,f1,f2);
double f1_new=interpreter.evaluate(diff1,k,dx,f,f1);
f += (k*dx-f*f)*dx;
f1=f1_new;
f2=f2_new;
}
Разумеется, интерфейс всех вызовов условный.
По-моему, я ничего не забыл. И не сделал ничего лишнего, всё, что написано, необходимо, чтобы передать в движок нужную информацию — по крайней мере, на уровне моей фантазии. И мне такое решение не кажется простым и безопасным — шансов «отстрелить ногу» слишком много.
Какой интерфейс символьного движка может дать более удобное и надёжное решение?
Если я правильно понял ваше затруднение с рядами Тейлора/Лорана в том смысле, что вы на бегу не увидели, умеет ли ginac дифференцировать, эта ссылка может помочь: www.ginac.de/tutorial/Symbolic-differentiation.html
Не смогу сейчас детально прокомментировать Ev3, сам с ним разбирался полдня. Первые четыре ссылки и правда только на парсеры готовых выражений, но зато они (по крайней мере некоторые) решают задачу «дано фиксированное выражение, быстро и много раз посчитай значение в точке». Можно символьные манипуляции проводить одним движком, а вычисление — другим. Согласен, выглядит странно — но другого готового решения я не знаю.
Чур меня, вот предлагать использовать Maple для этой задачи я точно не собирался. Его я использовал лишь как пример того, что задача счета N-й производной — одно из базовых свойств почти любого движка компьютерной алгебры.
В общем и целом, вы написали именно ту схему решения, которую я предлагал — я рад, что по крайней мере есть общее понимание. Необходимость подставок предопределена взаимозависимостью двух численных задач, от них никуда не уйти в сопряженных задачах (adjoint problems). Решение топикстартера тоже подразумевает подстановку исходной функции в ее производную, когда выполняется f += (update) c переменной f типа Derivable.
Что можно улучшить? Можно продифференцировать выражение один раз, объявить это одним символьным объектом, а в diff1 и diff2 заниматься только подстановками, причем сразу подставив значение неизменяемых констант (например, dx), до interpreter.evalute(). Так получится сократить объем выполняемого байт-кода в горячем цикле.
Вы вольны принимать решения о дизайне вашего кода. Моей задачей было рассказать об альтернативном решении. Считаете, что так вы больше рискуете — ну, ок, ваше решение. Мне же кажется, что так вы явно отделяете сущности и пишете код более близким к его математическому смыслу. Тем оно и удобнее, и безопаснее в смысле адаптации к новым задачам.
Не смогу сейчас детально прокомментировать Ev3, сам с ним разбирался полдня. Первые четыре ссылки и правда только на парсеры готовых выражений, но зато они (по крайней мере некоторые) решают задачу «дано фиксированное выражение, быстро и много раз посчитай значение в точке». Можно символьные манипуляции проводить одним движком, а вычисление — другим. Согласен, выглядит странно — но другого готового решения я не знаю.
Чур меня, вот предлагать использовать Maple для этой задачи я точно не собирался. Его я использовал лишь как пример того, что задача счета N-й производной — одно из базовых свойств почти любого движка компьютерной алгебры.
В общем и целом, вы написали именно ту схему решения, которую я предлагал — я рад, что по крайней мере есть общее понимание. Необходимость подставок предопределена взаимозависимостью двух численных задач, от них никуда не уйти в сопряженных задачах (adjoint problems). Решение топикстартера тоже подразумевает подстановку исходной функции в ее производную, когда выполняется f += (update) c переменной f типа Derivable.
Что можно улучшить? Можно продифференцировать выражение один раз, объявить это одним символьным объектом, а в diff1 и diff2 заниматься только подстановками, причем сразу подставив значение неизменяемых констант (например, dx), до interpreter.evalute(). Так получится сократить объем выполняемого байт-кода в горячем цикле.
Вы вольны принимать решения о дизайне вашего кода. Моей задачей было рассказать об альтернативном решении. Считаете, что так вы больше рискуете — ну, ок, ваше решение. Мне же кажется, что так вы явно отделяете сущности и пишете код более близким к его математическому смыслу. Тем оно и удобнее, и безопаснее в смысле адаптации к новым задачам.
Затруднений с рядами Тейлора у меня не было — наоборот, я указал, что они дают решение задачи, но оно лежит не в той плоскости, в которой мы работаем. По ссылке, которую вы привели, снова находится дифференцирование композиции известных функций по данной переменной (и последующее вычисление ответа в конкретной точке), и нет ответа на главные вопросы:
— как объяснить движку, что в выражении f+(k*dx-f*f)*dx величина f зависит от параметра a, по которому мы ведем дифференцирование?
— как заставить движок объяснить нам, что для корректного вычисления diff(f+(k*dx-f*f)*dx,a,a) требуется задать не только k,dx и f, но и df/da и ddf/da^2 в данной точке?
— Какой движок позволит передать эти параметры интерпретатору байт-кода так, чтобы они попали в правильные места, без сверхусилий с нашей стороны? А то, боюсь, окажется, что для решения задачи потребуется сначала глазами посмотреть на формулу, которая получилась, потом понять, как привести её к виду, пригодному для представления в байт-коде, и лишь потом запрограммировать всю эту цепочку, оставив движку самое лёгкое — само символьное дифференцирование. И теряем мы больше, чем можем выиграть.
— как объяснить движку, что в выражении f+(k*dx-f*f)*dx величина f зависит от параметра a, по которому мы ведем дифференцирование?
— как заставить движок объяснить нам, что для корректного вычисления diff(f+(k*dx-f*f)*dx,a,a) требуется задать не только k,dx и f, но и df/da и ddf/da^2 в данной точке?
— Какой движок позволит передать эти параметры интерпретатору байт-кода так, чтобы они попали в правильные места, без сверхусилий с нашей стороны? А то, боюсь, окажется, что для решения задачи потребуется сначала глазами посмотреть на формулу, которая получилась, потом понять, как привести её к виду, пригодному для представления в байт-коде, и лишь потом запрограммировать всю эту цепочку, оставив движку самое лёгкое — само символьное дифференцирование. И теряем мы больше, чем можем выиграть.
… и вам снова пришлось явно разделить функцию и производную, и работать с ними отдельно — вместо того, чтобы пользоваться единым объектом «функция с приращением». Все операции усложнятся в несколько раз: вместо арифметической операции над этим объектом пишется тройка — операция над функцией, генерация байт-кода, вызов байт-кода. Причём код функции-таки пишется в двух местах — в точке вызова и в генерации байт-кода. И связь через параметры тоже придётся отслеживать. При любом изменении кода функции…
И мы с размаху вляпываемся во все прелести Not Invented Here, когда нужно проверять, одинаково ли мы с разработчиками понимаем задачу, а кроме того, надеяться, что библиотека будет поддерживаться на тех платформах, на которых мы планируем поддерживать нашу программу… В сложных ситуациях или в программах «на один запуск» это оправдано. Но для производных?
При этом interpreter у вас уже есть и поддерживается другой командой разработчиков.
И мы с размаху вляпываемся во все прелести Not Invented Here, когда нужно проверять, одинаково ли мы с разработчиками понимаем задачу, а кроме того, надеяться, что библиотека будет поддерживаться на тех платформах, на которых мы планируем поддерживать нашу программу… В сложных ситуациях или в программах «на один запуск» это оправдано. Но для производных?
И это нормально, разделить две независимые сущности: функцию и ее производную. Ведь фактически вы строите два параллельных (но связанных) итерационных процесса. Если вам в будущем захочется рассмотреть третью сущность, вам так и так менять код.
Вы явно не хотите услышать простого: перегрузить арифметическую операцию — хорошее решение только тогда, когда у вас первая производная. Когда нужно вычислить сразу вторую — становится сложнее. Вам код переписывать нужно тогда, когда подход с символьным пакетом сработает без изменений.
Я вам предлагаю не готовый код, а подход. И вам решать, пойти по этому маршруту или нет. В обоих случаях есть свои минусы и плюсы. Вы явно от меня ждете готовых ответов там, где я предлагаю задуматься.
Насчет NIH, извините, это аргумент не к месту. Математика — она на всех одна. И тот же GiNaC будет вычислять так, как того математика требует. Проблемы будут с поддержкой JIT генерации байт-кода — тут да, это еще область исследований, готового подхода нет. Хотя вот Numba для Python'а уже показывает результаты, сравнимые с кодом на numpy (то есть на фортране, быстрее уже вряд ли выйдет), но это пример немного не отсюда.
Вы явно не хотите услышать простого: перегрузить арифметическую операцию — хорошее решение только тогда, когда у вас первая производная. Когда нужно вычислить сразу вторую — становится сложнее. Вам код переписывать нужно тогда, когда подход с символьным пакетом сработает без изменений.
И мы с размаху вляпываемся во все прелести Not Invented Here, когда нужно проверять, одинаково ли мы с разработчиками понимаем задачу, а кроме того, надеяться, что библиотека будет поддерживаться на тех платформах, на которых мы планируем поддерживать нашу программу…
Я вам предлагаю не готовый код, а подход. И вам решать, пойти по этому маршруту или нет. В обоих случаях есть свои минусы и плюсы. Вы явно от меня ждете готовых ответов там, где я предлагаю задуматься.
Насчет NIH, извините, это аргумент не к месту. Математика — она на всех одна. И тот же GiNaC будет вычислять так, как того математика требует. Проблемы будут с поддержкой JIT генерации байт-кода — тут да, это еще область исследований, готового подхода нет. Хотя вот Numba для Python'а уже показывает результаты, сравнимые с кодом на numpy (то есть на фортране, быстрее уже вряд ли выйдет), но это пример немного не отсюда.
Более того (хотя Mrrl привел тоже хороший пример), пусть вы, например, решаете какое-то уравнение
(В этом примере производную
f(x,a)=0 относительно x тем же методом Ньютона. И хотите получать не только корень, но и его производную по a. Символьно вы это никак не сделаете, т.к. условие выхода из цикла зависит от текущего (на итерации) значения x. (Ну или будете вынуждены фиксировать количество итераций). А с Derivable все просто:// решает уравнение f(x,a)=x*x-a=0 и возвращает корень и его производную по a
Derivable sqrtWithDeriv(double a) {
static const double EPS = 1e-6;
Derivable ad = Derivable::IndependendVariable(a);
Derivable x = 1;
while (1) {
Derivable thisF = x*x - ad;
if (fabs(thisF.getValue()) < EPS) return x;
x = (ad + x * x)/2/x;
}
}
(В этом примере производную
df/dx я посчитал руками, но можно было и сделать Derivable по двум переменным и написать цикл так же, как и в примере.) посмотрите пожалуйста мой камент выше habrahabr.ru/company/intel/blog/170729/#comment_5937583
Мне кажется, вы под символьными вычислениями сейчас поняли безкомпромисное решение «все целиком считать символьно, а когда данные распухнут, вот посчитаем в точке». Нет, я лишь предлагаю построить аналог вашего класса, который производную той же одной итерации посчитает однажды, но символьно. И по результату построит быстрый байт-код. Вы получите те же два параллельных численных процесса, каждый сходящийся к своему ответу, но символьные манипуляции вы получите гарантированно правильно.
Мне кажется, вы под символьными вычислениями сейчас поняли безкомпромисное решение «все целиком считать символьно, а когда данные распухнут, вот посчитаем в точке». Нет, я лишь предлагаю построить аналог вашего класса, который производную той же одной итерации посчитает однажды, но символьно. И по результату построит быстрый байт-код. Вы получите те же два параллельных численных процесса, каждый сходящийся к своему ответу, но символьные манипуляции вы получите гарантированно правильно.
Вы, похоже, основным достоинством символьных вычисления полагаете, что для них есть готовые пакеты. Для автоматического дифференцирования готовые пакеты тоже есть, я про это неоднократно говорю в тексте и указываю, где их можно взять. И если вы думаете, что я предлагаю всегда использовать «написанный на коленке код
Но мне кажется столь же ясным, что во многих приложениях проще написать код на 50 строчек и оттестировать его, чем подключать стороннюю библиотеку. На тестирование того, что вы правильно понимаете стороннюю библиотеку и на тестирование кода связки у вас уйдет не меньше времени, чем на тестирование корректности кода
Derivable», то вы ошибаетесь. Ясно, что в определенных приложениях лучше взять готовый и более надежный код — но вот вам пожалуйста, пройдите по паре ссылок из поста и найдете их, готовых, огромное количество.Но мне кажется столь же ясным, что во многих приложениях проще написать код на 50 строчек и оттестировать его, чем подключать стороннюю библиотеку. На тестирование того, что вы правильно понимаете стороннюю библиотеку и на тестирование кода связки у вас уйдет не меньше времени, чем на тестирование корректности кода
Derivable. И если ваша задача — считать производную для функции, записанной несколькими операторами присваивания с, возможно, использованием вспомогательных функций (типичный пример в комменте), то это будет работать, и не потребует (что, на мой взгляд, весьма важно) никакого изменения в промежуточный код, кроме замены double на Derivable.в моем исходном комментарии я задавался вопросом о том, какова постановка задачи: и математически, и с точки зрения дизайна кода. Я до сих пор четкого ответа не услышал, мне три комментатора сейчас отвечают так, словно я эту постановку задачи уже слышал, но не принял к сведению.
Если задача стоит «есть много готового кода, надо один раз добавить фичу», то все что я тут пишу, разумеется, нерелевантно. Если перед вами задача периодически добавлять математический функционал, причем функционал может быть сколь угодно сложным, то предлагаемый мной класс решений по крайней мере стоит рассмотреть. Вы в посте писали про сложное моделирование, это и послужило для меня мотивацией.
Основное достоинство символьного пакета не только в том что он есть, но и в том, что он заведомо реализует продвинутый математический функционал, который вам реализовывать вручную очень дорого. Это оправдывает его использование, если есть потребность в таком функционале. А вот есть потребность или нет — опять возвращаемся к постановке задачи, которую вы явно и не описали.
Если задача стоит «есть много готового кода, надо один раз добавить фичу», то все что я тут пишу, разумеется, нерелевантно. Если перед вами задача периодически добавлять математический функционал, причем функционал может быть сколь угодно сложным, то предлагаемый мной класс решений по крайней мере стоит рассмотреть. Вы в посте писали про сложное моделирование, это и послужило для меня мотивацией.
Основное достоинство символьного пакета не только в том что он есть, но и в том, что он заведомо реализует продвинутый математический функционал, который вам реализовывать вручную очень дорого. Это оправдывает его использование, если есть потребность в таком функционале. А вот есть потребность или нет — опять возвращаемся к постановке задачи, которую вы явно и не описали.
Задача стояла так:
Я даже выше привел условный, но хорошо передающий задачу, пример: habrahabr.ru/company/intel/blog/170729/#comment_5934705.
И да, в вашем исходном комментарии вы заявляли, что для
За ссылки на JIT-компиляторы и символьные пакеты для C++ спасибо, но он не отменяют осмысленности и применимости описанного мною подхода.
Казалось бы, все просто, но для этого надо уметь вычислять производную левой части уравнения. Эта левая часть у меня была довольно сложная — даже просто вычисление ее значений в программе было разбросано по нескольким функциям,
Я даже выше привел условный, но хорошо передающий задачу, пример: habrahabr.ru/company/intel/blog/170729/#comment_5934705.
И да, в вашем исходном комментарии вы заявляли, что для
sin(x^2) не сработает, а также утверждали, что все это не нужно вообще.За ссылки на JIT-компиляторы и символьные пакеты для C++ спасибо, но он не отменяют осмысленности и применимости описанного мною подхода.
Кстати,
Мне кажется, мой код по сути именно так и сработает. Компилятор c++, компилируя строку
именно что символьно (!) построит код вычисления производной, по результаты построит быстрый исполняемый код, и получит два параллельных численных процесса, сходящиеся к своему ответу. И в том, что компилятор произведет символьные вычисления гарантированно правильно, я не сомневаюсь.
Нет, я лишь предлагаю построить аналог вашего класса, который производную той же одной итерации посчитает однажды, но символьно.
Мне кажется, мой код по сути именно так и сработает. Компилятор c++, компилируя строку
x = (ad + x * x)/2/x;
именно что символьно (!) построит код вычисления производной, по результаты построит быстрый исполняемый код, и получит два параллельных численных процесса, сходящиеся к своему ответу. И в том, что компилятор произведет символьные вычисления гарантированно правильно, я не сомневаюсь.
компиллятор выполнит это вычисление, нет проблем. Как насчет интеграла от функции?
А какой из известных пакетов способен выполнить символьное интегрирование в рамках C++ — библиотеки?
It should be noted that GiNaC is not very good (yet?) at symbolically evaluating integrals. In fact, it can only integrate polynomials.
Да, но компиллятор и этого не может, у него задача совсем другая.
Кроме того, задача символьного интегрирования в общем случае не решена (да и как его определить, общий случай), и решение для каждого класса подинтегральных выражений строится по-разному. Нужно сперва определиться с классом подинтегральных функций, затем смотреть, решена эта задача в том или ином пакете.
Кроме того, задача символьного интегрирования в общем случае не решена (да и как его определить, общий случай), и решение для каждого класса подинтегральных выражений строится по-разному. Нужно сперва определиться с классом подинтегральных функций, затем смотреть, решена эта задача в том или ином пакете.
Я еще подумал, в этом коде будут довольно сильные погрешности; например, для a=1 код вернет производную 0, а не 1/2, как должно быть. (Из-за того, что и последовательность значений сходится к нужному значению, и последовательность производных тоже, но последняя на пару итераций медленнее.)
Правда, это решается добавлением нескольких дополнительных итераций, и столь же правда, что другие методы подсчета производной, думаю, будут давать схожие проблемы по тем же в точности причинам.
Поэтому, видимо, этот код стоит считать лежащим на грани применимости автоматического дифференцирования.
Правда, это решается добавлением нескольких дополнительных итераций, и столь же правда, что другие методы подсчета производной, думаю, будут давать схожие проблемы по тем же в точности причинам.
Поэтому, видимо, этот код стоит считать лежащим на грани применимости автоматического дифференцирования.
Sign up to leave a comment.
Автоматическое дифференцирование «на пальцах»