For Loops giving weird results

General Usage
1

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??

2

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

1 Like
3

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()

4

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

5

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?

6

Cross posted: multidimensional array - Why is the for loop giving me weird results in julia? - Stack Overflow

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

7

This is Infix operator definition syntax needs documentation · Issue #15483 · JuliaLang/julia · GitHub I believe.

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

1 Like
8

Thanks, yeah I’ve posted there as well.

9

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.