Imagine you have a series of n points randomly generated in a box in 3D space. You also have a list of distance bounds, e.g. points 5 and 3 should be between 1.0 and 2.0 Angstroms apart.There are various methods to get from random coordinates to a solution that fills some or all distance bounds. Below is an implementation of Stochastic Proximity Embedding, one such method, in Julia. I am looking for input on improving the speed.function run_spe(present_i, present_j, bounds_lower, bounds_upper, no_atoms)    squares_lower = bounds_lower .* bounds_lower    squares_upper = bounds_upper .* bounds_upper    learning_rate = 1.0    coord_new_i = zeros(3)    coord_new_j = zeros(3)    bound_near = 0.0    no_cycles = 1500    no_steps = 70000    no_bounds = length(present_i)    # Generate random starting coordinates    coords = rand!(zeros(3, no_atoms)) * 100.0    for cycle in 1:no_cycles        for step in 1:no_steps            # Pick a random bound from the list of bounds            rand_ind = rand(1:no_bounds)            i = present_i[rand_ind]            j = present_j[rand_ind]            # Compute square distance between atoms            sq_dist = 0.0            for m = 1:3                sq_dist += abs2(coords[m, i] - coords[m, j])            end            # Check if square distance is within bounds            # If it is, go to the next iteration            # If it is not, find the nearest bound            if sq_dist < squares_lower[rand_ind]                bound_near = bounds_lower[rand_ind]            elseif sq_dist > squares_upper[rand_ind]                bound_near = bounds_upper[rand_ind]            else                continue            end            # Adjust coordinates to bring them closer to the nearest bound            distance = sqrt(sq_dist)            prefactor = learning_rate * 0.5 * (bound_near - distance) / distance            for k in 1:3                coord_new_i[k] = coords[k, i] + prefactor * (coords[k, i] - coords[k, j])                coord_new_j[k] = coords[k, j] + prefactor * (coords[k, j] - coords[k, i])            end            for l in 1:3                coords[l, i] = coord_new_i[l]                coords[l, j] = coord_new_j[l]            end        end        # Decrease learning rate after each cycle        learning_rate -= 1 / no_cycles    end    return coordsendpresent_i/present_j are lists of indices. bounds_lower/bounds_upper are lists of distances. For example, present_i[10] = 5, present_j[10] = 3, bounds_lower[10] = 1.0, bounds_upper[10] = 2.0 would mean the tenth bound corresponds to the example bound above.I have done a bit of profiling and the major bottlenecks are the initial lines in the inner loop (from rand_ind up to continue). Ideally Julia tricks to speed individual steps up would be helpful, but if there is a better data structure than the lists I use that would also be good.Sample data can be found here for testing and read in using something quick like:# Set in_filepath as path to downloaded and uncompressed filepresent_i = Int[]present_j = Int[]bounds_lower = Float64[]bounds_upper = Float64[]no_atoms = 2035open(in_filepath, r) do in_file    for line in eachline(in_file)        cols = split(strip(line), '\t')        push!(present_i, int(cols[1]))        push!(present_j, int(cols[2]))        push!(bounds_lower, float(cols[3]))        push!(bounds_upper, float(cols[4]))    endend
Iterative procedure to get 3D coordinates from distance constraints
performance;iteration;numerical methods;julia