Решение проблем округления с плавающей точкой C ++

Question

Я разрабатываю научное приложение (симуляция движения хромосом в ядре клетки).Хромосомы разделены на маленькие фрагменты, которые вращаются вокруг случайной оси с использованием матриц вращения 4x4.

Проблема состоит в том, что моделирование выполняет сотни миллиардов вращений, поэтому ошибки округления с плавающей точкой складываются и растут экспоненциально,поэтому фрагменты имеют тенденцию «уплывать» и отделяться от остальной части хромосомы с течением времени.

Я использую двойную точность с C ++.Софт работает на CPU на данный момент, но будет перенесен на CUDA, и симуляции могут длиться не более 1 месяца.

Я понятия не имею, как я мог бы как-то перенормировать хромосому, потому чтовсе фрагменты связаны друг с другом (вы можете видеть это как двусвязный список), но я думаю, что это будет лучшая идея, если это возможно.

У вас есть какие-либо предложения?Я чувствую себя немного потерянным.

Большое спасибо,

H.

EDIT: добавлен упрощенный пример кода.Вы можете предположить, что все математические математике являются классическими реализациями.

// Rotate 1000000 times
for (int i = 0; i < 1000000; ++i)
{
    // Pick a random section start
    int istart = rand() % chromosome->length;

    // Pick the end 20 segments further (cyclic)
    int iend = (istart + 20) % chromosome->length;

    // Build rotation axis
    Vector4 axis = chromosome->segments[istart].position - chromosome->segments[iend].position;
    axis.normalize();

    // Build rotation matrix and translation vector
    Matrix4 rotm(axis, rand() / float(RAND_MAX));
    Vector4 oldpos = chromosome->segments[istart].position;

    // Rotate each segment between istart and iend using rotm
    for (int j = (istart + 1) % chromosome->length; j != iend; ++j, j %= chromosome->length)
    {
        chromosome->segments[j].position -= oldpos;
        chromosome->segments[j].position.transform(rotm);
        chromosome->segments[j].position += oldpos;
    }
}

Answer 1

Вам нужно найти какое-то ограничение для вашей системы и работать, чтобы удержать его в разумных пределах. Я провел несколько симуляций молекулярных столкновений, и в этих системах полная энергия сохраняется, поэтому каждый шаг я дважды проверяю общую энергию системы и, если она изменяется по некоторому порогу, то я знаю, что мой временной шаг был выбран неправильно (слишком большой или слишком маленький), и я выбираю новый временной шаг и перезапускаю его. Таким образом, я могу отслеживать, что происходит с системой в режиме реального времени.

Для этой симуляции я не знаю, какое у вас есть консервативное количество, но если оно у вас есть, вы можете попытаться сохранить его постоянным. Помните, что уменьшение вашего временного шага не всегда увеличивает точность, вам нужно оптимизировать размер шага с той точностью, которую вы имеете. Я провел численное моделирование в течение нескольких недель процессорного времени, и сохраненные величины всегда были в пределах 1 части в 10 ^ 8, так что возможно, вам просто нужно поэкспериментировать с некоторыми.

Кроме того, как сказал Томалак, возможно, старайтесь всегда ссылаться на вашу систему во время запуска, а не на предыдущий шаг. Поэтому вместо того, чтобы постоянно перемещать свои хромосомы, держите хромосомы на их начальном месте и сохраняйте с ними матрицу трансформации, которая приведет вас к текущему местоположению. Когда вы вычисляете свое новое вращение, просто измените матрицу преобразования. Это может показаться глупым, но иногда это работает хорошо, потому что ошибки в среднем равны 0.

Например, допустим, у меня есть частица, которая находится в точке (x, y), и на каждом шагу я вычисляю (dx, dy) и перемещаю частицу. Пошаговый способ сделал бы это

t0 (x0,y0)
t1 (x0,y0) + (dx,dy) -> (x1, y1)
t2 (x1,y1) + (dx,dy) -> (x2, y2)
t3 (x2,y2) + (dx,dy) -> (x3, y3)
t4 (x3,30) + (dx,dy) -> (x4, y4)
...

