Using a Framework to Compare Algorithm Performance

Preamble

In [1]:
# used to create block diagrams
%reload_ext xdiag_magic
%xdiag_output_format svg
    
import numpy as np                   # for multi-dimensional containers
import pandas as pd                  # for DataFrames
import plotly.graph_objects as go    # for data visualisation
import plotly.io as pio              # to set shahin plot layout
import platypus as plat              # multi-objective optimisation framework

pio.templates['shahin'] = pio.to_templated(go.Figure().update_layout(legend=dict(orientation="h",y=1.1, x=.5, xanchor='center'),margin=dict(t=0,r=0,b=40,l=40))).layout.template
pio.templates.default = 'shahin'

Introduction

When preparing to implement multi-objective optimisation experiments, it's often more convenient to use a ready-made framework/library instead of programming everything from scratch. There are many libraries and frameworks that have been implemented in many different programming languages. With our focus on multi-objective optimisation, our choice is an easy one. We will choose Platypus which has a focus on multi-objective problems and optimisation.

Platypus is a framework for evolutionary computing in Python with a focus on multiobjective evolutionary algorithms (MOEAs). It differs from existing optimization libraries, including PyGMO, Inspyred, DEAP, and Scipy, by providing optimization algorithms and analysis tools for multiobjective optimization.

In this section, we will use the Platypus framework to compare the performance of the Non-dominated Sorting Genetic Algorithm II (NSGA-II)1 and the Pareto Archived Evolution Strategy (PAES)2. To do this, we will use them to generate solutions to four problems in the DTLZ test suite3.

Because both of these algorithms are stochastic, meaning that they will produce different results every time they are executed, we will select a sufficient sample size of 30 per algorithm per test problem. We will also use the default configurations for all the test problems and algorithms employed in this comparison. We will use the Hypervolume Indicator (introduced in earlier sections) as our performance metric.

Executing an Experiment and Generating Results

In this section, we will using the Platypus implementation of NSGA-II and PAES to generate solutions for the DTLZ1, DTLZ2, DTLZ3, and DTLZ4 test problems.

First, we will create a list named problems where each element is a DTLZ test problem that we want to use.

In [2]:
problems = [plat.DTLZ1, plat.DTLZ2, plat.DTLZ3, plat.DTLZ4]

Similarly, we will create a list named algorithms where each element is an algorithm that we want to compare.

In [3]:
algorithms = [plat.NSGAII, plat.PAES]

Now we can execute an experiment, specifying the number of function evaluations, $nfe=10,000$, and the number of executions per problem, $seeds=30$. This may take some time to complete depending on your processor speed and the number of function evaluations.

Warning

Running the code below will take a long time to complete even if you have good hardware. To put things into perspective, you will be executing an optimisation process 30 times, per 4 test problems, per 2 algorithms. That's 240 executions of 10,000 function evaluations, totalling in at 2,400,000 function evaluations. You may prefer to change the value of seeds below to something smaller like 5 for now.

We are also using the ProcessPoolEvaluator in Platypus to speed things up.

In [4]:
with plat.ProcessPoolEvaluator(10) as evaluator:
    results = plat.experiment(algorithms, problems, seeds=30, nfe=10000, evaluator=evaluator)

Once the above execution has completed, we can initialize an instance of the hypervolume indicator provided by Platypus.

In [15]:
hyp = plat.Hypervolume(minimum=[0, 0, 0], maximum=[1, 1, 1])

Now we can use the calculate function provided by Platypus to calculate all our hypervolume indicator measurements for the results from our above experiment.

In [16]:
hyp_result = plat.calculate(results, hyp)

Finally, we can display these resultsu using the display function provided by Platypus.

