Looping in Julia

New to Julia
1

I have written the code displayed below to calculate the maximum distance which is working perfectly. I want to modify the code for another task but facing some challenges. I will be grateful for any help and guidance please. This is what I want to do:
I want to calculate the distance from each center (xs,ys) to all points (a,b). Thus if center is (0,0) I need the distance from this center to all points in (a,b). I also want to create a matrix such that if the distance between a center to the various points is less than or equal to 2(<=2), it should write one in the matrix but if its greater than 2 it should write zero in the matrix. This process should happen from each of the centers (xs,ys) to all points in (a,b) for me to get a matrix with ones and zeros after calculating the distance from each center to all points . I am particularly interested in displaying the matrix with ones and zeros after each iteration. I will be very grateful for help to get this done please. Thank you for helping me. Below is my code which I want to modify to accomplish this task.


xs=[0 2.5 -2.5 2 -2.5 -3 -4 1.5 -1 3.7] #x-cordinates of center
ys=[0 0.3 0.25 -2.5 -2.25 1.5 -1.7 2 -3 -2]# y-cordinates of center

a=[-1, -1, 1.2, 1.2, 1.2, 0.8, 2.5, 2.5, -3, -3, -1.5, 3, 0.9, -3.7, 3.7, 3, 1.8, -1.8, -3.7, -3.7, 0.5, 2.5, -2.2, 3.5, 0.5 ] #x cordinate of the points
b=[0.2, 0.8, -0.85, 1.2, 0.5, -1.5, -1.3, 1.4, 1.5, -1.3, -0.8, -2.7, -2.4, 0.5, 0.8, 2.8, 2.4, -2.6, -2.6, 2.6, 1.5, -2.5, -3.3, -0.8, 3.1 ] #y cordinate of the points
r=[]  

function distance(xs, ys, a, b) 
    r = Float64[]
    for (i,xv) in enumerate(xs)
        R = Float64[]
        for (j, jv) in enumerate(a)
           
                push!(R, ((xs[i]-a[j])^2 + (ys[i]-b[j])^2)^0.5)
        
         end
         if length(R) > 0
             push!(r, maximum(R))
         end
    end
    return r
end
r=distance(xs, ys, a, b)
2

How about Distances.pairwise?

using Distances

xs=[0 2.5 -2.5 2 -2.5 -3 -4 1.5 -1 3.7] #x-cordinates of center
ys=[0 0.3 0.25 -2.5 -2.25 1.5 -1.7 2 -3 -2]# y-cordinates of center

a=[-1, -1, 1.2, 1.2, 1.2, 0.8, 2.5, 2.5, -3, -3, -1.5, 3, 0.9, -3.7, 3.7, 3, 1.8, -1.8, -3.7, -3.7, 0.5, 2.5, -2.2, 3.5, 0.5 ] #x cordinate of the points
b=[0.2, 0.8, -0.85, 1.2, 0.5, -1.5, -1.3, 1.4, 1.5, -1.3, -0.8, -2.7, -2.4, 0.5, 0.8, 2.8, 2.4, -2.6, -2.6, 2.6, 1.5, -2.5, -3.3, -0.8, 3.1 ] #y cordinate of the points

xy = [xs ; ys]
ab = [a' ; b']
pairwise(Euclidean(), xy, ab)
maximum(pairwise(Euclidean(), xy, ab))

julia> xy = [xs ; ys]
2×10 Matrix{Float64}:
 0.0  2.5  -2.5    2.0  …  1.5  -1.0   3.7   
 0.0  0.3   0.25  -2.5     2.0  -3.0  -2.0   

julia> ab = [a' ; b']
2×25 Matrix{Float64}:
 -1.0  -1.0   1.2   …  -2.2   3.5  0.5       
  0.2   0.8  -0.85     -3.3  -0.8  3.1       