Если вы всегда ссылаетесь на t0, вы можете сделать это

t0 (x0, y0) (0, 0)
t1 (x0, y0) (0, 0) + (dx, dy) -> (x0, y0) (dx1, dy1)
t2 (x0, y0) (dx1, dy1) + (dx, dy) -> (x0, y0) (dx2, dy2)
t3 (x0, y0) (dx2, dy2) + (dx, dy) -> (x0, y0) (dx3, dy3)

Таким образом, в любое время, tn, чтобы получить реальную позицию, вы должны сделать (x0, y0) + (dxn, dyn)

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

Answer 2

Напишите свои формулы, чтобы данные для временного шага T не были получены исключительно из данных с плавающей запятой в временном шаге T-1. Постарайтесь убедиться, что производство ошибок с плавающей запятой ограничено одним временным шагом.

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

Answer 3

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

Вариант 1:

Найдите некоторый набор ограничений, таких, что (1) они должны всегда выполняться, (2) если они терпят неудачу, но только справедливо, легко настроить систему таким образом, (3) если они все держат, то ваш симуляция не сходит с ума, и (4) когда система начинает сходить с ума, ограничения начинают выходить из строя, но лишь незначительно. Например, возможно, расстояние между соседними битами хромосомы должно быть не более d, для некоторого d, а если несколько расстояний чуть больше d, то вы можете (например) пройти по хромосоме с одного конца, фиксируя любые расстояния, которые слишком велики при перемещении следующего фрагмента к его предшественнику вместе со всеми его преемниками. Или что-то.

Затем проверяйте ограничения достаточно часто, чтобы быть уверенным, что любое нарушение все равно будет небольшим при обнаружении; и когда вы поймаете нарушение, исправьте это. (Вам, вероятно, следует договориться о том, что когда вы исправляете ситуацию, вы «более чем удовлетворяете» ограничения.)

Если проверять ограничения все время дешево, то, конечно, вы можете это сделать. (Это также может позволить вам сделать исправление дешевле, например, если это означает, что любые нарушения всегда незначительны.)

Вариант 2:

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

Вариант 3:

Вместо того, чтобы думать о своих ограничениях как об ограничениях, предоставьте несколько небольших «восстанавливающих сил», чтобы, когда что-то выходит из-под контроля, оно тянулось назад, как и должно быть. Позаботьтесь о том, чтобы, когда все в порядке, восстанавливающие силы были равны нулю или, по крайней мере, очень незначительны, чтобы они не мешали вашим результатам хуже, чем исходные числовые ошибки.

Answer 4

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

1. Вместо того, чтобы записывать позицию, поскольку некоторая начальная позиция работала много раз, вы можете записать, какой оператор будет явно после N операций. Например, представьте, что у вас есть позиция x, и вы добавляете значение e (которое вы не можете точно представить.) Гораздо лучше, чем вычисление x + = e; большое количество раз будет вычислять x + EN; где EN - более точный способ представления того, что происходит с операцией после N раз. Вам следует подумать, есть ли у вас какой-то способ более точно представить действие многих вращений.
2. Чуть более искусственно взять новую найденную точку и спроецировать любые расхождения с ожидаемым радиусом от вашего центра вращения. Это гарантирует, что оно не сместится (но не обязательно гарантирует, что угол поворота является точным.)

Answer 5

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

В зависимости от обстоятельств вам может потребоваться несколько раз применить это ограничение во время основного цикла.

Answer 6

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

Например, с точностью до 4 десятичных знаков вы получите

значение с плавающей точкой -> значение типа int 1.0000 -> 10000 1.0001 -> 10001 0,9999 -> 09999

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

1,0001 * 1,0001 становится 10001 * 10001/10000

Answer 7

Я думаю, это зависит от используемого вами компилятора.

Компилятор Visual Studio поддерживает ключ / fp, который сообщает поведение операций с плавающей запятой

, вы можете читатьоб этом .В основном, / fp: строгий режим самый жесткий