# For Loops giving weird results

Julia is giving me weird results when I run the following code.

``````using Plots
plotly()
using PyCall
@pyimport numpy as np

"""function for calculating the corner frequency of the ijth subfault"""
function corner_f0ij(M0j::Float64,Nrj::Int64,N::Int64,β::Float64,σ::Float64)
σ = σ/100000 #input stress drop in pascal output in bars
β = β/1000 # input shear wave velocity in m/sec output in km/sec
#M0 in dynecm
return 4.9*10^6*Nrj^(-1/3)*N^(1/3)*β*(σ/M0j)^(1/3)
end

"""function for calculating the FAS of the source"""
function source(f::Float64,f0j::Float64,M0j::Float64) #input M0j in Nm
return (2*π*f)^2*M0j/(1+(f/f0j)^2)
end

"""function to account the geometrical spreading"""
function path(f::Float64,r::Float64,β::Float64)
if r <= 100000.
G = 1/r
else
G = 1/(10*sqrt(r))
end
Q = 333*f^0.74
return G*exp(-π*f*r/(β*Q))
end

"""function to account the constant"""
function cons(rad::Float64,ρ::Float64,β::Float64)
return rad*sqrt(2)/(4*π*ρ*β^3)
end

"""function for accounting the filter"""
function filter_Pj(Hj::Float64,N::Int64,f::Float64,f0::Float64,f0j::Float64)
return sqrt(N)/Hj*(1+(f/f0j)^2)/(1+(f/f0)^2)
end

"""function to account the local site(amplification or demaplification)"""
function local_site(f::Float64,k0::Float64)
x = [0.5,1.,2.,5.,10.]
y = [1,1.13,1.22,1.36,1.41]
amp = np.interp(f,x,y)
return amp*exp(-π*k0*f)
end

"""function for calculating the distance from sub-fault to observation point \n
Input these parameters
Δl = subfault length
Δw = subfault width
ϕ1 = fault strike
ϕ2 = azimuth to the observation point
δ1 = fault dip
i,j = subfault number
R = epicentral distance
"""
function r1(Δl::Float64,Δw::Float64,ϕ1::Float64,ϕ2::Float64,δ1::Float64,
h::Float64,i::Int64,j::Int64,R::Float64)
δ = 90.-δ1
return ((R*cosd(ϕ2-ϕ1)-(2*i-1)*Δl/2)^2+(R*sind(ϕ2-ϕ1)-(2*j-1)*Δw/2*sind(δ))^2 +
(h+(2*j-1)*Δw/2*cosd(δ))^2)^0.5
end

"""function for calculating the number of pulsing subfaults at a time ij \n
i0,j0 = location of the hypocenter
i,j = (i,j) the subfault at the ijth time
nw = number of fault rows
nl = number of fault columns
n0_of_effective_subfaults = number of faults once the ruptre propagates
"""
function pulsing_subfaults(i0::Int64,j0::Int64,i::Int64,j::Int64, nw::Int64,
nl::Int64,no_of_effective_subfaults::Float64)

Rmax = maximum([abs(i-i0)+1,abs(j-j0)+1])
Rmin = Rmax-no_of_effective_subfaults

if Rmin<0
Rmin = 0
end

n = 0
for jj in 1:nl
for ii in 1:nw
r = maximum([abs(ii-i0)+1,abs(jj-j0)+1])
if r>Rmin && r<Rmax
n = n+1
end
end
end
return n
end

### Now the problem starts in here inside this function

function FAS(M::Float64,R::Float64,dij::Array{Float64,2},l::Float64,w::Float64,i0::Int64,j0::Int64,
pulsing_percent::Float64,β::Float64,ρ::Float64,σ::Float64,rad::Float64,k0::Float64,
ϕ1::Float64,ϕ2::Float64,δ1::Float64,h::Float64)

