Я переношу процедуру вертикальной интерполяции из Fortran 77 в NumPy. Я стремлюсь к побитовой четности (средняя абсолютная ошибка (\около 10^{-16})) с выводом gfortran, который использует ДВОЙНУЮ ТОЧНОСТЬ. В моей реализации NumPy я сталкиваюсь с предупреждением RuntimeWarning: деление на ноль, обнаруженным в журнале. Я понимаю, что это связано с тем, что np.where быстро оценивает обе ветви x и y, попадая в заполнители 0,0 или -9999,0 в моих «неактивных» ячейках сетки перед применением маски.
Исходная логика F77
IF (ABS (PSFC (I, J) - PRES (K_IN)) .lt. 0.01) then
PMID = PRES (K_IN)
PUP = PRES (K_IN + 1)
SMID = SPRES (I, J, K_IN)
SUP = SPRES (I, J, K_IN + 1)
LNP1P2 = LNPU1P(K_IN)
LNP1P3 = LOG (PUP / PDWN)
LNP2P3 = LOG (PMID / PDWN)
Моя текущая попытка NumPy (вызывает предупреждение):
l13[active] = np.where(c1_3d, np.log(pup/pdwn), lnpu2p[k])[active]
Мое ограничение:
Я должен поддерживать (10^{-16}) MAE. Простое использование np.seterr(divide='ignore') неприемлемо. Я подозреваю, что вычисление журнала ненужных данных - даже если они отброшены - может незаметно повлиять на флаги состояния FPU или промежуточное округление в векторизованной среде, что может отклонить результаты от ссылки F77.
Вопрос:
Является ли логическое индексирование (показанное ниже) лучшим способом имитации поведения блока IF Фортрана с точки зрения точности и предотвращения предупреждений времени выполнения? Есть ли более «NumPy-tonic» способ добиться этого, более эффективный для больших 3D-сеток?
Из 168192 баллов я получаю только 4 предупреждения, что доказывает, что большая часть перезаписи
успеха.
python
mask_true = active & c1_3d
if np.any(mask_true):
l13[mask_true] = np.log(pup[mask_true] / pdwn[mask_true])
mask_false = active & (~c1_3d)
if np.any(mask_false):
l13[mask_false] = lnpu2p[k]
Подробнее здесь: https://stackoverflow.com/questions/798 ... arity-with