julia> pairwise(Euclidean(), xy, ab)
10×25 Matrix{Float64}:
 1.0198   1.28062  …  3.59026  3.14006       
 3.50143  3.53553     1.48661  3.44093       
 1.50083  1.59765     6.09118  4.13793       
 4.03609  4.45982     2.26716  5.79741       
 2.87272  3.3989      6.17272  6.13372       
 2.38537  2.11896  …  6.89493  3.84838       
 3.55106  3.90512     7.55381  6.57951       
 3.08058  2.77308     3.44093  1.48661       
 3.2      3.8         5.00899  6.28172       
 5.18941  5.47083     1.21655  6.0208        

julia> maximum(pairwise(Euclidean(), xy, ab))
8.713208364316786

For reprod your result, dims = 2 is included.

julia> maximum(pairwise(Euclidean(), xy, ab), dims = 2) # my code
10×1 Matrix{Float64}:
 4.522167621838006
 6.844705983459042
 6.241193796061776
 7.648529270389178
 7.466759672039807
 7.323933369440222
 8.32165848854662
 6.942621983083913
 7.045565981523415
 8.713208364316786

julia> r=distance(xs, ys, a, b) # your code
10-element Vector{Float64}:
 4.522167621838006
 6.844705983459042
 6.241193796061776
 7.648529270389178
 7.466759672039807
 7.323933369440222
 8.32165848854662
 6.942621983083913
 7.045565981523415
 8.713208364316786
1 Like
3

There are many ways to answer the question. I will give you one answer which is very close to your approach. Just to avoid over-optimisation and making it complicated.

The new aspect here is that while for vectors you can use push!, for matices you need
to create a matrix with the correct size from the beginning!
In general, it is always a good idea to initialise your arrays directly with the size they need.

function distance(xs, ys, a, b) 
    n = length(xs)
    m = length(a)
    r = zeros(n, m)  # create the matrix and fill it with zeros
    for j in 1:m
        for i in 1:n
            if sqrt( (xs[i]-a[j])^2 + (ys[i]-b[j])^2 ) < 2 
                r[i,j] = 1
                @show i j r      # display the matrix in each iteration
            end
         end
    end
    return r
end

Now, there are many ways how to make it faster and more elegant,
for example as @rmsmsgood showed you!

4

Here’s an approach built on broadcasting. It will return a logical array, a BitMatrix, instead of a Matrix{Float64}, which may make sense, or not be quite what you are looking for:

r = @. sqrt((xs - a)^2 + (ys - b)^2) < 2

The @. is a macro which puts broadcasting dots in all the right places. It’s the same as writing:

r = sqrt.((xs .- a).^2 .+ (ys .- b).^2) .< 2

If you would rather prefer a Matrix{Bool} or a Matrix{Float64}, etc., you can just pre-allocate r and then move the @. in front of the entire expression:

r = Matrix{Bool}(undef, length(a), length(xs))  # or Matrix{Float64}(...)
@. r = sqrt((xs - a)^2 + (ys - b)^2) < 2

which is the same as

r = Matrix{Bool}(undef, length(a), length(xs))
r .= sqrt.((xs .- a).^2 .+ (ys .- b).^2) .< 2

Broadcasting can really simplify some awkward expressions.

A couple of remarks on your code:

If you know the size of the array beforehand, it is much more efficient to pre-allocate it and then fill it, than to initialize an empty vector and push!ing repeatedly.

You should normally use sqrt(x) instead of x^0.5. I’m unsure if there are differences in accuracy, but sqrt is much, much faster.

5

When I run this code I get the error “distance (generic function with 2 methods)”

How do I make this work, please?

6

I just called the function and its working now. Thank you

7

@rmsmsgood .
This method is not giving me the expected results, please. I actually want to display a matrix with ones and zeros as I gave the conditions in my question. Thank you for you guidance. I have learnt a lot with the help. The code from @ SteffenPL gives an expected results. Highly appreciated please. Thank you