M0 = 10^((M+10.7)*(3./2.))*1.0e-7 #Nm
nw = size(dij)[1]
nl = size(dij)[2]
N = length(dij)
Δl = l/nl
Δw = w/nw
subfault_radius = sqrt(Δl*Δw/π)
no_effective_subfaults = nl*pulsing_percent/100
no_effective_subfaults = no_effective_subfaults/2
vrup = 0.8*β

if no_effective_subfaults < 1.
no_effective_subfaults = 1.
end

#declaration of variables
dij_sum = sum(dij)
t_end_max = 0.
t-arrive_min = 10000.
risetime_max = 0.0

#initialization of variables
no_active_subfaults = zeros(Int64,nw,nl)
f0ij = zeros(Float64,nw,nl)
M0ij = zeros(Float64,nw,nl)
risetimeij = zeros(Float64,nw,nl)
subfault_distance = zeros(Float64,nw,nl)
dur_subij = zeros(Float64,nw,nl)
delay = zeros(Float64,nw,nl)

for j in 1:nl
for i in 1:nw
delay[i,j] = i-i0 ## this is the problem if I type i0 or any values say 1 or 2 it gives me 1000
M0ij[i,j] = M0*dij[i,j]/dij_sum
no_active_subfaults[i,j] =  pulsing_subfaults(i0,j0,i,j,nl,nw,no_effective_subfaults)
if no_active_subfaults[i,j] == 0
no_active_subfaults[i,j] = 1
end
f0ij[i,j] = corner_f0ij(M0ij[i,j]*1.0e7,no_active_subfaults[i,j],N,β,σ)
subfault_distance[i,j] = r1(Δl,Δw,ϕ1,ϕ2,δ1,h,i,j,R)
risetimeij[i,j] = 1.0/f0ij[i,j]
if risetimeij[i,j] > risetime_max
risetime_max = risetimeij[i,j]
end

dur_subij[i,j] = risetimeij[i,j] + 0.05*subfault_distance[i,j]/1000.;

end
end
return delay
end

function main()
location="F:\\Books\\Thesis\\Paper\\Scripts\\displacements.txt"
dij = readdlm(location,',')
M=7.9
R = 60000.
l = 195000.
w = 150000.
i0 = 6
j0 = 6
pulsing_percent = 50.0
β = 3500.
ρ = 28000.
σ = 50*100000.
rad = 0.55
k0 = 0.05
ϕ1 = 293.
ϕ2 = 120.
h = 15000.
δ1 = 7.0
value = FAS(M,R,dij,l,w,i0,j0,pulsing_percent,β,ρ,σ,rad,k0,ϕ1,ϕ2,δ1,h)
return value
end

main()
``````

results

``````10×13 Array{Float64,2}:
10000.0  10000.0  10000.0  10000.0  …  10000.0  10000.0  10000.0  10000.0
10000.0  10000.0  10000.0  10000.0     10000.0  10000.0  10000.0  10000.0
10000.0  10000.0  10000.0  10000.0     10000.0  10000.0  10000.0  10000.0
10000.0  10000.0  10000.0  10000.0     10000.0  10000.0  10000.0  10000.0
10000.0  10000.0  10000.0  10000.0     10000.0  10000.0  10000.0  10000.0
10000.0  10000.0  10000.0  10000.0  …  10000.0  10000.0  10000.0  10000.0
10000.0  10000.0  10000.0  10000.0     10000.0  10000.0  10000.0  10000.0
10000.0  10000.0  10000.0  10000.0     10000.0  10000.0  10000.0  10000.0
10000.0  10000.0  10000.0  10000.0     10000.0  10000.0  10000.0  10000.0
10000.0  10000.0  10000.0  10000.0     10000.0  10000.0  10000.0  10000.0

``````

However if I do like this

``````delay[i,j] = i+(-1*i0)
``````

then the result is fine like below

