Numpy to Julia conversion questions

As an exercise I’m converting this Python Differential Evolution tutorial code to Julia (Differential Evolution from Scratch in Python)

Here’s enough of the python to get going:

 from numpy.random import rand
 from numpy.random import choice
 from numpy import asarray
 from numpy import clip
 from numpy import argmin
 from numpy import min
 from numpy import around
 
 # define population size
 pop_size = 10
 # define lower and upper bounds for every dimension
 bounds = asarray([(-5.0, 5.0), (-5.0, 5.0)])
 # define number of iterations
 iter = 100
# define scale factor for mutation
 F = 0.5
 # define crossover rate for recombination
 cr = 0.7
 
 # define objective function
 def obj(x):
     return x[0]**2.0 + x[1]**2.0
 
 # define mutation operation
 def mutation(x, F):
    return x[0] + F * (x[1] - x[2])
 
 
# define boundary check operation
 def check_bounds(mutated, bounds):
     mutated_bound = [clip(mutated[i], bounds[i, 0], bounds[i, 1]) for i in range(len(bounds))]
     return mutated_bound
 
 # define crossover operation
 def crossover(mutated, target, dims, cr):
     # generate a uniform random value for every dimension
     p = rand(dims)
     # generate trial vector by binomial crossover
     trial = [mutated[i] if p[i] < cr else target[i] for i in range(dims)]
     return trial
 
 def differential_evolution(pop_size, bounds, iter, F, cr):
     # initialise population of candidate solutions randomly within the specified bounds
     pop = bounds[:, 0] + (rand(pop_size, len(bounds)) * (bounds[:, 1] - bounds[:, 0]))
     # evaluate initial population of candidate solutions
     obj_all = [obj(ind) for ind in pop]
     # find the best performing vector of initial population
     best_vector = pop[argmin(obj_all)]
     best_obj = min(obj_all)
     prev_obj = best_obj
 

Let’s focus in on the pop assignment there in the differential_evolution function trying that out in Julia, breaking it up a bit to make it easier to examine what’s happening:

pop_size = 10
bounds = [ [-5.0 5.0]
           [-5.0 5.0]
         ]
iter = 100
F    = 0.5
cr   = 0.7

Using the REPL:

julia> diff = bounds[:,2] - bounds[:,1]
2-element Vector{Float64}:
 10.0
 10.0
julia> r = rand(pop_size, size(bounds)[1])
10×2 Matrix{Float64}:
 0.488923  0.0440984
 0.662779  0.82573
 0.846927  0.536482
 ...
 0.466651  0.657436
 0.236172  0.799999
julia> m = r * diff
10-element Vector{Float64}:
  5.330218765700669
 14.885082264670109
 13.83409268389118
 ...
 11.240868248935161
 10.361709221756865

#now putting it all together:
julia> pop = bounds[:, 1] + m
ERROR: DimensionMismatch("dimensions must match: a has dims (Base.OneTo(2),), b has dims (Base.OneTo(10),), mismatch at 1")

Ok, so numpy is doing something different here - IIRC numpy will figure out if the dimensions mismatch and then use broadcasting. So what I finally came up with to get all the dimensions to match up was this:

julia> pop = (bounds[:, 1] .+ (rand(size(bounds)[1], pop_size)) .* (bounds[:, 2] - bounds[:,1]))'
10×2 adjoint(::Matrix{Float64}) with eltype Float64:
  4.72087   -3.63459
  4.31493    2.73267
...
 -0.169094  -2.07064
 -4.84554    4.71972
julia> typeof(pop)
LinearAlgebra.Adjoint{Float64, Matrix{Float64}}

Q1: Why is the type of pop LinearAlgebra.Adjoint?

Ok, I think that does what I want, now going on:

julia> obj_all = [obj(ind) for ind in pop]
10×2 Matrix{Float64}:
 22.2866     13.2102
 18.6186      7.4675
...
  0.0285929   4.28755
 23.4793     22.2758
julia> best_vector = pop[argmin(obj_all)]
-2.5181339556744553

This was just directly what the Python/numpy code was, but here I was expecting to get a vector with 2 floats, but only got one. For example, in the Python REPL:

>>> best_vector = pop[argmin(obj_all)]
>>> best_vector
array([-1.62902733, -0.44214267])

Looking at the output of that argmin back in the Julia REPL:

julia> argmin(obj_all)
CartesianIndex(9, 1)

So it’s returning a CartesianIndex(9,1) which seems to be why I’m only getting a scalar value back. Looks like I need to pull out the first item of that index to get what I want:

julia> best_vector = pop[argmin(obj_all)[1],:]
2-element Vector{Float64}:
 -2.5181339556744553
  3.2965151490022073
julia> best_obj = minimum(obj_all)
0.028592946239433044

Ok, that’s better, best_vector is actually a vector now.

Q2: Am I doing this right? It seems a bit painful to have to do things in both the Python REPL and the Julia REPL to figure out the differences. Is there a better way? Is there some kind of numpy compatibility mode out there somewhere?

The difference is that both bounds[:,1] and bounds[:, 1] in Julia always produces a 1d array, which is treated as a “column” vector, whereas in numpy the former produces a 2d “row vectormatrix”.

It sounds like you just want

pop = bounds[:, 1] + m

Also,

bounds = [ [-5.0 5.0]
           [-5.0 5.0]
         ]

is not what you want, I think. You declare a 2d array in Julia with

bounds = [ -5.0 5.0
           -5.0 5.0  ]

whereas you have an array of arrays. (Unlike in Python, which doesn’t have 2d arrays natively so numpy’s 2d arrays act like a “list of lists”.)

PS. Quote code blocks using ```, not using '''

5 Likes