# 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]))'
4.72087   -3.63459
4.31493    2.73267
...
-0.169094  -2.07064
-4.84554    4.71972
julia> typeof(pop)
``````

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