In [17]:
plat.display(hyp_result, ndigits=3)
NSGAII
    DTLZ1
        Hypervolume : [0.802, 0.209, 0.748, 0.794, 0.13, 0.177, 0.334, 0.837, 0.449, 0.002, 0.01, 0.0, 0.039, 0.063, 0.18, 0.44, 0.794, 0.0, 0.602, 0.259, 0.348, 0.0, 0.613, 0.0, 0.835, 0.14, 0.27, 0.442, 0.057, 0.0]
    DTLZ2
        Hypervolume : [0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.208, 0.209, 0.209, 0.209, 0.209, 0.208, 0.209, 0.209, 0.208]
    DTLZ3
        Hypervolume : [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    DTLZ4
        Hypervolume : [0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.0, 0.209, 0.209, 0.209, 0.0, 0.209, 0.209, 0.209, 0.0, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.209, 0.0, 0.209]
PAES
    DTLZ1
        Hypervolume : [0.811, 0.806, 0.718, 0.773, 0.746, 0.619, 0.676, 0.843, 0.851, 0.814, 0.843, 0.723, 0.76, 0.0, 0.662, 0.721, 0.704, 0.691, 0.847, 0.733, 0.577, 0.599, 0.846, 0.619, 0.833, 0.838, 0.786, 0.709, 0.706, 0.832]
    DTLZ2
        Hypervolume : [0.188, 0.175, 0.159, 0.199, 0.195, 0.189, 0.189, 0.197, 0.19, 0.186, 0.182, 0.198, 0.194, 0.197, 0.194, 0.198, 0.182, 0.184, 0.184, 0.191, 0.196, 0.192, 0.19, 0.192, 0.198, 0.189, 0.181, 0.198, 0.179, 0.171]
    DTLZ3
        Hypervolume : [0.0, 0.0, 0.018, 0.0, 0.0, 0.0, 0.0, 0.061, 0.0, 0.0, 0.0, 0.0, 0.01, 0.0, 0.036, 0.0, 0.079, 0.003, 0.0, 0.0, 0.055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.037]
    DTLZ4
        Hypervolume : [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.196, 0.0, 0.0, 0.0, 0.055, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

Statistical Comparison of the Hypervolume Results

Now that we have a data structure that has been populated with results from each execution of the algorithms, we can do a quick statistical comparison to give us some indication as to which algorithmn (NSGA-II or PAES) performs better on each problem.

We can see in the output of display above that the data structure is organised as follows:

  • Algorithm name (e.g. NSGAII)
    • Problem name (e.g. DTLZ1)
      • Performance metric (e.g. Hypervolume)
        • The score for each run (e.g. 30 individual scores).

As a quick test, let's try and get the hypervolume indicator score for the first execution of NSGA-II on DTLZ1.

In [18]:
hyp_result['NSGAII']['DTLZ1']['Hypervolume'][0]
Out[18]:
0.8023119078249699

To further demonstrate how this works, let's also get the hypervolume indicator score for the sixth execution of NSGA-II on DTLZ1.

In [19]:
hyp_result['NSGAII']['DTLZ1']['Hypervolume'][5]
Out[19]:
0.17719629788983282

Finally, let's get the hypervolume indicator scores for all executions of NSGA-II on DTLZ1.

In [20]:
hyp_result['NSGAII']['DTLZ1']['Hypervolume']
Out[20]:
[0.8023119078249699,
 0.20895064088416979,
 0.7482486285062558,
 0.7939235818872432,
 0.13009540700503586,
 0.17719629788983282,
 0.3341790359353015,
 0.8374402507405787,
 0.4490952632404714,
 0.0019587776569011616,
 0.009744344770431082,
 0.0,
 0.038642954483156106,
 0.06344866450427324,
 0.18014094660025184,
 0.4395868015525359,
 0.7941322398615742,
 0.0,
 0.6023653045244652,
 0.2587345548431213,
 0.3484347435364232,
 0.0,
 0.6131583263721389,
 0.0,
 0.8353137888564315,
 0.1400975085584478,
 0.26975429683836716,
 0.4419755380562997,
 0.05717678680365499,
 0.0]

Perfect. Now we can use numpy to calculate the mean hypervolume indicator value for all of our executions of NSGA-II on DTLZ1.

In [21]:
np.mean(hyp_result['NSGAII']['DTLZ1']['Hypervolume'])
Out[21]:
0.31920355305774445

Let's do the same for PAES.

In [22]:
np.mean(hyp_result['PAES']['DTLZ1']['Hypervolume'])
Out[22]:
0.722838349335488

We can see that the mean hypervolume indicator value for PAES on DTLZ1 is higher than that of NSGA-II on DTLZ1. A higher hypervolume indicator value indicates better performance, so we can tentatively say that PAES outperforms NSGA-II on our configuration of DTLZ1 according to the hypervolume indicator. Of course, we haven't determined if this result is statistically significant.

Let's create a DataFrame where each column refers to the mean hypervolume indicator values for the test problems DTLZ1, DTLZ2, DTLZ3, and DTLZ4, and each row represent the performance of an algorithm (in this case, PAES and NSGA-II).

In [23]:
df_hyp_results = pd.DataFrame(index = hyp_result.keys())

for key_algorithm, algorithm in hyp_result.items():
    for key_problem, problem in algorithm.items():
        for hypervolumes in problem['Hypervolume']:
            df_hyp_results.loc[key_algorithm,key_problem] = np.mean(hypervolumes)
            
df_hyp_results
Out[23]:
DTLZ1 DTLZ2 DTLZ3 DTLZ4
NSGAII 0.000000 0.208399 0.00000 0.208582
PAES 0.831512 0.170900 0.03721 0.000000

Now we have an overview of how our selected algorithms performed on the selected test problems according to the hypervolume indicator. Personally, I find it easier to compare algorithm performance when each column represents a different algorithm rather than a problem.

In [24]:
df_hyp_results.transpose()
Out[24]:
NSGAII PAES
DTLZ1 0.000000 0.831512
DTLZ2 0.208399 0.170900
DTLZ3 0.000000 0.037210
DTLZ4 0.208582 0.000000

Without consideration for statistical significance, which algorithm performs best on each test problem?

Conclusion

In this section we have demonstrated how we can compare two popular multi-objective evolutionary algorithms on a selection of four test problems using the hypervolume indicator to measure their performance. In this case, we simply compared the mean of a sample of executions per problem per algorithm without consideration for statistical significance, however, this it is important to take this into account to ensure that any differences haven't occurred by chance.

Exercise

Create your own experiment, but this time include different algorithms and problems and determine which algorithm performs the best on each problem.


  1. Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. A. M. T. (2002). A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE transactions on evolutionary computation, 6(2), 182-197. 

  2. Knowles, J., & Corne, D. (1999, July). The pareto archived evolution strategy: A new baseline algorithm for pareto multiobjective optimisation. In Congress on Evolutionary Computation (CEC99) (Vol. 1, pp. 98-105). 

  3. Deb, K., Thiele, L., Laumanns, M., & Zitzler, E. (2002, May). Scalable multi-objective optimization test problems. In Proceedings of the 2002 Congress on Evolutionary Computation. CEC'02 (Cat. No. 02TH8600) (Vol. 1, pp. 825-830). IEEE. 

Using a Framework to Generate Results

Preamble

In [1]:
# used to create block diagrams
%reload_ext xdiag_magic
%xdiag_output_format svg
    
import numpy as np                   # for multi-dimensional containers
import pandas as pd                  # for DataFrames
import plotly.graph_objects as go    # for data visualisation
import plotly.io as pio              # to set shahin plot layout
import platypus as plat              # multi-objective optimisation framework

pio.templates['shahin'] = pio.to_templated(go.Figure().update_layout(legend=dict(orientation="h",y=1.1, x=.5, xanchor='center'),margin=dict(t=0,r=0,b=40,l=40))).layout.template
pio.templates.default = 'shahin'

Introduction

When preparing to implement multi-objective optimisation experiments, it's often more convenient to use a ready-made framework/library instead of programming everything from scratch. Many libraries and frameworks have been implemented in many different programming languages. With our focus on multi-objective optimisation, our choice is an easy one. We will choose Platypus which has a focus on multi-objective problems and optimisation.

Platypus is a framework for evolutionary computing in Python with a focus on multiobjective evolutionary algorithms (MOEAs). It differs from existing optimization libraries, including PyGMO, Inspyred, DEAP, and Scipy, by providing optimization algorithms and analysis tools for multiobjective optimization.

In this section, we will use the Platypus framework to optimise solutions to an test function using the Non-dominated Sorting Genetic Algorithm II (NSGA-II)1, a fast and elitist multi-objective genetic algorithm.

Approximating DTLZ2 Solutions using NSGA-II

In this section, we will using the Platypus implementation of NSGA-II to generate solutions for the DTLZ2 test problem2. The DTLZ2 test problem comes from the DTLZ test suite, which consists of box-constrained continuous and scalable multi-objective problems.

Using Platypus, we will instantiate a DTLZ2 test problem object using its constructor with no parameters, this means our DTLZ2 problem will using the default configurations.

In [2]:
problem = plat.DTLZ2()

Next, we want to create an instance of the NSGA-II optimiser, and we'll pass in our problem object as a parameter for its configurations.

In [3]:
algorithm = plat.NSGAII(problem)

Now we can execute the optimisation process. Let's give the optimiser a budget of 10,000 function evaluations as a termination criterion. This may take some time to complete depending on your processor speed and the number of function evaluations.

In [4]:
algorithm.run(10000)

Finally, we can display the results. In this case, we will be printing the objective values for each solution in the final population of the above execution.

In [5]:
for solution in algorithm.result:
    print(solution.objectives)
[1.0030951795341247, 5.903936310321151e-06]
[2.5934802815571933e-07, 1.0006403308986553]
[0.9792060155755691, 0.20497867423640143]
[0.8587886056380178, 0.5133791902361882]
[0.9717227988149256, 0.23829158525297897]
[0.051072900089373104, 0.99894640280302]
[0.12431790052429573, 0.9930427096025656]
[0.08032105919012233, 0.9974658987117496]
[0.28482907214559194, 0.9588849501419295]
[0.9976463918569725, 0.07657291729104365]
[0.37242012568865257, 0.9293617945640071]
[0.9844008593905644, 0.1775795721482782]
[0.9948418896682405, 0.10523628728375913]
[0.30867289793095154, 0.9515746224870116]
[0.6932392178420209, 0.7223033599935098]
[0.9172026928698491, 0.40019049047411226]
[0.9319132650514887, 0.3647235656545425]
[0.8733168191155475, 0.4884072898908119]
[0.9270459532575037, 0.3803926907237519]
[0.3561772683505977, 0.9348420209447224]
[0.5738814811444926, 0.8199344160989686]
[0.4817185053464776, 0.8780431260280246]
[0.1015337630553928, 0.9950918657624807]
[0.8994944708848679, 0.43757915511888695]
[0.410532080386212, 0.9121224161633411]
[0.8332512300123446, 0.5543351420456132]
[0.93983914858739, 0.34259015569394285]
[0.6733469498489381, 0.739999541145456]
[0.6199946011074543, 0.797012459902574]
[0.46433413947300767, 0.8876859920003506]
[0.9508490191572935, 0.31223725430232113]
[0.9872121485490931, 0.1603539093761116]
[0.8425062845863036, 0.53923038628592]
[0.5938240313367141, 0.8048857562401942]
[0.7021603149056569, 0.7178645594075399]
[0.027385609798355002, 0.9999250723438612]
[0.49584426225763983, 0.870272072410772]
[0.15028597157189252, 0.9889820125483831]
[0.990655840054075, 0.13929543412778958]
[0.3389923042542389, 0.9452433727595364]
[1.000128138034012, 0.021291039769240103]
[0.5324957443903208, 0.8472304106988626]
[0.8931824963745468, 0.4520543478221232]
[0.9571251936820613, 0.2912840321949625]
[0.8655486443299237, 0.5013900284094054]
[0.9993451639545191, 0.056098645817105545]
[0.5187540816531349, 0.8554554144812964]
[0.2062940443994896, 0.9801965492912985]
[0.9932921437977046, 0.12451816611150038]
[0.551669929402978, 0.8344449521807499]
[0.5086549818901388, 0.8622236606935018]
[0.9646687340279186, 0.26744387225973926]
[0.8834909492703003, 0.4690551674438777]
[0.17816297905714362, 0.984474776093553]
[0.9682825140953177, 0.2526742892109481]
[0.7269805055137187, 0.6873888694642003]
[0.43170143426663726, 0.9057554637618657]
[0.7525299719551056, 0.6589728085310002]
[0.9999571742683206, 0.03957153636673865]
[0.7152175092536589, 0.7004185660991041]
[0.8064630407109789, 0.5920667055359434]
[0.39391108272493, 0.9197043766949209]
[0.7999027837468777, 0.6010554147687575]
[0.7796760739485277, 0.6267042198308959]
[0.9106622777774218, 0.42120895418319276]
[0.7417202606373894, 0.6716167945690665]
[0.8145085470986322, 0.5809058102253095]
[0.7700482859219662, 0.6392441068031863]
[0.2281748888692325, 0.9741399780726558]
[0.9602141302228093, 0.2816229964682574]
[0.8254947622684579, 0.5716030505342371]
[0.5445222235432747, 0.8407582760681087]
[0.24833295326598553, 0.9688662152627027]
[0.7371649392970931, 0.6772689523072278]
[0.26719727820351075, 0.9641025961709829]
[0.644962173109967, 0.7646863917078949]
[0.16155465051863602, 0.9876239102454765]
[0.008546188630680702, 1.0002334463830467]
[0.2104248125964668, 0.9778834333435223]
[0.6691117548441301, 0.744674651953048]
[0.9107992768150537, 0.4149380869292202]
[0.8862464638848153, 0.4643450636538065]
[0.45109118717734475, 0.892913130993857]
[0.7583512783776537, 0.6526106608932809]
[0.6352227222947703, 0.7742703529899837]
[0.23761154772962634, 0.9727693743175875]
[0.5826879942516339, 0.8163002078601734]
[0.9462522772480144, 0.32579445375801847]
[0.8419582647266676, 0.5407927345929345]
[1.0009542178225288, 5.91323534837545e-06]
[0.39639110717449894, 0.918638015272697]
[0.43769304741927517, 0.9033022593261911]
[0.4455930080304075, 0.8955379997841134]
[0.26123977151089983, 0.965712363634355]
[0.6309897379482855, 0.7763307824072477]
[0.944671546145947, 0.33066994390004645]
[0.6590880040993868, 0.7535741318802496]
[0.32622335743809966, 0.9487338831352291]
[0.588589885991485, 0.8097323554356354]
[0.6516972023490352, 0.7605102585693663]

Visualising Solutions in Objective Space

In the last section, we concluded by printing the objective values for every solution. This information will be easier to digest using a plot, so let's quickly put the data into a Pandas DataFrame and then use Plotly to create a scatterplot. Let's start by moving our Platypus data structure to a DataFrame.

In [6]:
objective_values = np.empty((0, 2))

for solution in algorithm.result:
    y = solution.objectives
    objective_values = np.vstack([objective_values, y])
    
# convert to DataFrame
objective_values = pd.DataFrame(objective_values, columns=['f1','f2'])

All done, we can have a peek at our DataFrame to make sure we've not made any obvious mistakes.

In [7]:
objective_values
Out[7]:
f1 f2
0 1.003095e+00 0.000006
1 2.593480e-07 1.000640
2 9.792060e-01 0.204979
3 8.587886e-01 0.513379
4 9.717228e-01 0.238292
... ... ...
95 9.446715e-01 0.330670
96 6.590880e-01 0.753574
97 3.262234e-01 0.948734
98 5.885899e-01 0.809732
99 6.516972e-01 0.760510

100 rows × 2 columns

Looking good. All that we have left is the scatter plot.

In [9]:
fig = go.Figure()
fig.add_scatter(x=objective_values.f1, y=objective_values.f2, mode='markers')
fig.show()

Great! If you search for the true Pareto-optimal front for DTLZ2, you can see that our approximation is looking as expected.

Conclusion

In this section we have demonstrated how we can use a popular multi-objective optimisation algorithm, NSGA-II, to approximate multiple trade-off solutions to the DTLZ2 test problem. We did this using the Platypus framework, and also moved the results into a DataFrame in order to easily create a scatterplot.

Exercise

Repeat the same experiment in this section, but this time select a different algorithm that has been implemented in Platypus.


  1. Deb, K., Pratap, A., Agarwal, S., & Meyarivan, T. A. M. T. (2002). A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE transactions on evolutionary computation, 6(2), 182-197. 

  2. Deb, K., Thiele, L., Laumanns, M., & Zitzler, E. (2002, May). Scalable multi-objective optimization test problems. In Proceedings of the 2002 Congress on Evolutionary Computation. CEC'02 (Cat. No. 02TH8600) (Vol. 1, pp. 825-830). IEEE.. 

Docker Toolbox, Windows 7, and Shared Volumes

Introduction

Docker is a container platform with so much to offer, and I've already published a quick guide on how to set it up in an earlier article. However, if you don't meet the system requirements for the newer versions then you will have to rely on Docker Toolbox, which is described as a legacy desktop solution:

Docker Toolbox is for older Mac and Windows systems that do not meet the requirements of Docker for Mac and Docker for Windows. We recommend updating to the newer applications, if possible.

Docker Toolbox

To operate, Docker Toolbox installs a VirtualBox virtual machine named default, running the boot2docker Linux distribution. Because of this additional layer of virtualisation, additional steps are needed when sharing volumes between the host operating system and a Docker container. For example, a computer lab that I use for teaching runs Windows 7 with restricted student access. I have convinced our IT team to install Docker Toolbox, but my students still need some way of sharing a folder on their student drive (herein referred to as the H:\ drive) with the Docker container.

Preparing VirtualBox to Share Folders

To achieve this, the folder (e.g. H:\work) needs to be shared with the VirtualBox virtual machine named default, before it can be mounted as a volume within a Docker container. There may be other approaches, but here are the steps for our particular set-up:

  1. Open "Docker Quickstart Terminal".
  2. Once Docker is running, type docker-machine stop default.
  3. Open the Command Line from the start menu (search for cmd.exe).
  4. Navigate to the VirtualBox folder using cd C:/Program Files/Oracle/VirtualBox
  5. Type VBoxManage.exe sharedfolder add default --name "h/work" --hostpath "\\?\h:\work" --automount
  6. Type VBoxManage.exe setextradata default VBoxInternal2/SharedFoldersEnableSymlinksCreate/v-root 1
  7. Type VBoxManage.exe setextradata default VBoxInternal2/SharedFoldersEnableSymlinksCreate/h/work 1

The above steps will prepare VirtualBox for sharing H:\work. Now we need to mount the directory inside our docker-machine:

  1. Open "Docker Quickstart Terminal".
  2. Once Docker is running, type docker-machine ssh default
  3. Type sudo mkdir --parents /h/work
  4. Type sudo mount -t vboxsf h/work /h/work/
  5. Type exit

You can now create Docker container instances and share the H:\work folder with them as a mounted volume. Let's demonstrate this with jupyter/scipy-notebook

Sharing Folders with a Docker Container

To create a Docker container from the jupyter/scipy-notebook image, type the following command and wait for it to complete execution: docker run --name="scipy" --user root -v /h/work:/home/jovyan -d -e GRANT_SUDO=yes -p 8888:8888 jupyter/scipy-notebook start-notebook.sh --NotebookApp.token=''

This may take some time, as it will need to download and extract the image. Once it's finished, you should be able to access the Jupyter notebook using 127.0.0.1:8888. I hope this helps you get up and running with Docker Toolbox and shared folders. Of course, the process is typically easier when using the non-legacy Docker solutions.

Comments from the Old Blog

Docker Comments

Quick Guide to Docker

Introduction

Docker is a fantastic platform which has changed my workflow for the better. Before Docker, I was using VMWare images and linked-clones to create isolated environments for my experiments. Sharing these was not easy - the file sizes were too large and not everyone has a VMWare license. Performance wasn't great either, and VMWare updates would sometimes cause problems with the image. Before VMWare, I was using MATLAB and directly sharing the scripts for my experiments. This was an unpleasant experience, as each of the many versions of MATLAB (e.g. R2011a, R2011b, R2012a, R2012b, etc.) would break some functionality from one of the toolboxes. Again, not everyone could afford a MATLAB license, and I only had mine as part of my institution.

With Docker, I've moved to using Python's scientific stack and libraries such as TensorFlow and Keras. It took a while to get used to living without the MATLAB IDE, but there was a similar IDE called Spyder, which is a Scientific Python Development EnviRonment. However, I've gotten over the MATLAB IDE since discovering Jupyter. With all these technologies, I can share an environment identical to the one I executed my experiments in, and allow a colleague or reviewer to do the same. My students can also get up and running with the environments needed by my taught units without having to configure everything themselves. The only difficulty I've encountered so far is when someone has a Windows installation which doesn't support Docker for Windows, as Docker Toolbox isn't always a smooth experience.

I've created this video as a brief introduction to Docker, complete with instructions on how to install Docker on macOS, Windows 10, and Ubuntu 16.04.

The aim is to get someone started and tinkering with Docker quickly.

Common Problems

Kernel appears to have died

If your kernel dies when training a network using TensorFlow, then you might need to increase the amount of RAM allocated to Docker.

Kernel Restarting

You can do this in the Docker Preferences:

Increase Memory

Minimum requirements not met

I have Windows 10 Pro/Enterprise/Education, but the installer says I don't meet the minimum requirements

You need to be on Windows 10 build 10586 or newer. You can check your Windows build number by clicking the start-bar and winver, and pressing the enter key:

Windows Build

Contributing Hypervolume Indicator

Preamble

In [1]:
import numpy as np                   # for multi-dimensional containers
import pandas as pd                  # for DataFrames
import plotly.graph_objects as go    # for data visualisation
import plotly.io as pio              # to set shahin plot layout
import platypus as plat              # multi-objective optimisation framework
import pygmo as pg                   # multi-objective optimisation framework
import plotly.express as px

pio.templates['shahin'] = pio.to_templated(go.Figure().update_layout(legend=dict(orientation="h",y=1.1, x=.5, xanchor='center'),margin=dict(t=0,r=0,b=40,l=40))).layout.template
pio.templates.default = 'shahin'

Introduction

In this section, we're going to take a look at how to calculate the contributing hypervolume indicator values for a population of solutions.

The Contributing Hypervolume (CHV) indicator is a population sorting mechanism based on an adaptation of the hypervolume indicator. The hypervolume indicator works by calculating the size of the objective space that has been dominated by an entire approximation set with respect to a specified reference point, whereas the CHV indicator assigns each solution in an approximation set with the size of the space that has been dominated by each solution exclusively. With this information, the population can be sorted by the most dominant and diverse solutions. This has been illustrated below in two-dimensional space with a population of three solutions.

Contributing Hypervolume Indicator

Calculating the exact CHV indicator is attractively simple. The method begins by first calculating the hypervolume indicator quality of a population $\mathrm{X}$, and then for each solution in the population, the solution is removed and the hypervolume indicator quality is again calculated for the new population. The new hypervolume indicator value is then subtracted from the hypervolume indicator value of the whole population, which results in the CHV indicator value of the solution which was removed. It is then possible to calculate the CHV indicator values of all the solutions in the population, order them by descending value so that they are ordered by the greatest explicit hypervolume indicator contribution, and select the better solutions to form the next parent population. This approach has been listed in Algorithm below.

$$ \begin{aligned} &\textbf{CHVIndicator}(f^{ref}, X)\\ &\;\;\;\;X_{HV} \leftarrow HV(f^{ref},X)\\ &\;\;\;\;for\;\;{n = 1 : \lambda}\;\;do\\ &\;\;\;\;X_t \leftarrow X \backslash {X_n}\\ &HV_n \leftarrow HV(f^{ref},X_t)\\ &CHV_n \leftarrow X_{HV} - HV_n\\ &\textbf{return}\;\;CHV \end{aligned} $$

Although many optimisation algorithms use the CHV as a sorting criterion for selection, its calculation becomes computationally infeasible as the number of objectives being considered increase. Monte Carlo approximations have been used to speed up the calculation of the CHV in which through empirical experiments has shown that the method does not impair the quality of the approximation set. However, the speed increase provided by the Monte Carlo approximation method still results in an infeasibility of the CHV indicator on problems consisting of five objectives or more.

This particular measure of diversity preservation can also be used to reduce the size of a final approximation set produced by an optimiser, to a size that will not overwhelm and confuse a decision-maker.

CHV Dependence

It is important to note the dependence of the CHV indicator on the population on which it is being calculated. To demonstrate this a population initially consisting of three solutions in two-objective space has been created, these solutions are $A:(1,10)$, $B:(5,3)$, $C:(6,2)$. The figure below illustrates the explicit HV which is contributed by each solution in the population, and it can be observed that in this population solution B contributes the most to the overall HV indicator value. This can be seen clearly by the number of squares shaded by the corresponding colour, and in this case, solution B covers seven squares.

CHV with three solutions

In order to properly demonstrate this dependence, a solution has been inserted into the population used in the previous example. The solutions for the following example are now $A:(1,10)$, $B:(5,3)$, $C:(6,2)$, $D:(4,4)$. The figure below illustrates the explicit HV which is contributed by each solution in this updated population, and it can be observed that in this population solution B now contributes the least to the overall HV indicator value. This can again be seen clearly by the number of squares shaded by the corresponding colour, and in this case, solution B now covers a single square. This is because the new solution, D, offers dominance over new areas of the objective space as well as objective space which is mutually dominated by solution B. As such, the objective which is mutually dominated by both of these solutions is no longer considered an explicit contribution by any one of these solutions.

CHV with four solutions

Calculating the Contributing Hypervolume Indicator with Platypus

Let's define some necessary variables before invoking the Platypus implementation of the hypervolume indicator algorithm to calculate the CHV of each solution. We will use the ZDT1 test function with the number of design variables $\mathrm{D}=30$ throughout this example, and with the population size $\mathrm{N}=100$.

In [2]:
problem = plat.ZDT1()
D = 30
N = 100

With these variables defined, we will now move onto generating our initial population. We will be using Platypus Solution objects for this, which we will initialise with random problem variables, evaluate, and then append to a list named solutions.

In [3]:
solutions = []

for i in range(N):
    solution = plat.Solution(problem)
    solution.variables = np.random.rand(D)
    solution.evaluate()
    solutions.append(solution)

Let's print out the variables and objectives for the first solution in this list to see what they look like.

In [4]:
print(f"Design variables:\n {solutions[0].variables}")
print(f"Objective values:\n {solutions[0].objectives}")
Design variables:
 [0.15075328 0.67733039 0.28828798 0.77994058 0.45665193 0.57820741
 0.56201372 0.97177937 0.27651144 0.71249518 0.84085362 0.4039491
 0.64843207 0.62018162 0.55052643 0.15276781 0.95331894 0.12605556
 0.33930358 0.66852502 0.8936662  0.58917854 0.23363381 0.14102791
 0.62399114 0.96264146 0.67394996 0.47540925 0.86808249 0.12328518]
Objective values:
 [0.15075328001905353, 5.072052879064982]

Now that we have a population of solutions stored in the solutions variable, we can prepare an instance of the Platypus.indicators.Hypervolume() object with the desired reference point for the hypervolume indicator calculation. For ZDT1, the reference point is typically $\langle11,11\rangle$.

In [5]:
hyp = plat.indicators.Hypervolume(minimum=[0, 0], maximum=[11, 11])

We can now use this hyp object to calculate the hypervolume indicator for any population.

Note

The Platypus implementation of the hypervolume indicator requires either a minimum and a maximum point, or a reference_set (not the same as the hypervolume reference point). Normally, a hypervolume indicator algorithm would only require a single vector that defines the reference point $f ^{ref}$. In the case of Platypus, $f ^{ref}$ actually corresponds to maximum, but Platypus also forces us to provide a vector for minimum, which we have set to $\langle0, 0\rangle$

Let's calculate the hypervolume indicator value for the population of solutions we created above and named solution. We'll store this in HV.

In [6]:
HV = hyp.calculate(solutions)
print(f"Hypervolume indicator value: {HV}")
Hypervolume indicator value: 0.7598049263738476

We now have this single hypervolume indicator value that we will use to calculate the explicit hypervolume contribution of each solution in the population.

In [7]:
CHV = np.empty((0, 1))

for i in range(N):
    solutions_subset = solutions[:i] + solutions[i+1:]
    CHV = np.vstack([CHV, HV - hyp.calculate(solutions_subset)])
    print(f"CHV of solution {i}:\t{CHV[i][0]}")
CHV of solution 0:	0.0
CHV of solution 1:	0.0
CHV of solution 2:	1.1102230246251565e-16
CHV of solution 3:	0.0002092519790036773
CHV of solution 4:	0.0
CHV of solution 5:	0.0
CHV of solution 6:	0.0
CHV of solution 7:	0.0
CHV of solution 8:	0.0
CHV of solution 9:	0.0
CHV of solution 10:	0.0
CHV of solution 11:	0.0
CHV of solution 12:	0.0
CHV of solution 13:	0.0
CHV of solution 14:	0.0
CHV of solution 15:	0.0
CHV of solution 16:	0.0
CHV of solution 17:	0.0
CHV of solution 18:	5.576591804534736e-05
CHV of solution 19:	7.796603141430047e-06
CHV of solution 20:	0.0
CHV of solution 21:	-1.1102230246251565e-16
CHV of solution 22:	0.0
CHV of solution 23:	1.1102230246251565e-16
CHV of solution 24:	0.0
CHV of solution 25:	0.0
CHV of solution 26:	5.1443529762384976e-05
CHV of solution 27:	0.0
CHV of solution 28:	0.0
CHV of solution 29:	-1.1102230246251565e-16
CHV of solution 30:	0.0
CHV of solution 31:	0.0
CHV of solution 32:	0.0
CHV of solution 33:	0.00020671258311477647
CHV of solution 34:	0.0
CHV of solution 35:	0.0
CHV of solution 36:	1.1102230246251565e-16
CHV of solution 37:	0.0
CHV of solution 38:	0.0012109326868405823
CHV of solution 39:	0.0
CHV of solution 40:	0.0
CHV of solution 41:	0.0
CHV of solution 42:	0.0
CHV of solution 43:	0.0
CHV of solution 44:	0.0
CHV of solution 45:	0.0
CHV of solution 46:	0.0
CHV of solution 47:	0.0
CHV of solution 48:	-1.1102230246251565e-16
CHV of solution 49:	1.1102230246251565e-16
CHV of solution 50:	0.0
CHV of solution 51:	0.0
CHV of solution 52:	0.0
CHV of solution 53:	0.00016943666158597548
CHV of solution 54:	-1.1102230246251565e-16
CHV of solution 55:	0.0
CHV of solution 56:	1.1102230246251565e-16
CHV of solution 57:	0.0
CHV of solution 58:	-1.1102230246251565e-16
CHV of solution 59:	1.1102230246251565e-16
CHV of solution 60:	0.0
CHV of solution 61:	0.0
CHV of solution 62:	0.0
CHV of solution 63:	0.0
CHV of solution 64:	0.0
CHV of solution 65:	0.0
CHV of solution 66:	0.0
CHV of solution 67:	0.0
CHV of solution 68:	1.1102230246251565e-16
CHV of solution 69:	0.0
CHV of solution 70:	0.0
CHV of solution 71:	0.0
CHV of solution 72:	1.1102230246251565e-16
CHV of solution 73:	0.0
CHV of solution 74:	1.1102230246251565e-16
CHV of solution 75:	0.0
CHV of solution 76:	0.002428824169973476
CHV of solution 77:	0.0
CHV of solution 78:	0.0
CHV of solution 79:	0.0
CHV of solution 80:	-1.1102230246251565e-16
CHV of solution 81:	0.0
CHV of solution 82:	0.0
CHV of solution 83:	0.0
CHV of solution 84:	-1.1102230246251565e-16
CHV of solution 85:	1.1102230246251565e-16
CHV of solution 86:	-1.1102230246251565e-16
CHV of solution 87:	0.0
CHV of solution 88:	0.0
CHV of solution 89:	1.1102230246251565e-16
CHV of solution 90:	0.0
CHV of solution 91:	0.0
CHV of solution 92:	0.0
CHV of solution 93:	0.0
CHV of solution 94:	0.0
CHV of solution 95:	0.0
CHV of solution 96:	0.0
CHV of solution 97:	0.0
CHV of solution 98:	1.1102230246251565e-16
CHV of solution 99:	0.0

With the above list, we can consider anything with a CHV of $0$ or around the precision of $16$ positions behind the decimal as a solution that does not contribute to the hypervolume quality of the population.

Calculating the Contributing Hypervolume Indicator with PyGMO

We can also use a different framework's implementation of the hypervolume indicator on our population. We should be expecting the same values, but this is a good exercise to learn how to use a different framework, and perhaps to check that they do indeed arrive at the same value.

This time we will use the PyGMO framework. With PyGMO, we are lucky enough to have a built-in function to calculate the CHV of each solution. PyGMO's hypervolume indicator function can work with a few different data-types, including numpy.array(). We have previously moved our Platypus solutions to a pandas.DataFrame (which can easily be output as a numpy.array()). Let's begin by creating a new DataFrame with the columns f1 and f2 which will be used to store our objective values for each solution.

In [8]:
solutions_df = pd.DataFrame(index=range(N),columns=['f1','f2']).astype(float)
solutions_df.head()
Out[8]:
f1 f2
0 NaN NaN
1 NaN NaN
2 NaN NaN
3 NaN NaN
4 NaN NaN

We can see that we've also defined an index range that covers the number of solutions in our population, $\mathrm{N}=100$. This means we have $100$ rows ready, but their values are initialised to NaN (Not A Number), which in this case simply indicates missing data.

Let's use a loop to iterate through our solutions list of $100$ solutions and assign the desired values to the corresponding row in our solutions_df DataFrame

In [9]:
for i in range(N):
    solutions_df.loc[i].f1 = solutions[i].objectives[0]
    solutions_df.loc[i].f2 = solutions[i].objectives[1]
        
solutions_df.head()
Out[9]:
f1 f2
0 0.150753 5.072053
1 0.197567 5.161087
2 0.291905 3.832387
3 0.002193 6.197517
4 0.577217 3.750645

We can see our DataFrame now contains the desired values. We can now easily get this data as a numpy.array() to feed into PyGMO's hypervolume indicator object constructor.

In [10]:
hyp = pg.hypervolume(solutions_df[['f1','f2']].values)

Now we can invoke the exclusive() method on our hypervolume object and pass in the reference point to calculate the CHV values.

In [11]:
for i in range(N):
    print(f"CHV of solution {i}:\t{hyp.exclusive(i, [11, 11]) / np.prod([11, 11])}")
CHV of solution 0:	0.0
CHV of solution 1:	0.0
CHV of solution 2:	0.0
CHV of solution 3:	0.00020925197900382504
CHV of solution 4:	0.0
CHV of solution 5:	0.0
CHV of solution 6:	0.0
CHV of solution 7:	0.0
CHV of solution 8:	0.0
CHV of solution 9:	0.0
CHV of solution 10:	0.0
CHV of solution 11:	0.0
CHV of solution 12:	0.0
CHV of solution 13:	0.0
CHV of solution 14:	0.0
CHV of solution 15:	0.0
CHV of solution 16:	0.0
CHV of solution 17:	0.0
CHV of solution 18:	5.576591804538957e-05
CHV of solution 19:	7.796603141359396e-06
CHV of solution 20:	0.0
CHV of solution 21:	0.0
CHV of solution 22:	0.0
CHV of solution 23:	0.0
CHV of solution 24:	0.0
CHV of solution 25:	0.0
CHV of solution 26:	5.1443529762516185e-05
CHV of solution 27:	0.0
CHV of solution 28:	0.0
CHV of solution 29:	0.0
CHV of solution 30:	0.0
CHV of solution 31:	0.0
CHV of solution 32:	0.0
CHV of solution 33:	0.0002067125831148416
CHV of solution 34:	0.0
CHV of solution 35:	0.0
CHV of solution 36:	0.0
CHV of solution 37:	0.0
CHV of solution 38:	0.0012109326868405346
CHV of solution 39:	0.0
CHV of solution 40:	0.0
CHV of solution 41:	0.0
CHV of solution 42:	0.0
CHV of solution 43:	0.0
CHV of solution 44:	0.0
CHV of solution 45:	0.0
CHV of solution 46:	0.0
CHV of solution 47:	0.0
CHV of solution 48:	0.0
CHV of solution 49:	0.0
CHV of solution 50:	0.0
CHV of solution 51:	0.0
CHV of solution 52:	0.0
CHV of solution 53:	0.00016943666158564976
CHV of solution 54:	0.0
CHV of solution 55:	0.0
CHV of solution 56:	0.0
CHV of solution 57:	0.0
CHV of solution 58:	0.0
CHV of solution 59:	0.0
CHV of solution 60:	0.0
CHV of solution 61:	0.0
CHV of solution 62:	0.0
CHV of solution 63:	0.0
CHV of solution 64:	0.0
CHV of solution 65:	0.0
CHV of solution 66:	0.0
CHV of solution 67:	0.0
CHV of solution 68:	0.0
CHV of solution 69:	0.0
CHV of solution 70:	0.0
CHV of solution 71:	0.0
CHV of solution 72:	0.0
CHV of solution 73:	0.0
CHV of solution 74:	0.0
CHV of solution 75:	0.0
CHV of solution 76:	0.0024288241699734446
CHV of solution 77:	0.0
CHV of solution 78:	0.0
CHV of solution 79:	0.0
CHV of solution 80:	0.0
CHV of solution 81:	0.0
CHV of solution 82:	0.0
CHV of solution 83:	0.0
CHV of solution 84:	0.0
CHV of solution 85:	0.0
CHV of solution 86:	0.0
CHV of solution 87:	0.0
CHV of solution 88:	0.0
CHV of solution 89:	0.0
CHV of solution 90:	0.0
CHV of solution 91:	0.0
CHV of solution 92:	0.0
CHV of solution 93:	0.0
CHV of solution 94:	0.0
CHV of solution 95:	0.0
CHV of solution 96:	0.0
CHV of solution 97:	0.0
CHV of solution 98:	0.0
CHV of solution 99:	0.0

With the knowledge that Platypus normalises the HV values, we have already normalised PyGMO's results by dividing by the product of the reference point. Now we can see that both frameworks have arrived at the same hypervolume indicator value, although PyGMO has handled the precision issue differently.

Conclusion

In this section, we have introduced the contributing hypervolume indicator as a criterion that can be used in the selection of solutions from a population. With CHV values assigned to solutions in a population, we can, for example, sort them by descending order and select the top 10 solutions to get the 10 solutions that contribute most to the population. We also demonstrated the application of two implementations of the contributing hypervolume indicator, one in Platypus, and one in PyGMO.

Exercise

Try updating the code for the Platypus CHV calculation to address the precision issue.

Exercise

Try using the CHV values for each solution to create a visualisation which highlights the top $T$ (e.g. $T=5$) solutions in a non-dominated front.


  1. E. Zitzler, S. K ̈unzli, Indicator-based selection in multiobjective search, in: Parallel Problem Solving from Nature-PPSNVIII, Springer, 2004, pp. 832–842 

Hypervolume Indicator

Preamble

In [1]:
import numpy as np                   # for multi-dimensional containers
import pandas as pd                  # for DataFrames
import plotly.graph_objects as go    # for data visualisation
import plotly.io as pio              # to set shahin plot layout
import platypus as plat              # multi-objective optimisation framework
import pygmo as pg                   # multi-objective optimisation framework
import plotly.express as px

pio.templates['shahin'] = pio.to_templated(go.Figure().update_layout(legend=dict(orientation="h",y=1.1, x=.5, xanchor='center'),margin=dict(t=0,r=0,b=40,l=40))).layout.template
pio.templates.default = 'shahin'

Introduction

In this section we're going to take a look at how to calculate the hypervolume indicator value for a population of solutions. The hypervolume indicator (or $s$-metric) is a performance metric for indicating the quality of a non-dominated approximation set, introduced by [1] where it is described as the "size of the space covered or size of dominated space". It can be defined as:

$$ HV\big(f ^{ref}, \mathrm{X}\big) = \Lambda \left( \bigcup_{\mathrm{X}_n \in \mathrm{X} } \Big[ f_{1}(\mathrm{X}_n), f_{1}^{ref} \Big] \times \dots \times \Big[ f_{m}(\mathrm{X}_n), f_{m}^{ref} \Big] \right) $$

where $HV\left(f ^{ref}, \mathrm{X}\right)$ resolves the size of the space covered by an approximation set $\mathrm{X}$, $f^{ref} \in \mathbb{R}$ refers to a chosen reference point, and $\Lambda \left(.\right)$ refers to the Lebesgue measure. This can be illustrated in two-dimensional objective space (to allow for easy visualisation) with a population of three solutions.

Hypervolume Indicator

The hypervolume indicator is appealing because it is compatible with any number of problem objectives and requires no prior knowledge of the true Pareto-optimal front, this is important when working with real-world problems which have not yet been solved. The hypervolume indicator is currently used in the field of multi-objective optimisation as both a proximity and diversity performance metric, as well as in the decision-making process.

Unlike dominance-based criteria which require only two solutions for performing a comparison (which can be used on an ordinal scale), a reference vector is required to calculate the HV indicator value (i.e. it requires the objective to be measured on an interval scale). When used for pairwise or multiple comparison of optimisation algorithms, this reference vector must be the same, otherwise, the resulting HV indicator values are not comparable. This reference vector can be approximated as large values for each problem objective in order for all objective values in any approximation set to be within the reference vector.

Calculating the Hypervolume Indicator with Platypus

Let's define some necessary variables before invoking the Platypus implementation of the hypervolume indicator algorithm. We will use the ZDT1 test function with the number of design variables $\mathrm{D}=30$ throughout this example, and with the population size $\mathrm{N}=100$.

In [2]:
problem = plat.ZDT1()
D = 30
N = 100

With these variables defined, we will now move onto generating our initial population. We will be using Platypus Solution objects for this, which we will initialise with random problem variables, evaluate, and then append to a list named solutions.

In [3]:
solutions = []

for i in range(N):
    solution = plat.Solution(problem)
    solution.variables = np.random.rand(D)
    solution.evaluate()
    solutions.append(solution)

Let's print out the variables and objectives for the first solution in this list to see what they look like.

In [4]:
print(f"Design variables:\n {solutions[0].variables}")
print(f"Objective values:\n {solutions[0].objectives}")
Design variables:
 [0.60334361 0.87787542 0.03142044 0.14487709 0.92481308 0.95661311
 0.33349403 0.37519477 0.03984425 0.01261005 0.67984533 0.2307125
 0.96553172 0.15568727 0.50446299 0.51060793 0.56343829 0.59170951
 0.32728624 0.02539182 0.92936501 0.32965463 0.95077326 0.80808745
 0.10041791 0.9403275  0.61565891 0.39527713 0.56914812 0.59565309]
Objective values:
 [0.6033436096584304, 3.674672751755749]

Now that we have a population of solutions stored in the solutions variable, we can prepare an instance of the Platypus.indicators.Hypervolume() object with the desired reference point for the hypervolume indicator calculation. For ZDT1, the reference point is typically $\langle11,11\rangle$.

In [5]:
hyp = plat.indicators.Hypervolume(minimum=[0, 0], maximum=[11, 11])

We can now use this hyp object to calculate the hypervolume indicator for any population.

Note

The Platypus implementation of the hypervolume indicator requires either a minimum and a maximum point, or a reference_set (not the same as the hypervolume reference point). Normally, a hypervolume indicator algorithm would only require a single vector that defines the reference point $f ^{ref}$. In the case of Platypus, $f ^{ref}$ actually corresponds to maximum, but Platypus also forces us to provide a vector for minimum, which we have set to $\langle0, 0\rangle$

Let's calculate the hypervolume indicator value for the population of solutions we created above and named solution.

In [6]:
print(f"Hypervolume indicator value: {hyp.calculate(solutions)}")
Hypervolume indicator value: 0.7554161688220257

We now have this single hypervolume indicator value that we can use to gauge and compare the performance of our population. The higher the value is, the "better" the hypervolume quality.

Calculating the Hypervolume Indicator with PyGMO

We can also use a different framework's implementation of the hypervolume indicator on our population. We should be expecting the same value, but this is a good exercise to learn how to use a different framework, and perhaps to check that they do indeed arrive at the same value.

This time we will use the PyGMO framework. PyGMO's hypervolume indicator function can work with a few different data-types, including numpy.array(). We have previously moved our Platypus solutions to a pandas.DataFrame (which can easily be output as a numpy.array()). Let's begin by creating a new DataFrame with the columns f1 and f2 which will be used to store our objective values for each solution.

In [7]:
solutions_df = pd.DataFrame(index=range(N),columns=['f1','f2']).astype(float)
solutions_df.head()
Out[7]:
f1 f2
0 NaN NaN
1 NaN NaN
2 NaN NaN
3 NaN NaN
4 NaN NaN

We can see that we've also defined an index range that covers the number of solutions in our population, $\mathrm{N}=100$. This means we have $100$ rows ready, but their values are initialised to NaN (Not A Number), which in this case simply indicates missing data.

Let's use a loop to iterate through our solutions list of $100$ solutions and assign the desired values to the corresponding row in our solutions_df DataFrame

In [8]:
for i in range(N):
    solutions_df.loc[i].f1 = solutions[i].objectives[0]
    solutions_df.loc[i].f2 = solutions[i].objectives[1]
        
solutions_df.head()
Out[8]:
f1 f2
0 0.603344 3.674673
1 0.806312 2.628754
2 0.059896 5.468507
3 0.322516 4.536322
4 0.928396 2.653431

We can see our DataFrame now contains the desired values. We can now easily get this data as a numpy.array() to feed into PyGMO's hypervolume indicator object constructor.

In [9]:
hyp = pg.hypervolume(solutions_df[['f1','f2']].values)

Now we can invoke the compute() method on our hypervolume object and pass in the reference point to calculate the value.

In [11]:
hyp.compute([11, 11])
Out[11]:
91.40535642746507

Upon first inspection, this value seems to be different to the one calculated by Platypus so we may decide that we've either done something wrong, or there is a mistake in at least one of the implementations. However, in this case, it is simply the case that one of them (Platypus) has normalised the output. We can do the same with the output from PyGMO by dividing the hypervolume indicator value by the product of the reference point.

In [13]:
hyp.compute([11, 11]) / np.prod([11, 11])
Out[13]:
0.7554161688220253

Now we can see that both frameworks have arrived at the same hypervolume indicator.

Conclusion

In this section, we have introduced the hypervolume indicator as a criterion that can be used in the selection of a population. We also demonstrated the application of two implementations of the hypervolume indicator, one in Platypus, and one in PyGMO.

Exercise

Experiment with the hypervolume indicator: try using different reference points, different population sizes, and different problems.


  1. E. Zitzler, S. K ̈unzli, Indicator-based selection in multiobjective search, in: Parallel Problem Solving from Nature-PPSNVIII, Springer, 2004, pp. 832–842 

Non-Dominated Sorting

Preamble

In [1]:
# used to create block diagrams
%reload_ext xdiag_magic
%xdiag_output_format svg
    
import numpy as np                   # for multi-dimensional containers
import pandas as pd                  # for DataFrames
import plotly.graph_objects as go    # for data visualisation
import plotly.io as pio              # to set shahin plot layout
import platypus as plat              # multi-objective optimisation framework
import plotly.express as px

pio.templates['shahin'] = pio.to_templated(go.Figure().update_layout(legend=dict(orientation="h",y=1.1, x=.5, xanchor='center'),margin=dict(t=0,r=0,b=40,l=40))).layout.template
pio.templates.default = 'shahin'

Introduction

In this section, we're going to take a look at how to use the Platypus framework to apply the non-dominated sorting algorithm to a population of randomly generated solutions. Non-dominated sorting is important during the selection stage of an Evolutionary Algorithm because it allows us to prioritise the selection of solutions based on their dominance relations with respect to the rest of the population. For example, we may wish to only select 25 solutions from a population of 50, with the intention to use these 25 solutions in the variation stage that follows.

Non-Dominated Sorting with Platypus

Let's define some necessary variables before invoking the Platypus implementation of the non-dominated sorting algorithm. We will use the ZDT1 test function with the number design variables $\mathrm{D}=30$ throughout this example, and with the population size $\mathrm{N}=50$.

In [2]:
problem = plat.ZDT1()
D = 30
N = 50

With these variables defined, we will now move onto generating our initial population. We will be using Platypus Solution objects for this, which we will initialise with random problem variables, evaluate, and then append to a list named solutions.

In [3]:
solutions = []

for i in range(N):
    solution = plat.Solution(problem)
    solution.variables = np.random.rand(D)
    solution.evaluate()
    solutions.append(solution)

Let's print out the variables and objectives for the first solution in this list to see what they look like.

In [4]:
print(f"Design variables:\n {solutions[0].variables}")
print(f"Objective values:\n {solutions[0].objectives}")
Design variables:
 [0.14100709 0.77146247 0.99317335 0.24497978 0.29323577 0.13889485
 0.08392614 0.78460659 0.15222921 0.16373011 0.21783006 0.95824192
 0.31035052 0.03555406 0.84666011 0.48333918 0.43851082 0.95912164
 0.77199224 0.20026315 0.54436438 0.28541276 0.94191481 0.87752042
 0.84130242 0.85358657 0.8290964  0.33160395 0.71150884 0.78528698]
Objective values:
 [0.1410070894407488, 5.005306193676968]

Now that we have a population of solutions stored in the solutions variable, we can pass this into the Platypus.nondominated_sort() function that will assign each solution a rank parameter, starting at 0 for the first non-dominated front and incrementing for every additional front that is detected.

In [5]:
plat.nondominated_sort(solutions)

Let's check to see if the first solution in our solutions list has one of these newly calculated rank parameters.

In [6]:
solutions[0].rank
Out[6]:
1

We now have a population of solutions that have been assigned their design variables, their objective values, and their ranks according to the non-dominated sorting algorithm.

Visualising the Non-Dominated Fronts

It would be useful to visualise these non-dominated fronts, according to their newly assigned ranks, using a scatterplot for further investigation. Before we do this, we will migrate our solution variables from their Platypus data-structure to a DataFrame. We will do this for the simple reason that our visualisation tools, e.g. Plotly, work quite well with the DataFrame format.

Let's begin by creating a new DataFrame with the columns f1, f2, and front_rank. f1 and f2 will be used to store our objective values for each solution, which is sufficient for the two-objective ZDT1 problem, and rank_front will be used to store the non-dominated sorting ranks that were calculated earlier.

In [7]:
solutions_df = pd.DataFrame(index=range(N),columns=['f1','f2','front_rank'])
solutions_df.head()
Out[7]:
f1 f2 front_rank
0 NaN NaN NaN
1 NaN NaN NaN
2 NaN NaN NaN
3 NaN NaN NaN
4 NaN NaN NaN

We can see that we've also defined an index range that covers the number of solutions in our population, $\mathrm{N}=50$. This means we have $50$ rows ready, but their values are initialised to NaN (Not A Number), which in this case simply indicates missing data.

Let's use a loop to iterate through our solutions list of $50$ solutions and assign the desired values to the corresponding row in our solutions_df DataFrame

In [8]:
for i in range(N):
    solutions_df.loc[i].f1 = solutions[i].objectives[0]
    solutions_df.loc[i].f2 = solutions[i].objectives[1]
    solutions_df.loc[i].front_rank = solutions[i].rank
        
solutions_df.head()
Out[8]:
f1 f2 front_rank
0 0.141007 5.00531 1
1 0.174235 3.88209 0
2 0.522915 3.86185 2
3 0.31953 4.66639 3
4 0.571165 3.44048 3

We can see our DataFrame now contains the desired values. We're now ready to go with some visualisation, so let's start with a simple scatterplot for the entire population in the objective space.

In [14]:
fig = go.Figure(layout=dict(xaxis=dict(title='f1'), 
                            yaxis=dict(title='f2')))

fig.add_scatter(x=solutions_df.f1, y=solutions_df.f2, mode='markers')

fig.show()

This visualisation gives us an idea of how our solutions look in the objective space. However, we can do better and visualise each front using a different colour. This way, we can visually distinguish between the different non-dominated fronts in our population.

Before producing the visualisation, we may fish to determine the number of fronts found.

In [10]:
solutions_df.front_rank.nunique()
Out[10]:
7

We will then need to produce a sorted vector containing the rank of each front.

In [11]:
fronts = sorted(solutions_df.front_rank.unique())
print(fronts)
[0, 1, 2, 3, 4, 5, 6]

With this vector, we can now visualise each front on the same scatterplot, using different colours to distinguish between their ranking. If you click on items in the legend, you can also disable/enable each front to examine them independently.

In [15]:
fig = go.Figure(layout=dict(
    xaxis=dict(title='f1'),
    yaxis=dict(title='f2')))

for front in fronts:
    front_solutions_df = solutions_df.loc[solutions_df.front_rank == front]
    fig.add_scatter(x=front_solutions_df.f1, y=front_solutions_df.f2, 
                    name=f'front {front}',  
                    mode='markers',
                    marker=dict(color = px.colors.qualitative.Plotly[front]))

fig.show()

The above visualising is great and the code we used to generate it gave us some finer control over the order in which they are plotted, the colours, the axis labels, and so on. However, I wanted to also share with you an approach to get a very similar plot using the plotly.express. As you can see below, we have achieved a very similar plot in only a single line of code.

In [16]:
px.scatter(solutions_df, x="f1", y="f2", color="front_rank")

Conclusion

In this section, we have demonstrated how we can use a software package such as Platypus to calculate the non-dominated sorting ranks for a population, which is a useful technique in the selection stage of an Evolutionary Algorithm. We then used this information to create a scatterplot that distinguishes between each non-dominated front using different colours for every rank. This is useful when visualising a population either during runtime or after a search has concluded.