Fortran code adaptation for python [RuntimeWarning: overflow encountered in double scalars]
I tried adapting the following code in Fortran to python:
PROGRAM Escape_tent
IMPLICIT NONE
! Declare local variables
INTEGER :: i,j,tfinal,P,k,rede
Real,dimension(0:2800) :: x,xn
real, external :: f
Real :: sigma,pi,soma,r
pi=4.*atan(1.0)
!
! Parametros e CIs
rede=100
do j=0,75
x(j)=0.3+0.1*sin((j*2*pi)/rede)
enddo
do j=50,rede
x(j)=1.1+0.1*sin((j*2*pi)/rede+pi)
enddo
tfinal=100
sigma=0.38
r=0.32
P=int(rede*r)
!
! Mapa logistico
!
do i=1,tfinal
!atualiza os vetores com o mapa
do j=0,rede
x(j)=f(x(j))
enddo
!
!
! Finite Range Coupling
!
do j=0,rede
soma=0.0
do k=-P,P
soma=soma+x(modulo(j+K,rede+1))-x(j)
enddo
xn(j)=x(j)+(sigma/(2.*P))*soma
enddo
do j = 0,rede
x(j)=xn(j)
enddo
! Write out values of the map
write(1,100) (x(k),k=0,rede)
enddo
END PROGRAM Escape_tent
real function f(x)
implicit none
real, intent(in) :: x
real, parameter :: r=2.78
f=r*x(1 - x)
end function f
When I adapted it for python it looked like this:
import numpy as np
from time import time
# Variáveis do modelo
x, xn = [], []
r = 0.32 # 'raio' de acoplamento r = P/N
pi = np.pi
soma = 0.0
# Parâmetros para condições iniciais
rede = 100 # tamanho da rede
P = int(rede * r)
sigma = 0.38
t = np.arange(0, 1000) # tempo de iteração
#Condições iniciais da rede
for j in range(0, int(rede / 2)):
# determina as condições iniciais de metade da rede
x.insert(j, 0.3 + 0.1 * np.sin((j * 2 * pi) / (rede)))
for j in range(int(rede / 2), rede):
# determina as condições iniciais da outra metade da rede
x.insert(j, 1.1 + 0.1 * np.sin((j * 2 * pi) / (rede + pi)))
# função usada na iteração
def f(x):
return 3.8 * x * (1 - x)
t1 = time()
#iterando a rede
for i in t:
for j in range(0, rede):
x[j] = f(x[j]) # aplica a função em todos os elementos da rede
for j in range(0, rede): # para cada elemento j da rede
for k in range(0, rede): # temos um elemento k que esta a uma distância do j
soma = 0.0
if k != j: # k não interage com k
# calcula a distância do k até o j
soma = x[k] - x[j]
xn.insert(j, x[j] + (sigma / 2 * P) * soma) # add no vetor x
for j in range(0, rede):
x[j] = xn[j] % (1.0) #Aqui eu limito o valor de Xn[j] para não tender ao infinito.
print("T {}\n X :{}".format(i,x))
t2 = time()
print("Tempo decorrido: {}".format(t2 - t1))
After the ninth iteration the exception "RuntimeWarning: overflow encountered in double_scalars" and "RuntimeWarning: invalid value encountered in double_scalars" occurs in the PyCharm terminal and the data added to the X vector is 'nan'. I would like to know what is happening, because in fortran it normally runs on the python I have this problem. Thanks.
1 answers
Not a number (nan)
Your code is generating nan in the list you name from x (I suggest you watch this video), and what is a nan? In a very objective and sensible way, nan is the acronym for "Not a number" , that is, any mathematical operation that results in a value that cannot be expressed in a number, the representation of this result will be nan, usually an exception is generated when trying to perform an operation that results in a nan, for example example if you try to divide by zero, the exception ZeroDivisionError will be raised.
a = 5
b = 0
c = a/b
ZeroDivisionError: division by zero
But not always the exception is raised, unless you deal with this, it is difficult to purposely create a nan, but there is an operation with inf that is classic for an example:
a = float('Inf')
b = float('Inf')
c = a - b
print(c)
nan
What happens is that your code in python (I will not compare with the code in Fortran, because I do not know the language), is generating values nan in the list x, so when you send the value to the function f() the value nan is returned, because any operation involving a nan will result in a nan, the message RuntimeWarning: overflow encountered in double_scalars most of the time is related to operations involving values nan, find out because these values are being generated and you will solve the problem, I will make some suggestions for you to debug, but you can]}
1) change the function f() to intercept the value nan in x, as follows shape:
# função usada na iteração
def f(x):
result = 3.8 * x * (1 - x)
if np.isnan(result):
line = '='*22
print(line,'x recebido na função f: ', x,line, sep='\n')
return result
2) comment the output of x in the iteration loop, just so that the debug output is cleaner, also before the procedure that calls the function f() check if the element j in x that you are sending to the function f is not nan, if it is present the variables involved, the iteration loop should look like this:
t1 = time()
#iterando a rede
for i in t:
for j in range(0, rede):
if np.isnan(x[j]):
print('valor de j na iteracao, quando x[j] é nan', j)
print('valor de x na iteracao, quando x[j] é nan', x)
x[j] = f(x[j]) # aplica a função em todos os elementos da rede
for j in range(0, rede): # para cada elemento j da rede
for k in range(0, rede): # temos um elemento k que esta a uma distância do j
soma = 0.0
if k != j: # k não interage com k
# calcula a distância do k até o j
soma = x[k] - x[j]
xn.insert(j, x[j] + (sigma / 2 * P) * soma) # add no vetor x
for j in range(0, rede):
x[j] = xn[j]
# print("T {}\n X :{}".format(i,x))
t2 = time()
print("Tempo decorrido: {}".format(t2 - t1))
After these changes the result should be something close to:
valor de j na iteracao, quando x[j] é nan 0
valor de x na iteracao, quando x[j] é nan [nan, nan, nan,... 1.0350165585196508e+289,... ]
======================
x recebido na função f:
nan
======================
...
The output above is only one fragment, run the code with the changes to see the actual output, you can (and should) use print in the variables involved in the calculations to generate the elements in x to find out what is generating the values nan.