``````10×13 Array{Float64,2}:
-5.0  -5.0  -5.0  -5.0  -5.0  -5.0  -5.0  -5.0  -5.0  -5.0  -5.0  -5.0  -5.0
-4.0  -4.0  -4.0  -4.0  -4.0  -4.0  -4.0  -4.0  -4.0  -4.0  -4.0  -4.0  -4.0
-3.0  -3.0  -3.0  -3.0  -3.0  -3.0  -3.0  -3.0  -3.0  -3.0  -3.0  -3.0  -3.0
-2.0  -2.0  -2.0  -2.0  -2.0  -2.0  -2.0  -2.0  -2.0  -2.0  -2.0  -2.0  -2.0
-1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0
1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0
2.0   2.0   2.0   2.0   2.0   2.0   2.0   2.0   2.0   2.0   2.0   2.0   2.0
3.0   3.0   3.0   3.0   3.0   3.0   3.0   3.0   3.0   3.0   3.0   3.0   3.0
4.0   4.0   4.0   4.0   4.0   4.0   4.0   4.0   4.0   4.0   4.0   4.0   4.0
``````

Can some one please figure out what is going on??

Can you remove all code that is not needed to reproduce the problem?

1 Like

The original code contains these arguments

function FAS(M::Float64,R::Float64,dij::Array{Float64,2},l::Float64,w::Float64,i0::Int64,j0::Int64,
pulsing_percent::Float64,β::Float64,ρ::Float64,σ::Float64,rad::Float64,k0::Float64,
ϕ1::Float64,ϕ2::Float64,δ1::Float64,h::Float64)

``````M0 = 10^((M+10.7)*(3./2.))*1.0e-7 #Nm
nw = size(dij)[1]
nl = size(dij)[2]
N = length(dij)
Δl = l/nl
Δw = w/nw
subfault_radius = sqrt(Δl*Δw/π)
no_effective_subfaults = nl*pulsing_percent/100
no_effective_subfaults = no_effective_subfaults/2
vrup = 0.8*β

if no_effective_subfaults < 1.
no_effective_subfaults = 1.
end

#declaration of variables
dij_sum = sum(dij)
t_end_max = 0.
t-arrive_min = 10000.
risetime_max = 0.0

#initialization of variables

delay = zeros(Float64,nw,nl)

for j in 1:nl
for i in 1:nw

delay[i,j] = i-i0  #the problem is here, if i do

end
end
return delay
``````

end

function main()
location=“F:\Books\Thesis\Paper\Scripts\displacements.txt”
dij = readdlm(location,‘,’)
M=7.9
R = 60000.
l = 195000.
w = 150000.
i0 = 6
j0 = 6
pulsing_percent = 50.0
β = 3500.
ρ = 28000.
σ = 50*100000.
rad = 0.55
k0 = 0.05
ϕ1 = 293.
ϕ2 = 120.
h = 15000.
δ1 = 7.0
value = FAS(M,R,dij,l,w,i0,j0,pulsing_percent,β,ρ,σ,rad,k0,ϕ1,ϕ2,δ1,h)
return value
end

main()

And removing anything of that code (for example the declaration of all the unused variables) makes the problem go away?

Thanks man finally I was able to resolve the issue. But I think Julia could’ve reported error

instead of setting t_arrive_min = 100000.
so inside the nested for loops whenever “-” operator occurred the value was set to 10000.

It should have reported error instead of giving this weird result. What is the reason behind this?

You’ll want to make sure the answer gets put over there as well.

Quite obscure bug. Some people have argued that this type of function declarations should be removed.

1 Like

Thanks, yeah I’ve posted there as well.

This issue with Julia comes up from time to time on this discourse site, and each time it comes up, I complain again. The issue is: it is too easy to accidentally redefine an operator in Base via a tiny typo, leaving a program with a bug extremely difficult to locate.

I was also bitten by it a few years ago and opened another issue. The core developers have responded to my posts that they have no intention of fixing this issue because it is part of the design of the language. But as far as I know, there is not a single example of anyone ever intentionally redefining an Base operator inside a function using this syntax. In other words, this language feature has been used so far only to create bugs.