Numpy to Julia conversion questions

General Usage
1

As an exercise I’m converting this Python Differential Evolution tutorial code to Julia (https://machinelearningmastery.com/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?

2

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