LINUX.ORG.RU

иключить деление на ноль

 


0

1

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

$ cat ttest.py 
#!/usr/bin/python
import numpy 
import sympy
x,y=sympy.symbols('x y')
G=1./(x-y)
F=sympy.diff(G,x)
sympy.pprint(F)
def FunWithoutDivZero(sFunct=None):
    def lmbd(a,b):
        result=numpy.zeros_like(a-b)
        gd=(a!=b)
        #a[numpy.logical_not(gd)]+=1.
        result[gd]=sympy.lambdify((x,y),sFunct,modules='numpy')(a,b)[gd]
        return result
    return lmbd 
KF=FunWithoutDivZero(sFunct=F)
t=numpy.array([1.,2.,3.,4.,5.])
u=numpy.array([1.,2.,3.,4.,5.]).reshape(-1,1)
print KF(t,u)
$ ./ttest.py 
   -1   
────────
       2
(x - y) 
Warning: divide by zero encountered in true_divide
[[ 0.         -1.         -0.25       -0.11111111 -0.0625    ]
 [-1.          0.         -1.         -0.25       -0.11111111]
 [-0.25       -1.          0.         -1.         -0.25      ]
 [-0.11111111 -0.25       -1.          0.         -1.        ]
 [-0.0625     -0.11111111 -0.25       -1.          0.        ]]

p.s. если бы можно было определить символьную функцию как G= при x!=y вернуть 1/(x-y), а при x=y вернуть 0.0, то все было бы ok.

p.s. 2 поэлементное решение слишком медленное (работает в десятки раз медлее):

def withoutZero(sFunct=None,eps=0.0001):
    def lmbd(a,b):
        if numpy.abs(a-b)>eps: 
            return sympy.lambdify((x,y),sFunct,modules='numpy')(a,b) 
        else: return 0.0
    return numpy.vectorize(lmbd)

1. а в задаче IRL тоже такие ограничения или то артефакты численного решения, может имеет смысл просто «переформулировать» задачу?

2. если идти не включая моск - раздели задачу на 2 этапа, на первом обработай все ситуации x=y (только не делай так: if numpy.abs(a-b)>eps), на втором обработай все ситуации требующие вычисления

ну и скастую ка я пожалуй AIv :)

shty ★★★★★
()
Ответ на: комментарий от shty

А я с sympy и numpy не работал никогда, что вообще тут происходит то? А чем не устраивает inf в качестве результата?

AIv ★★★★★
()
Ответ на: комментарий от AIv

А чем не устраивает inf в качестве результата?

как минимум консоль засоряется сообщениями

Warning: invalid value encountered in true_divide
Warning: invalid value encountered in true_divide
Warning: invalid value encountered in true_divide
....

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

dmitry-d67
() автор топика
Ответ на: комментарий от shty

1.

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

2.

разделить задачу можно (хоть и некрасиво), например в таком случае, когда x и y вектора строки, в этом случае подправляем числа, приводящие к делению на ноль:

gd=(a!=b)
a[numpy.logical_not(gd)]+=1.

Но в примере выше a-строка, b-вектор, а результат - матрица ... ;

наверное придется подавать матрицы на вход ...

dmitry-d67
() автор топика
Ответ на: комментарий от dmitry-d67

пока такое, не очень красивое решение

def FunWithoutDivZero(sFunct=None):
    def lmbd(a,b):
        A,B=numpy.meshgrid(a,b)
        result=numpy.zeros_like(A)
        gd=(A!=B)
        A[numpy.logical_not(gd)]+=1.
        result[gd]=sympy.lambdify((x,y),sFunct,modules='numpy')(A,B)[gd]
        return result
    return lmbd 

может можно как-то вообще отказаться от вычислений функции в точках с плохими индексами numpy.logical_not(gd)?

dmitry-d67
() автор топика
Ответ на: комментарий от dmitry-d67

яхз, я бы вообще сделал так - по известной формуле (пусть она получена с помощью sympy) сгенерировал бы сишный файл с расчётами, который умеет читать ветро параметров из stdin, законпелял бы его при помощи gcc и запустил новый процесс, в stdin которому бы подсунул вектор значений, будет самое быстрое решение имхо

shty ★★★★★
()
Ответ на: комментарий от shty

Категорически поддерживаю! Слова о том, что дескать хочется увернуться от лишних вычислений при x==y дабы ускорится выглядят весьма наивно - ускорение тут будет какие то проценты, а вот переход на С даст реальное ускорение. А sympy умеет кодогенерацию?

Я, когда хочу увернуться от деления на ноль, делаю в общем как ТС - т.е. проверяю abs(x-y)<e где e чиселка определяемая спецификой задачи. Обычно это происходит при вычислении ф-ии которая при x-y==0 дает ноль (если ее скажем пролапиталить), но при малых x-y стреляет ошибка округления (одно маленькое число делим на другое маленькое число получаем произвольное и совсем не маленькое число).

Если x и y формируются руками - самое правильное написать цикл так, что бы ситуации x==y просто не возникали. Еще вариант - юзать тернарный оператор ? или его аналоги, в гнуплоте это обычно выглядит так

plot x-y==0 ? 0 : ...

т.е. поверх полученной из sympy ф-ии навернуть такую проверку - наверное sympy позволяет?

AIv ★★★★★
()
Ответ на: комментарий от AIv

поверх полученной из sympy ф-ии навернуть такую проверку - наверное sympy позволяет?

вот этого я в симпай не нашел ...

dmitry-d67
() автор топика
Ответ на: комментарий от AIv

Категорически поддерживаю! Слова о том, что дескать хочется уернуться от лишних вычислений при x==y дабы ускорится выглядят весьма наивно - ускорение тут будет какие то проценты, а вот переход на С даст реальное ускорение.

numpy по скорости не уступает чистому Си, если удается решать задачу в матричном виде. В последнем примере я говорил о некрасивом решении (а код очень даже быстр :-)) и хотел бы вывернутся в рамках numpy или sympy.

В сложных случаях, когда в матричном виде выражение совсем запутанное, я пишу функции на Си (+ ctypes в Питоне)

dmitry-d67
() автор топика
Ответ на: комментарий от dmitry-d67

Это не совсем так, точнее совсем не так, потому что «чистых С» бывает ну очень много разных (я говорю про разные реализации одного и того же алгоритма). Даже в самых простых случаях скорость разных реализаций может различаться на порядки.

AIv ★★★★★
()
Ответ на: комментарий от true_admin

Насколько я понимаю, sympy раскручивает выражение в AST (я сам так делаю;-)), так что нужна некая симпаевская конструкция этому тернарному оператору отвечающая.

AIv ★★★★★
()
Ответ на: комментарий от AIv

согласен, но ничего подходящего не могу найти/придумать.

P.S. тернарный оператор отрабатывает как оператор Питон

>>> F=1/(x-y) if x!=y else 0.0
>>> sympy.pprint(F)
  1  
─────
x - y
>>> sympy.pprint(F.subs(x,y))
∞

dmitry-d67
() автор топика
Ответ на: комментарий от dmitry-d67

ес-но он возвращает 1/(x-y) поскольку x!=y (это ж символьные симпаевские экземпляры классов). Да и как иначе - тернарный оператор не перегружается... У симпая должна быть или ф-я с тремя аргументами, или что то в этом роде, реализующая символьный тернарный оператор. Если нет - я буду сильно удивлен.

AIv ★★★★★
()
Ответ на: комментарий от AIv

нашел вот, но что-то не работает

>>> from sympy import *
>>> from sympy.abc import x,y
>>> f=1./(x-y)
>>> F=Piecewise((0,x==y),(f,True))
>>> F.subs(x,y)
oo
>>> F.subs({x:1,y:1})
-oo

даже совсем в простом случае

>>> F=Piecewise((1/x,x!=0),(0,True))
>>> F.subs(x,0)
oo

P.S. http://nullege.com/codes/search/sympy.Piecewise

dmitry-d67
() автор топика
Ответ на: комментарий от dmitry-d67

к сожалению :-)

$ cat ttests.py 
#!/usr/bin/python
import numpy 
import sympy
x,y=sympy.symbols('x y')
G=1./(x-y)
F=sympy.Piecewise((0.0,sympy.abs(x-y)<0.0001), (sympy.diff(G,x),True))
sympy.pprint(F)
print F.subs({x:1,y:1})
print F.subs({x:3,y:1})
def FunWithoutDivZero(sFunct=None):
    def lmbd(a,b):
        result=sympy.lambdify((x,y),sFunct,modules='numpy')(a,b)
        return result
    return lmbd 
KF=FunWithoutDivZero(sFunct=F)
t=numpy.array([1.,2.,3.,4.,5.])
u=numpy.array([1.,2.,3.,4.,5.]).reshape(-1,1)
print KF(t,u)

$ ./ttests.py 
⎧   0      for │x - y│ < 0.0001
⎪                              
⎪   -1                         
⎨────────       otherwise      
⎪       2                      
⎪(x - y)                       
⎩                              
0
-1/4
Traceback (most recent call last):
  File "./ttests.py", line 18, in <module>
    print KF(t,u)
  File "./ttests.py", line 12, in lmbd
    result=sympy.lambdify((x,y),sFunct,modules='numpy')(a,b)
  File "<string>", line 1, in <lambda>
NameError: global name 'Piecewise' is not defined

dmitry-d67
() автор топика
Ответ на: комментарий от dmitry-d67

дело пошло, просто сравнение похоже численное

А Вы какого ожидали???

AIv ★★★★★
()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.