xuebaunion@vip.163.com

3551 Trousdale Rkwy, University Park, Los Angeles, CA

留学生论文指导和课程辅导

无忧GPA：https://www.essaygpa.com

工作时间：全年无休-早上8点到凌晨3点

微信客服：xiaoxionga100

微信客服：ITCS521

r代写-MATH3609

时间：2021-04-11

MATH3609 Optimization, Networks and Graphs:

Coursework 1 and Coursework 2

Dr Nathan Broomhead and Dr Julian Stander

Module Taught in Semester 2, Academic Year 2020/21

1 General Information about Coursework 1 and

Coursework 2

Please read the following points before attempting the courseworks:

• The deadline for both Coursework 1 and Coursework 2 is 10 am on Thursday, 22nd

April, 2022. You should submit Coursework 1 and Coursework 2 separately through

the MATH3609 Optimization, Networks and Graphs DLE site. Your submission will be

marked anonymously.

• Coursework 1 and Coursework 2 are Group Courseworks. Please work in

a self-assigned group of up to four people. You will need to register your

group on the DLE in the usual way.

Within each coursework, each member of the group will receive the same mark, unless

any member chooses to make use of the Peer Assessment option. You should keep notes

of all your group meetings to use as evidence in case you choose to make use of the Peer

Assessment option. If you wish to make use of the Peer Assessment option, you will need

to make an appointment to see Dr Nathan Broomhead or Dr Julian Stander.

• You should give due consideration to your personal time management to ensure that

your work is submitted in plenty of time prior to the deadline. For further information,

please consult the Student Handbook

https://www.plymouth.ac.uk/your-university/governance/student-handbook,

https://www.plymouth.ac.uk/student-life/your-studies/essential-information/exams/exam-rules-and-

regulations/extenuating-circumstances,

https://www.plymouth.ac.uk/uploads/production/document/path/18/18487/Extenuating_Circumstances_

Policy_and_Procedures_20-21_v.3.pdf and

https://liveplymouthac.sharepoint.com/sites/x70/SitePages/Student-QA.aspx#q-a-session-1-safety-net-

video%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B.

• Each member of the group can expect to spend up to around 15 hours on Coursework 1

and up to around 15 hours on Coursework 2. These suggested times do not include the

preparation work that you will need to do.

• Each coursework counts for 15% = 12 × 30% of your final mark on this module. Marks

will be assigned according to the marking grids on pages 26 and 30.

1

• Marked scripts will be returned within 20 working days of the submission date, and,

in particular, you will get full feedback on your work by Friday, 21st May, 2021, unless

extensions are used.

• You are reminded of the University’s Academic Regulations:

Academic offences occur when activity is undertaken which could confer an unfair

advantage to any candidate(s) in assessment. The University recognises the following

(including any attempt to carry out the actions described) as academic offences, regardless

of intent:

a. Plagiarism, which is copying or paraphrasing of other people’s work or ideas into a

submitted assessment without full acknowledgement. More information on plagiarism

is available here:

https://www.plymouth.ac.uk/student-life/your-studies/essential-information/regulations/plagiarism

b. Collusion, which is unauthorised collaboration of students (or others) in producing a

submitted assessment. The offence of collusion occurs if a student copies any part of

another student’s work, or allows their own work to be copied. Collusion also occurs

if other people contribute significantly to work that a student submits as their own.

The complete list of regulations can be found here:

https://www.plymouth.ac.uk/student-life/your-studies/essential-information/regulations

All group members automatically confirm that they have understood the University’s

policy on plagiarism and collusion when they submit each coursework.

We now state the relevant MATH3609 Optimization, Networks and Graphs Assessed Learning

Outcomes (ALOs) for this assignment.

At the end of the module the learner will be expected to be able to:

ALO1 Develop and analyse methods for the solution of optimisation problems;

ALO3 Demonstrate an understanding of the underlying concepts of graph theory;

ALO4 Use computer software for the solution of optimisation problems and report the

results.

You should keep these ALOs in mind when working on these courseworks.

2

2 Image Restoration: Preparation Work

You should study this section carefully. There has also been a lecture on it.

2.1 The Basic Idea

Let us begin by reading in an image and displaying it.

# Set the working directory

# This is where the image R_image_binary.csv is stored

#

# Please note that your working directory will be different

# from ours (only one of which is needed)

#

## setwd("/Users/julianstander/Dropbox/MATH3609_CW_2020_2021")

#

setwd("/Users/nbroomhead/Dropbox/MATH3609_CW_2020_2021")

#

R_logo <- as.matrix(read.table("R_image_binary.csv",

sep = ",")) # The image is thought of as a matrix

#

image(R_logo,

col = c("black", "white"))

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

Please note that this image is a binary image. It comprises two ‘colours’, black and white. In

Coursework 1 and Coursework 2 we work exclusively with binary images for simplicity.

We will discuss the axes labels later.

3

Unfortunately, often images are received in a degraded or noisy form such as

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

This could be because they are observed and then transmitted by a satellite and so are

subject to atmospheric distortions. Similarly, medical images are often distorted in some way;

just think of an ultrasound image of a baby or even an x-ray or magnetic resonance image

when there has been some patient movement.

4

The aim of image restoration is to try to restore the first image

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

from the degraded image

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

5

2.2 Some Mathematical Notation

Images are represented as matrices. Here is a small example to fix ideas:

x <- matrix(c(0, 0, 0,

0, 0, 1,

1, 1, 1,

0, 0, 1),

nrow = 4,

byrow = TRUE)

#

x

## [,1] [,2] [,3]

## [1,] 0 0 0

## [2,] 0 0 1

## [3,] 1 1 1

## [4,] 0 0 1

#

# Dimension

#

dim(x)

## [1] 4 3

nrow(x)

## [1] 4

ncol(x)

## [1] 3

#

# Specific

#

x[1, 2]

## [1] 0

x[3, 1]

## [1] 1

So the image that is stored in the R object x comprises 4 rows and 3 columns. The value

taken by the image in row 1, column 2 is 0 corresponding to black, while the value taken

by the image in row 3, column 1 is 1 corresponding to white. We say that the above image

comprises 4× 3 = 12 pixels. We shall use i to index the pixels: here i = 1, . . . , n = 12, where n

is the total number of pixels. The value taken at pixel i is denoted xi, where xi can take

the values 0 for black or 1 for white: xi ∈ {0,1}. This enables us to write the whole image as

the vector x = (x1, . . . , xn).

6

We label the pixels column by column:

Columns

R

ow

s

i = 1 or ( 1 , 1 )

i = 2 or ( 2 , 1 )

i = 3 or ( 3 , 1 )

i = 4 or ( 4 , 1 )

i = 5 or ( 1 , 2 )

i = 6 or ( 2 , 2 )

i = 7 or ( 3 , 2 )

i = 8 or ( 4 , 2 )

i = 9 or ( 1 , 3 )

i = 10 or ( 2 , 3 )

i = 11 or ( 3 , 3 )

i = 12 or ( 4 , 3 )

1 2 3

1

2

3

4

So that

x

## [,1] [,2] [,3]

## [1,] 0 0 0

## [2,] 0 0 1

## [3,] 1 1 1

## [4,] 0 0 1

is the vector

as.vector(x)

## [1] 0 0 1 0 0 0 1 0 0 1 1 1

or x = (x1, x2, . . . , x12) = (0,0,1,0,0,0,1,0,0,1,1,1).

7

Please note that in

Columns

R

ow

s

i = 1 or ( 1 , 1 )

i = 2 or ( 2 , 1 )

i = 3 or ( 3 , 1 )

i = 4 or ( 4 , 1 )

i = 5 or ( 1 , 2 )

i = 6 or ( 2 , 2 )

i = 7 or ( 3 , 2 )

i = 8 or ( 4 , 2 )

i = 9 or ( 1 , 3 )

i = 10 or ( 2 , 3 )

i = 11 or ( 3 , 3 )

i = 12 or ( 4 , 3 )

1 2 3

1

2

3

4

we also show the row and column numbers in brackets. We will not usually use these, but

they are added here for information. A point to notice is that the column labels increase

from left to right, while the row labels increase from top to bottom which is different from

standard Cartesian co-ordinates.

8

Let’s consider how R draws an image using its image function

x

## [,1] [,2] [,3]

## [1,] 0 0 0

## [2,] 0 0 1

## [3,] 1 1 1

## [4,] 0 0 1

#

image(x,

col = c("black", "white"))

0.0 0.2 0.4 0.6 0.8 1.0

−

0.

2

0.

2

0.

6

1.

0

#

# Or with better labels

#

n_row <- nrow(x) # Number of rows in x

n_row

## [1] 4

#

n_col <- ncol(x) # Number of columns in x

n_col

## [1] 3

#

image(1:n_row, # Label the rows

1:n_col, # Label the columns

x,

col = c("black", "white"),

xlab = "Row number",

9

ylab = "Column number",

axes = FALSE)

#

axis(1, 1:n_row)

axis(2, 1:n_col)

Row number

Co

lu

m

n

nu

m

be

r

1 2 3 4

1

2

3

So R thinks of the horizontal axis as representing row number, and the vertical axis as

representing column number, just like standard Cartesian co-ordinates: (row i, column j) =(horizontal co-ordinate x,vertical co-ordinate y). This means that the first column of the

matrix x is plotted horizontally at the bottom, the middle column of x is plotted horizontally

in the middle, and the final column of x is plotted horizontally at the top.

These plotting details are not important for what will come. However, we give them here for

information.

10

2.3 A Little More Mathematical Notation

We have seen that in practice we observe a degraded version of a ‘true’ image such as:

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

We refer to this observed or degraded image using the vector y = (y1, . . . , yn), where n is the

number of pixels. Here yi ∈ {0,1}, i = 1, . . . , n.

11

2.4 Some Literature and Other Irrelevant Information

In scientific writing, we should always give credit where credit is due. Many people have

contributed to statistical image restoration. The following are but two references among a

great many:

Besag, J. E. (1986) On the statistical analysis of dirty pictures (with discussion). Journal of

the Royal Statistical Society B, 48, 259–302.

Geman, D. and Geman, S. (1984) Stochastic relaxation, Gibbs distributions and the Bayesian

restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6,

721–41.

Coursework 1 and Coursework 2 are based on the work of

Greig, D. M., Porteous, B. T. and Seheult, A. H. (1989) Maximum a posteriori estimation

for binary images. Journal of the Royal Statistical Society B, 51, 271–279.

Greig, Porteous and Seheult (1989) assume black is coded as 1 and white as 0. In Coursework

1 and Coursework 2, we use take 1 to represent white and 0 to represent black, as this is

more natural and usual.

The work of Greig, Porteous and Seheult (1989) was expanded on in these works:

Jubb, M. D. (1989) Image Reconstruction. PhD thesis. University of Bath.

Stander, J. (1992) Some Topics in Statistical Image Analysis. PhD thesis. University of Bath.

Stander, J. and Silverman, B. W. (1994) Temperature schedules for simulated annealing.

Statistics and Computing, 4, 21–32.

Franconi, L. and Jennison, C. (1997) Comparison of a genetric algorithm and simulated

annealing in an application to statistical image reconstruction. Statistics and Computing, 7,

193–207.

No knowledge of any of this literature is assumed in either Coursework 1 or

Coursework 2.

We end this section by remarking that images are not usually binary. Often each pixel can

take a grey level so that xi ∈ {0,1, . . . ,255} where 0 is black and 255 is white for example,

or a colour, which is more complicated to represent. We work with binary images to keep

life relatively simple, as this may be the first time that you have met image reconstruction

problems.

12

2.5 Statistical Models

Image restoration proceeds by specifying statistical models for possible true images x =(x1, . . . , xn) and for degraded images y = (y1, . . . , yn).

Please note that this is not a statistics module, so not all the details that we present here are

needed for you to be able to do Coursework 1 or Coursework 2. We include these details

to explain how the image restoration problem can be expressed as an optimization

problem. It is with that optimization problem that we shall work.

2.6 A Statistical Model for y given x

We will assume the following statistical model for yi given xi:

Pr(yi = 0 ∣xi = 0) = 1 − p (no error made, no swap)

Pr(yi = 1 ∣xi = 0) = p (error made, swap from 0 to 1)

Pr(yi = 0 ∣xi = 1) = p (error made, swap from 1 to 0)

Pr(yi = 1 ∣xi = 1) = 1 − p (no error made, no swap)

Essentially, we are assuming that the colour at pixel i is swapped from black to white, or

from white to black, with probability p ∈ [0,1]; otherwise, it is not changed. For simplicity,

we assume that p is known.

We assume that the error mechanism at each pixel is independent of the error mechanisms at

every other pixel. This enables us to write down Pr(y ∣x):

Pr(y ∣x) = n∏

i=1 Pr(yi ∣xi).

Let’s recall that the degraded image y is observed, while the ‘true’ image x is unknown and

needs to be estimated. Thought of as a function of the unknown image vector x, Pr(y ∣x)

is the likelihood function that you may have met in MATH2602 Statistical Inference and

Regression or MATH3613 Data Modelling. We shall write L(y ∣x) = Pr(y ∣x).

13

Now let’s write f(yi ∣xi) = Pr(yi ∣xi), so that f(1 ∣0) = Pr(yi = 1 ∣xi = 0) for example. We

have

L(y ∣x) = Pr(y ∣x)

= n∏

i=1 Pr(yi ∣xi)

= n∏

i=1 f(yi ∣1)xif(yi ∣0)1−xi ,

since f(yi ∣1)xif(yi ∣0)1−xi = f(yi ∣1) if xi = 1 and f(yi ∣0) if xi = 0.

From this it follows that

L(y ∣x) = n∏

i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}

xi

. (1)

Later, we will use λi, which is defined as follows:

λi = log{f(yi ∣1)

f(yi ∣0)} . (2)

Please note that λi only depends on the observed yi and not on xi. Please also note that, if

yi = 0, then

λi = log{f(yi ∣1)

f(yi ∣0)} = log{f(0 ∣1)f(0 ∣0)} = log ( p1 − p) ,

while, if yi = 1, then

λi = log{f(yi ∣1)

f(yi ∣0)} = log{f(1 ∣1)f(1 ∣0)} = log (1 − pp ) .

14

2.7 A Statistical Model for x

To define a statistical model for x, we need to explain what it means for two pixels to be

neighbours. This can be done by the following picture

Columns

R

ow

s

i = 1 or ( 1 , 1 )

i = 2 or ( 2 , 1 )

i = 3 or ( 3 , 1 )

i = 4 or ( 4 , 1 )

i = 5 or ( 1 , 2 )

i = 6 or ( 2 , 2 )

i = 7 or ( 3 , 2 )

i = 8 or ( 4 , 2 )

i = 9 or ( 1 , 3 )

i = 10 or ( 2 , 3 )

i = 11 or ( 3 , 3 )

i = 12 or ( 4 , 3 )

1 2 3

1

2

3

4

The red and blue lines indicate neighbouring pixles; we use blue to indicate horizontal

neighbours and red to indicate vertical neighbours, although this distinction is not important.

As an example we see that the neighbours of pixel i = 6 are 5, 10, 7 and 2. It turns out

that it simpler to assume that the image is wrapped on a torus, so that, for example the

neighbours of pixel i = 1 are 4 (wrapped top to bottom), 5, 2 and 9 (wrapped left to right).

This assumption has very little practical effect for large images, but it does mean that every

pixel has four neighbours.

15

These neighbourhood relationships lead to a symmetric adjacency matrix A, which here

will be of dimension n × n = 12 × 12:

## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]

## [1,] 0 1 0 1 1 0 0 0 1 0 0 0

## [2,] 1 0 1 0 0 1 0 0 0 1 0 0

## [3,] 0 1 0 1 0 0 1 0 0 0 1 0

## [4,] 1 0 1 0 0 0 0 1 0 0 0 1

## [5,] 1 0 0 0 0 1 0 1 1 0 0 0

## [6,] 0 1 0 0 1 0 1 0 0 1 0 0

## [7,] 0 0 1 0 0 1 0 1 0 0 1 0

## [8,] 0 0 0 1 1 0 1 0 0 0 0 1

## [9,] 1 0 0 0 1 0 0 0 0 1 0 1

## [10,] 0 1 0 0 0 1 0 0 1 0 1 0

## [11,] 0 0 1 0 0 0 1 0 0 1 0 1

## [12,] 0 0 0 1 0 0 0 1 1 0 1 0

Do you understand why?

16

We can now define a probability mass function on images x as

Pr(x) = exp (−β × number of discrepant neighbouring pairs in x)∑images x′ exp (−β × number of discrepant neighbouring pairs in x′) ,

in which β ≥ 0 is a constant that we will discuss later. The sum over all images x′ comprises

2n terms, where n is the number of pixels, and can be very hard to compute when n is large.

As this sum does not depend on the argument x of Pr(x), we can discard it, if we replace

the equality with a ‘proportional to’ sign. So we write

Pr(x)∝ exp (−β × number of discrepant neighbouring pairs in x) .

For the image x that we have been working with:

## [,1] [,2] [,3]

## [1,] 0 0 0

## [2,] 0 0 1

## [3,] 1 1 1

## [4,] 0 0 1

the number of discrepant neighbouring pairs is 10:

Columns

R

ow

s

0

0

1

0

0

0

1

0

0

1

1

1

1 2 3

1

2

3

4

Hence, the probability of this image is Pr(x) = C exp(−β × 10), where C is a constant.

17

Similarly, the number of discrepant neighbouring pixels in the all black and the all white

images is, unsurprisingly, 0:

Columns

R

ow

s

0

0

0

0

0

0

0

0

0

0

0

0

1 2 3

1

2

3

4

Columns

R

ow

s

1

1

1

1

1

1

1

1

1

1

1

1

1 2 3

1

2

3

4

This means that the all black and the all white images are the most probable with probabilities

C exp(β × 0) = C each.

18

The chess-board image is the least probable, as it has the highest number of discrepant

neighbouring pixels:

Columns

R

ow

s

1

0

1

0

0

1

0

1

1

0

1

0

1 2 3

1

2

3

4

There are 20 discrepant pairs, meaning that the probability of the chess board image is

C exp(−β × 20).

You may be asking: what does β ≥ 0 do? The larger the value of β, the higher the relative

probabity given to the all black and to the all white image. In other words, higher β values

encourage ‘smoother’ images. Because of this, sometimes β is referred to as a smoothing

parameter.

Finally, we remark that our probability model for images x,

Pr(x)∝ exp (−β × number of discrepant neighbouring pairs in x) ,

is essentially the same as the Ising model used a lot in physics and statistical mechanics

to represent substances such as magnets. To quote Federico Camia: ‘The Ising model,

introduced by Lenz in 1920 to describe ferromagnetism, is one of the most studied models

of statistical mechanics. Its two dimensional version has played a special role in rigorous

statistical mechanics since Peierls’ proof of a phase transition in 1936 and Onsager’s derivation

of the free energy in 1944.’

19

It can be shown that the

number of discrepant neighbouring pairs in x = 12 n∑i=1 n∑j=1aij (xi − xj)2 ,

in which x = (x1, . . . , xn) is an image and aij is the (i, j)-th element of the adjacence matrix

A such that

aij = {1 if pixels i and j are neighbours0 if pixels i and j are not neighbours.

A consequence of this result is that

Pr(x) ∝ exp (−β × number of discrepant neighbouring pairs in x)

= exp(−β × 12 n∑i=1 n∑j=1aij (xi − xj)2)

= exp(−β2 n∑i=1 n∑j=1aij (xi − xj)2) . (3)

Another irrelevant fact is that more general probability mass functions can be defined by

replacing (xi − xj)2 with φ (xi − xj) where the function φ satisfies certain properties.

20

2.8 A Statistical Model for x given y

It makes sense to consider Pr(x ∣ y), as we are given the observed image y. Pr(x ∣ y) can be

found using Bayes’ theorem as follows:

Pr(x ∣ y) ∝ L(y ∣x) ×Pr(x)

∝ n∏

i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}

xi ×Pr(x) from (1)

∝ n∏

i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}

xi × exp(−β2 n∑i=1 n∑j=1aij (xi − xj)2) from (3) (4)

In the language of Bayesian statistical inference, Pr(x) is the prior probability mass

function of all images x before observing the degraded image y. Pr(x) expresses our belief

about x before seeing the data in terms of probabilities. Pr(x ∣ y) is our belief about x after

seeing the data y and is referred to as our posterior probability mass function.

Please note that, as mentioned, it does not matter if you do not follow all these details,

although they are not hard. The important part is that image restoration leads to an

optimization problem, as we shall see.

21

2.9 The Restoration of a Binary Image as an Optimization Problem

The restored image is defined as the most probable image under the posterior Pr(x ∣ y).

That is, it is the image x that maximizes Pr(x ∣ y). For this reason it is referred to as the

maximum a posteriori or MAP image.

From (4) we see that the restored image maximizes

Pr(x ∣ y)∝ n∏

i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}

xi × exp(−β2 n∑i=1 n∑j=1aij (xi − xj)2)

over all images x.

We know that it’s very often easier to work with the logarithm of a function, rather than

the function itself, when maximizing. This is because the logarithm of a function and the

function itself are maximized at the same point, and also because taking logarithms turns

products into sums, which are easier to deal with. So, the restored image maximizes

log {Pr(x ∣ y)} = additive constant D + n∑

i=1 log [f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}

xi] − β2 n∑i=1 n∑j=1aij (xi − xj)2

= D + n∑

i=1 log {f(yi ∣0)} + n∑i=1 log [{f(yi ∣1)f(yi ∣0)}

xi] − β2 n∑i=1 n∑j=1aij (xi − xj)2

= D + n∑

i=1 log {f(yi ∣0)} + n∑i=1 xi log{f(yi ∣1)f(yi ∣0)} − β2 n∑i=1 n∑j=1aij (xi − xj)2

= D + n∑

i=1 log {f(yi ∣0)} + n∑i=1 xiλi − β2 n∑i=1 n∑j=1aij (xi − xj)2 (5)

by the definition of λi in (2).

Please note that ∑ni=1 log {f(yi ∣0)} does not depend on x, and so can be considered a constant.

It therefore follows that we need to maximize

n∑

i=1 xiλi − β2 n∑i=1 n∑j=1aij (xi − xj)2

over x = (x1, . . . , xn), or, equivalently, to minimize

− n∑

i=1 xiλi + β2 n∑i=1 n∑j=1aij (xi − xj)2

over all images x.

22

So, we have finally arrived at our image restoration optimization problem. Given a degraded

image y = (y1, . . . , yn), we can find our restored image by solving the problem:

minimize: − n∑

i=1 xiλi + β2 n∑i=1 n∑j=1aij (xi − xj)2 over images x = (x1, . . . , xn). (6)

It can be shown that the minimiation problem given by (6) is equivalent to the problem

minimize:

n∑

i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2 (7)

over images x, in which max(a, b) means the maximum of a and b, so that, for example,

max(0,1) = 1 and max(0,−1) = 0.

This quantity is never negative.

23

2.10 Solving the Minimization Problem using the Ford-Fulkerson

Algorithm

We will now see how Greig, Porteous and Seheult (1989) solved the minimization problem

(7), namely,

minimize:

n∑

i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2

over images x, using the Ford-Fulkerson algorithm. To do this, Greig, Porteous and Seheult

(1989) defined a network with capacities as follows:

2.10.1 Definition of the Network with Capacities

• There are n + 2 vertices:

– A source s;

– A sink t;

– A vertex for each of the n pixels, indexed by i or j below.

• If λi > 0, there is an arc (directed edge) (s, i) from the source s to i with capactiy

csi = λi; otherwise

• If λi < 0, there is an arc (i, t) from i to the sink t with capactiy cit = −λi;

• There is an undirected edge {i, j} between two internal vertices i and j (that corre-

spond to pixels i and j) with capacities cij = β aij. This can be thought of as two arcs,

one between i and j, and the other between j and i, both of which have capacities β aij .

• There are no other arcs.

24

2.10.2 Definition of a Cut

Let us say that x = (x1, . . . , xn) is an image that is being proposed as a reconstruction of the

observed y = (y1, . . . , yn).

• Partitioning the Network

We can use x to partition the network into two sets:

W (x) = {s} ∪ {i ∶ xi = 1}, source s and all white pixels

B(x) = {i ∶ xi = 0} ∪ {t}, all black pixels and the sink t.

This partition of the vertices in the network into the sets W (x) and B(x) defines a cut.

• The Cut Capacity

The associated cut capacity is

C(x) ∶= C(∂+(W (x))) = ∑

kl∈∂+(W (x)) ckl. (8)

Notice that the arcs of ∂+(W (x)) are of three types: arcs starting at s; arcs ending at t; arcs

between intermediate vertices.

The cut capacity (8) can also be written as

C(x) = ∑

k∈W (x) ∑l∈B(x) ckl.

It turns out that

C(x) = n∑

i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2 ,

which is exactly the function that we wanted to minimize in (7) to restore the image from y.,

So, the image restoration problem becomes a max-flow/min-cut problem: we find the cut

that minimizes the cut capacity C(x) and this determines x. We can do this using the

Ford-Fulkerson algorithm (or modern improvements of it). ,

We will expand on and explore the above in Coursework 1 and Coursework 2.

25

3 Coursework 1

Marks will be assigned for Coursework 1 according to this marking grid.

MATH3609 Optimization, Networks and Graphs:

Coursework 1 Marking Grid

Assessment Area Maximum Awarded Feedback

Mark Mark

Task 1: Correctness and va-

lidity of argument. Read-

ability of explanation.

15

Task 2: Correctness and va-

lidity of argument. Read-

ability of explanation.

20

Task 3: Correctness and va-

lidity of argument. Read-

ability of explanation.

10

Task4: Correctness and va-

lidity of argument. Read-

ability of explanation.

25

Task 5: Correctness and va-

lidity of argument. Read-

ability of explanation.

15

Task 6: Correctness and va-

lidity of argument. Read-

ability of explanation.

15

Total1 100

1This mark is provisional. Like all marks it is considered by an External Examiner and an Assessment Panel.

26

3.1 Coursework 1 Tasks

Coursework 1 comprises a number of tasks. You need to produce a document containing

your solutions to these tasks. These may be hand-written and then scanned or typed up

using, for example, LATEX. It is your responsibility to make sure that they are readable.

The page limit is fifteen pages. Please note that this is the limit and should not be

considered as the target. Well written, clear and concise mathematics will receive high credit.

Please arrange your report under subheadings referring to the tasks, that is, Task 1,

Task 2, . . .. This will help you to respond systematically to all parts of Coursework 1.

The marks for each task are shown in the grid on page 26.

Task 1: Let G = (V,E) be a graph with vertex set V = {1,2, . . . , n} and let A be the

adjacency matrix of G with (i, j)-th entry aij. Suppose x∶V → {0,1} is a function which

assigns to each vertex i a number xi, which is either 0 or 1. A pair of neighbouring vertices i

and j are called discrepant if xi ≠ xj.

Explain why the number of discrepant neighbouring (unordered) pairs is given by the formula:

1

2

n∑

i=1

n∑

j=1aij (xi − xj)2 .

Task 2: Explain why the minimization problem given by (6) on page 23 is equivalent to the

problem

minimize:

n∑

i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2

over images x, in which max(a, b) means the maximum of a and b, so that, for example,

max(0,1) = 1 and max(0,−1) = 0.

Task 3: Prove that

n∑

i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2

is never negative.

Task 4: Prove that

C(∂+(W (x))) = n∑

i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2 .

Hint: The arcs of ∂+(W (x)) can be split into three types, namely, arcs that start at s, arcs

that end at t and arcs between intermediate vertices. Consider the sums of the capacities

arcs of each type separately.

27

Task 5: Show that every cut K of the network can be written as

K = ∂+(W (x))

for some image x, and deduce using this and the previous tasks that the Ford-Fulkerson

algorithm can be used to find a restored image.

Task 6: Is there always a unique solution x to the minimization problem (i.e. a unique

restored image)? Illustrate your answer with an example.

28

3.2 Coursework 1: What You Need to Submit

Your group needs to submit the following file electronically using the DLE:

• A Portable Document Format file MATH3609_Coursework_1_2020_2021.pdf containing

your report for Coursework 1.

All members of the group should take responsibility for seeing that your work

is submitted before the deadline.

It is also the responsibility of all members to check that the correct file has

been uploaded and that it is readable.

If anything is unclear, you should ask without delay.

29

4 Coursework 2

Marks will be assigned for Coursework 2 according to this marking grid.

MATH3609 Optimization, Networks and Graphs:

Coursework 2 Marking Grid

Assessment Area Maximum Awarded Feedback

Mark Mark

Image Restoration

Tasks: Defining the set-

up. Correctness, clarity

and presentation of code.

25

Image Restoration

Tasks: Optimizing,

presenting the results,

choosing β. Depth of

understanding and quality

of description of insights

gained.

50

Report: display results

clearly, concisely, system-

atically and professionally,

using correct spelling and

grammar. This is important

for employability.

25

Total2 100

2This mark is provisional. Like all marks it is considered by an External Examiner and an Assessment Panel.

30

4.1 Coursework 2 Tasks

Coursework 2 comprises a number of tasks. You need to produce a report of your work

following the instructions below.

Your report should contain well presented and annotated R code.

The page limit is twenty pages, including R code and figures. Please see page 39

for more details. Please do not submit an additional appendix as it will not be considered.

Reports that contain irrelevant or uninteresting discussion or R code will be penalized. Please

note that twenty pages is the limit and should not be considered as the target. Well written,

concise reports will receive high credit.

You may also omit lengthy, technical R code, if you feel that it does not add to

the argument that you are making. Essential code should, however, be shown.

Please arrange your report under subheadings referring to the tasks, that is, Task 1,

Task 2, . . .. This will help you to respond systematically to all parts of Coursework 2.

31

4.1.1 Implementing the Ford-Fulkerson Algorithm to Restore the

Image in R: Getting Started

Following on from the information in the section ‘Image Restoration: Preparation Work’ we

will now produce an image restoration from observed data y using the Ford-Fulkerson

algorithm implemented in R.

First, we need to read y into R. The data that we will use are supplied in the file

R_degraded_image_0.2.csv. We can read in and display these data in the same was

as we did at the beginning of this document. Here is the code again to make sure that you

get started properly:

y <- as.matrix(read.table("R_degraded_image_0.2.csv",

sep = ","))

#

image(y,

col = c("black", "white"))

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

#

n_row <- nrow(y)

n_col <- ncol(y)

#

c(n_row, n_col)

## [1] 23 28

#

y_vec <- as.vector(y)

#

length(y_vec)

## [1] 644

32

We can label the axes in a possibly nicer way:

image(1:n_row,

1:n_col,

y,

xlab = "Row number",

ylab = "Column number",

main = "Data: degraded image y",

col = c("black", "white"))

5 10 15 20

5

10

15

20

25

Data: degraded image y

Row number

Co

lu

m

n

nu

m

be

r

• It is not hard to imagine that y is a degraded version of the image in R_logo seen above.

• It is known that p = 0.2. (In practice, p and any other unknown parameter or parameters

can also be estimated form the data by assigning prior distributions to them and sampling

from the resulting posterior using Markov chain Monte Carlo methods as discussed in

MATH3613 Data Modelling. This is beyond the scope of Coursework 2.)

• We will assume for now that β = 1.

33

4.1.2 Defining the Vector of λ Values

Task 1:

Using the function ifelse or otherwise, define the vector of λ values (λ1, . . . , λn) in R and

show a few λ values. Identify and record in R which λi values are positive and which λi values

are negative using the which function, storing the results in vectors which_lambda_positive

and which_lambda_negative respectively.

4.1.3 Defining the Network with Capacities

We need to define a data frame comprising three variables or columns called from, to and

capacity. This data frame defines the capacities in the network. If, for example, we wish to

specify ckl, then the variable from would take the value k, the variable to would take the

value l and the variable capacity would take the value ckl.

One way to define the required data frame of capacities is to construct three data

frames df_source, df_sink and df_internal, and to combine df_source, df_sink and

df_internal by binding their rows into one data frame df using rbind.

The data frame df_source records the capacities of arcs coming from the source s.

Task 2: Define df_source in R. The square brackets [] may help you here. In order to

check this, the first few rows of the result are given

## from to capacity

## 1 source 2 1.386294

## 2 source 3 1.386294

## 3 source 4 1.386294

## 4 source 6 1.386294

## 5 source 7 1.386294

## 6 source 8 1.386294

The data frame df_sink records the capacities of arcs going to the sink t.

Task 3: Define df_sink in R and show a part of it.

34

Before we define df_internal, we need to obtain the adjacency matrix A. This can be

done by using the function get_neighbourhood_matrix with which you are supplied in the

file of helper functions helper_functions.R; you do not need to understand the details of

these helper functions. Please note that get_neighbourhood_matrix is defined using other

functions, so that you’ll have to make sure that all the functions in the file of helper functions

helper_functions.R are loaded into R. We also need to extract the rows and columns of A

corresonding to neighbouring pixels:

A <- get_neighbourhood_matrix(y, n_row)

#

A[1:5, 1:5]

## [,1] [,2] [,3] [,4] [,5]

## [1,] 0 1 0 0 0

## [2,] 1 0 1 0 0

## [3,] 0 1 0 1 0

## [4,] 0 0 1 0 1

## [5,] 0 0 0 1 0

#

neighbouring_pixels <- which(A == 1,

arr.ind = TRUE)

#

head(neighbouring_pixels)

## row col

## [1,] 2 1

## [2,] 23 1

## [3,] 24 1

## [4,] 622 1

## [5,] 1 2

## [6,] 3 2

35

The data frame df_internal records cij where i and j correspond to pixels i and j.

Task 4: Define df_internal in R and show a part of it. Define the data frame df using

rbind. Use the function graph_from_data_frame from the igraph package to create an

igraph graph called g from df. Produce a plot of your capacitated network, as shown here;

please note, that it’s very big (but rather nice)! Please note that your plot may be displayed in

a different way from the one here. You may find https://igraph.org/r/doc/plot.common.html

helpful.

source

2

23

24

622

1

3

25

623

4

26

624

5

27

625

6

28

626

7

29

627

8

30

628

9

31

629

10

32

630

11

33

631

12

34

632

13

35633

14 36

634

15

37

635

16

38

636

17

39 637

18

40

638

19

41

639

20

42

640

21

43

64122

44

642

45

643

46

64447

48

49

50

51

52

53

54

55

56

57

58

5960

6162

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

8283

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124 125

126

127

128

129

130

131132

133

134

135

136

137

138

139

140

141

142

143

144

145

146 147

148

149

150

151

152153

154

155

156

157

158

159

160

161 162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195196

197198

199

200

201

202

203

204

205

206

207 208

209

210

211

212

213

214

215

216

217218219220

221

222

223

224

225

226

227

228

229

230 231

232 233

234

235

236

237

238

239

240

241242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264265

266

267

268

269

270

271

272

273

274 275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300 301

302

303

304

305

306

307

308

309

310

311

312313314

315

316

317

318

319

320

321 322

323

324

325

326

327

328

329330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384385

386

387

388

389

390

391

392

393

394

395

396

397

398

399400

401

402

40

404

405

406

407

408409

410

411

412

413

414

415

416

417

418

419

420

421

422423

424

425

426

427

428

429

430

431

432

433

434

435

436 437

438

439

440

441

442

443

444

445446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

47

472

473

474

475

476

477

478

79

480

481

482

483

484

485486487

488

489

490

491 492

493

494

495

496

497498

499

500

501

502

503

504

505

506

507

508

509

510

511

512513

514

515

516

517

518

519

520521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

540

541

542543

544

545

546

547

548

549550

551

552553

554

555

556

557

558

559

560

561

562

563

564

565566567

568

569

570

571

572573574

575

576

577

578

579

580

581

582

583

584

585

586 587

588589590

591

592

593

594

595

596

597

598

599

600601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

sink

36

4.1.4 Producing the Image Restoration by Maximizing Flow or Minimizing Cut Capacity

The function max_flow provided by the igraph package finds the maximum flow in a graph

or equivalently the minimum cut capacity. max_flow effectively implements an enhanced

modern version of the Ford-Fulkerson algorithm. Amongst other output, max_flow

returns a partition of the vertices into two sets that it calls partition1 and partition2.

partition1 is the set W (x) and partition2 is the set B(x) defined above, corresponding to

the ‘maximum flow’/‘minimum cut’ x. partition1 and partition2 are quite complicated

R objects. However, the part of interest to us can be extracted using the names function;

a character vector is given by names. Once names has been applied and any unnecessary

element has been removed, the resulting character vector can be turned into pixels numbers

using the as.numeric function.

Task 5: Use max_flow to maximize the flow in g. Report the value of the maximum flow.

Working with partition1 find the numerical values of the white pixels in the set W (x),

corresponding to the max-flow/min-cut, and report the first six of these. Similarly, find the

numerical values of the black pixels in the set B(x) and report the first six.

Task 6: Use these pixel values to define the reconstructed image. An easy way (but not the

only way) to do this is to use rep to define a vector, x_map_vec say, of the same length as

y_vec comprising NAs, and then to assign the value 0 to pixels that have been identified as

black and the value 1 to pixels that have been identified as white. Then turn the vector into

a matrix, remembering that pixels are labelled column by column; please see the help file for

matrix, if necessary. Display your image. Here is the one that we obtained.

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

37

4.1.5 Choosing β

This reconstruction x obtained from the observed image y was obtained using β = 1. In fact,

the choice of β is in the hands of the user.

Task 7:

Discuss and illustrate how you could choose β.

There is no right answer here. What we are looking for is sensible suggestions for

the choice of β, discussed clearly and critically. Your discussion should be illus-

trated with useful figures, as appropriate. Concise answers comprising sensible

and well illustrated suggestions will receive high credit.

38

4.2 Coursework 2: Report Production

The ability to report technical results clearly and concisely is an important skill.

You should write a report that, as a minimum:

• discusses in detail the above tasks and

• contains well presented R code developed in responding to the tasks.

You are advised, but not obliged, to use R Markdown. We will certainly help you to get

started with R Markdown, if you ask. You can use LATEX from within R Markdown, but you

are under no obligation to use LATEX.

You should use clearly readable fonts and maintain sensible margins and line spacing through-

out. Your work should not look squashed. The default settings of R Markdown are completely

acceptable.

4.3 Coursework 2: What You Need to Submit

Your group needs to submit the following file electronically using the DLE:

• A Portable Document Format file MATH3609_Coursework_2_2020_2021.pdf containing

your report for Coursework 2.

All members of the group should take responsibility for seeing that your work

is submitted before the deadline.

It is also the responsibility of all members to check that the correct file has

been uploaded.

There have been occasions when R code has been submitted instead of the

required report. This error is easy to make, but will have severe consequences.

If anything is unclear, you should ask without delay.

39

学霸联盟

Coursework 1 and Coursework 2

Dr Nathan Broomhead and Dr Julian Stander

Module Taught in Semester 2, Academic Year 2020/21

1 General Information about Coursework 1 and

Coursework 2

Please read the following points before attempting the courseworks:

• The deadline for both Coursework 1 and Coursework 2 is 10 am on Thursday, 22nd

April, 2022. You should submit Coursework 1 and Coursework 2 separately through

the MATH3609 Optimization, Networks and Graphs DLE site. Your submission will be

marked anonymously.

• Coursework 1 and Coursework 2 are Group Courseworks. Please work in

a self-assigned group of up to four people. You will need to register your

group on the DLE in the usual way.

Within each coursework, each member of the group will receive the same mark, unless

any member chooses to make use of the Peer Assessment option. You should keep notes

of all your group meetings to use as evidence in case you choose to make use of the Peer

Assessment option. If you wish to make use of the Peer Assessment option, you will need

to make an appointment to see Dr Nathan Broomhead or Dr Julian Stander.

• You should give due consideration to your personal time management to ensure that

your work is submitted in plenty of time prior to the deadline. For further information,

please consult the Student Handbook

https://www.plymouth.ac.uk/your-university/governance/student-handbook,

https://www.plymouth.ac.uk/student-life/your-studies/essential-information/exams/exam-rules-and-

regulations/extenuating-circumstances,

https://www.plymouth.ac.uk/uploads/production/document/path/18/18487/Extenuating_Circumstances_

Policy_and_Procedures_20-21_v.3.pdf and

https://liveplymouthac.sharepoint.com/sites/x70/SitePages/Student-QA.aspx#q-a-session-1-safety-net-

video%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B%E2%80%8B.

• Each member of the group can expect to spend up to around 15 hours on Coursework 1

and up to around 15 hours on Coursework 2. These suggested times do not include the

preparation work that you will need to do.

• Each coursework counts for 15% = 12 × 30% of your final mark on this module. Marks

will be assigned according to the marking grids on pages 26 and 30.

1

• Marked scripts will be returned within 20 working days of the submission date, and,

in particular, you will get full feedback on your work by Friday, 21st May, 2021, unless

extensions are used.

• You are reminded of the University’s Academic Regulations:

Academic offences occur when activity is undertaken which could confer an unfair

advantage to any candidate(s) in assessment. The University recognises the following

(including any attempt to carry out the actions described) as academic offences, regardless

of intent:

a. Plagiarism, which is copying or paraphrasing of other people’s work or ideas into a

submitted assessment without full acknowledgement. More information on plagiarism

is available here:

https://www.plymouth.ac.uk/student-life/your-studies/essential-information/regulations/plagiarism

b. Collusion, which is unauthorised collaboration of students (or others) in producing a

submitted assessment. The offence of collusion occurs if a student copies any part of

another student’s work, or allows their own work to be copied. Collusion also occurs

if other people contribute significantly to work that a student submits as their own.

The complete list of regulations can be found here:

https://www.plymouth.ac.uk/student-life/your-studies/essential-information/regulations

All group members automatically confirm that they have understood the University’s

policy on plagiarism and collusion when they submit each coursework.

We now state the relevant MATH3609 Optimization, Networks and Graphs Assessed Learning

Outcomes (ALOs) for this assignment.

At the end of the module the learner will be expected to be able to:

ALO1 Develop and analyse methods for the solution of optimisation problems;

ALO3 Demonstrate an understanding of the underlying concepts of graph theory;

ALO4 Use computer software for the solution of optimisation problems and report the

results.

You should keep these ALOs in mind when working on these courseworks.

2

2 Image Restoration: Preparation Work

You should study this section carefully. There has also been a lecture on it.

2.1 The Basic Idea

Let us begin by reading in an image and displaying it.

# Set the working directory

# This is where the image R_image_binary.csv is stored

#

# Please note that your working directory will be different

# from ours (only one of which is needed)

#

## setwd("/Users/julianstander/Dropbox/MATH3609_CW_2020_2021")

#

setwd("/Users/nbroomhead/Dropbox/MATH3609_CW_2020_2021")

#

R_logo <- as.matrix(read.table("R_image_binary.csv",

sep = ",")) # The image is thought of as a matrix

#

image(R_logo,

col = c("black", "white"))

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

Please note that this image is a binary image. It comprises two ‘colours’, black and white. In

Coursework 1 and Coursework 2 we work exclusively with binary images for simplicity.

We will discuss the axes labels later.

3

Unfortunately, often images are received in a degraded or noisy form such as

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

This could be because they are observed and then transmitted by a satellite and so are

subject to atmospheric distortions. Similarly, medical images are often distorted in some way;

just think of an ultrasound image of a baby or even an x-ray or magnetic resonance image

when there has been some patient movement.

4

The aim of image restoration is to try to restore the first image

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

from the degraded image

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

5

2.2 Some Mathematical Notation

Images are represented as matrices. Here is a small example to fix ideas:

x <- matrix(c(0, 0, 0,

0, 0, 1,

1, 1, 1,

0, 0, 1),

nrow = 4,

byrow = TRUE)

#

x

## [,1] [,2] [,3]

## [1,] 0 0 0

## [2,] 0 0 1

## [3,] 1 1 1

## [4,] 0 0 1

#

# Dimension

#

dim(x)

## [1] 4 3

nrow(x)

## [1] 4

ncol(x)

## [1] 3

#

# Specific

#

x[1, 2]

## [1] 0

x[3, 1]

## [1] 1

So the image that is stored in the R object x comprises 4 rows and 3 columns. The value

taken by the image in row 1, column 2 is 0 corresponding to black, while the value taken

by the image in row 3, column 1 is 1 corresponding to white. We say that the above image

comprises 4× 3 = 12 pixels. We shall use i to index the pixels: here i = 1, . . . , n = 12, where n

is the total number of pixels. The value taken at pixel i is denoted xi, where xi can take

the values 0 for black or 1 for white: xi ∈ {0,1}. This enables us to write the whole image as

the vector x = (x1, . . . , xn).

6

We label the pixels column by column:

Columns

R

ow

s

i = 1 or ( 1 , 1 )

i = 2 or ( 2 , 1 )

i = 3 or ( 3 , 1 )

i = 4 or ( 4 , 1 )

i = 5 or ( 1 , 2 )

i = 6 or ( 2 , 2 )

i = 7 or ( 3 , 2 )

i = 8 or ( 4 , 2 )

i = 9 or ( 1 , 3 )

i = 10 or ( 2 , 3 )

i = 11 or ( 3 , 3 )

i = 12 or ( 4 , 3 )

1 2 3

1

2

3

4

So that

x

## [,1] [,2] [,3]

## [1,] 0 0 0

## [2,] 0 0 1

## [3,] 1 1 1

## [4,] 0 0 1

is the vector

as.vector(x)

## [1] 0 0 1 0 0 0 1 0 0 1 1 1

or x = (x1, x2, . . . , x12) = (0,0,1,0,0,0,1,0,0,1,1,1).

7

Please note that in

Columns

R

ow

s

i = 1 or ( 1 , 1 )

i = 2 or ( 2 , 1 )

i = 3 or ( 3 , 1 )

i = 4 or ( 4 , 1 )

i = 5 or ( 1 , 2 )

i = 6 or ( 2 , 2 )

i = 7 or ( 3 , 2 )

i = 8 or ( 4 , 2 )

i = 9 or ( 1 , 3 )

i = 10 or ( 2 , 3 )

i = 11 or ( 3 , 3 )

i = 12 or ( 4 , 3 )

1 2 3

1

2

3

4

we also show the row and column numbers in brackets. We will not usually use these, but

they are added here for information. A point to notice is that the column labels increase

from left to right, while the row labels increase from top to bottom which is different from

standard Cartesian co-ordinates.

8

Let’s consider how R draws an image using its image function

x

## [,1] [,2] [,3]

## [1,] 0 0 0

## [2,] 0 0 1

## [3,] 1 1 1

## [4,] 0 0 1

#

image(x,

col = c("black", "white"))

0.0 0.2 0.4 0.6 0.8 1.0

−

0.

2

0.

2

0.

6

1.

0

#

# Or with better labels

#

n_row <- nrow(x) # Number of rows in x

n_row

## [1] 4

#

n_col <- ncol(x) # Number of columns in x

n_col

## [1] 3

#

image(1:n_row, # Label the rows

1:n_col, # Label the columns

x,

col = c("black", "white"),

xlab = "Row number",

9

ylab = "Column number",

axes = FALSE)

#

axis(1, 1:n_row)

axis(2, 1:n_col)

Row number

Co

lu

m

n

nu

m

be

r

1 2 3 4

1

2

3

So R thinks of the horizontal axis as representing row number, and the vertical axis as

representing column number, just like standard Cartesian co-ordinates: (row i, column j) =(horizontal co-ordinate x,vertical co-ordinate y). This means that the first column of the

matrix x is plotted horizontally at the bottom, the middle column of x is plotted horizontally

in the middle, and the final column of x is plotted horizontally at the top.

These plotting details are not important for what will come. However, we give them here for

information.

10

2.3 A Little More Mathematical Notation

We have seen that in practice we observe a degraded version of a ‘true’ image such as:

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

We refer to this observed or degraded image using the vector y = (y1, . . . , yn), where n is the

number of pixels. Here yi ∈ {0,1}, i = 1, . . . , n.

11

2.4 Some Literature and Other Irrelevant Information

In scientific writing, we should always give credit where credit is due. Many people have

contributed to statistical image restoration. The following are but two references among a

great many:

Besag, J. E. (1986) On the statistical analysis of dirty pictures (with discussion). Journal of

the Royal Statistical Society B, 48, 259–302.

Geman, D. and Geman, S. (1984) Stochastic relaxation, Gibbs distributions and the Bayesian

restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6,

721–41.

Coursework 1 and Coursework 2 are based on the work of

Greig, D. M., Porteous, B. T. and Seheult, A. H. (1989) Maximum a posteriori estimation

for binary images. Journal of the Royal Statistical Society B, 51, 271–279.

Greig, Porteous and Seheult (1989) assume black is coded as 1 and white as 0. In Coursework

1 and Coursework 2, we use take 1 to represent white and 0 to represent black, as this is

more natural and usual.

The work of Greig, Porteous and Seheult (1989) was expanded on in these works:

Jubb, M. D. (1989) Image Reconstruction. PhD thesis. University of Bath.

Stander, J. (1992) Some Topics in Statistical Image Analysis. PhD thesis. University of Bath.

Stander, J. and Silverman, B. W. (1994) Temperature schedules for simulated annealing.

Statistics and Computing, 4, 21–32.

Franconi, L. and Jennison, C. (1997) Comparison of a genetric algorithm and simulated

annealing in an application to statistical image reconstruction. Statistics and Computing, 7,

193–207.

No knowledge of any of this literature is assumed in either Coursework 1 or

Coursework 2.

We end this section by remarking that images are not usually binary. Often each pixel can

take a grey level so that xi ∈ {0,1, . . . ,255} where 0 is black and 255 is white for example,

or a colour, which is more complicated to represent. We work with binary images to keep

life relatively simple, as this may be the first time that you have met image reconstruction

problems.

12

2.5 Statistical Models

Image restoration proceeds by specifying statistical models for possible true images x =(x1, . . . , xn) and for degraded images y = (y1, . . . , yn).

Please note that this is not a statistics module, so not all the details that we present here are

needed for you to be able to do Coursework 1 or Coursework 2. We include these details

to explain how the image restoration problem can be expressed as an optimization

problem. It is with that optimization problem that we shall work.

2.6 A Statistical Model for y given x

We will assume the following statistical model for yi given xi:

Pr(yi = 0 ∣xi = 0) = 1 − p (no error made, no swap)

Pr(yi = 1 ∣xi = 0) = p (error made, swap from 0 to 1)

Pr(yi = 0 ∣xi = 1) = p (error made, swap from 1 to 0)

Pr(yi = 1 ∣xi = 1) = 1 − p (no error made, no swap)

Essentially, we are assuming that the colour at pixel i is swapped from black to white, or

from white to black, with probability p ∈ [0,1]; otherwise, it is not changed. For simplicity,

we assume that p is known.

We assume that the error mechanism at each pixel is independent of the error mechanisms at

every other pixel. This enables us to write down Pr(y ∣x):

Pr(y ∣x) = n∏

i=1 Pr(yi ∣xi).

Let’s recall that the degraded image y is observed, while the ‘true’ image x is unknown and

needs to be estimated. Thought of as a function of the unknown image vector x, Pr(y ∣x)

is the likelihood function that you may have met in MATH2602 Statistical Inference and

Regression or MATH3613 Data Modelling. We shall write L(y ∣x) = Pr(y ∣x).

13

Now let’s write f(yi ∣xi) = Pr(yi ∣xi), so that f(1 ∣0) = Pr(yi = 1 ∣xi = 0) for example. We

have

L(y ∣x) = Pr(y ∣x)

= n∏

i=1 Pr(yi ∣xi)

= n∏

i=1 f(yi ∣1)xif(yi ∣0)1−xi ,

since f(yi ∣1)xif(yi ∣0)1−xi = f(yi ∣1) if xi = 1 and f(yi ∣0) if xi = 0.

From this it follows that

L(y ∣x) = n∏

i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}

xi

. (1)

Later, we will use λi, which is defined as follows:

λi = log{f(yi ∣1)

f(yi ∣0)} . (2)

Please note that λi only depends on the observed yi and not on xi. Please also note that, if

yi = 0, then

λi = log{f(yi ∣1)

f(yi ∣0)} = log{f(0 ∣1)f(0 ∣0)} = log ( p1 − p) ,

while, if yi = 1, then

λi = log{f(yi ∣1)

f(yi ∣0)} = log{f(1 ∣1)f(1 ∣0)} = log (1 − pp ) .

14

2.7 A Statistical Model for x

To define a statistical model for x, we need to explain what it means for two pixels to be

neighbours. This can be done by the following picture

Columns

R

ow

s

i = 1 or ( 1 , 1 )

i = 2 or ( 2 , 1 )

i = 3 or ( 3 , 1 )

i = 4 or ( 4 , 1 )

i = 5 or ( 1 , 2 )

i = 6 or ( 2 , 2 )

i = 7 or ( 3 , 2 )

i = 8 or ( 4 , 2 )

i = 9 or ( 1 , 3 )

i = 10 or ( 2 , 3 )

i = 11 or ( 3 , 3 )

i = 12 or ( 4 , 3 )

1 2 3

1

2

3

4

The red and blue lines indicate neighbouring pixles; we use blue to indicate horizontal

neighbours and red to indicate vertical neighbours, although this distinction is not important.

As an example we see that the neighbours of pixel i = 6 are 5, 10, 7 and 2. It turns out

that it simpler to assume that the image is wrapped on a torus, so that, for example the

neighbours of pixel i = 1 are 4 (wrapped top to bottom), 5, 2 and 9 (wrapped left to right).

This assumption has very little practical effect for large images, but it does mean that every

pixel has four neighbours.

15

These neighbourhood relationships lead to a symmetric adjacency matrix A, which here

will be of dimension n × n = 12 × 12:

## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]

## [1,] 0 1 0 1 1 0 0 0 1 0 0 0

## [2,] 1 0 1 0 0 1 0 0 0 1 0 0

## [3,] 0 1 0 1 0 0 1 0 0 0 1 0

## [4,] 1 0 1 0 0 0 0 1 0 0 0 1

## [5,] 1 0 0 0 0 1 0 1 1 0 0 0

## [6,] 0 1 0 0 1 0 1 0 0 1 0 0

## [7,] 0 0 1 0 0 1 0 1 0 0 1 0

## [8,] 0 0 0 1 1 0 1 0 0 0 0 1

## [9,] 1 0 0 0 1 0 0 0 0 1 0 1

## [10,] 0 1 0 0 0 1 0 0 1 0 1 0

## [11,] 0 0 1 0 0 0 1 0 0 1 0 1

## [12,] 0 0 0 1 0 0 0 1 1 0 1 0

Do you understand why?

16

We can now define a probability mass function on images x as

Pr(x) = exp (−β × number of discrepant neighbouring pairs in x)∑images x′ exp (−β × number of discrepant neighbouring pairs in x′) ,

in which β ≥ 0 is a constant that we will discuss later. The sum over all images x′ comprises

2n terms, where n is the number of pixels, and can be very hard to compute when n is large.

As this sum does not depend on the argument x of Pr(x), we can discard it, if we replace

the equality with a ‘proportional to’ sign. So we write

Pr(x)∝ exp (−β × number of discrepant neighbouring pairs in x) .

For the image x that we have been working with:

## [,1] [,2] [,3]

## [1,] 0 0 0

## [2,] 0 0 1

## [3,] 1 1 1

## [4,] 0 0 1

the number of discrepant neighbouring pairs is 10:

Columns

R

ow

s

0

0

1

0

0

0

1

0

0

1

1

1

1 2 3

1

2

3

4

Hence, the probability of this image is Pr(x) = C exp(−β × 10), where C is a constant.

17

Similarly, the number of discrepant neighbouring pixels in the all black and the all white

images is, unsurprisingly, 0:

Columns

R

ow

s

0

0

0

0

0

0

0

0

0

0

0

0

1 2 3

1

2

3

4

Columns

R

ow

s

1

1

1

1

1

1

1

1

1

1

1

1

1 2 3

1

2

3

4

This means that the all black and the all white images are the most probable with probabilities

C exp(β × 0) = C each.

18

The chess-board image is the least probable, as it has the highest number of discrepant

neighbouring pixels:

Columns

R

ow

s

1

0

1

0

0

1

0

1

1

0

1

0

1 2 3

1

2

3

4

There are 20 discrepant pairs, meaning that the probability of the chess board image is

C exp(−β × 20).

You may be asking: what does β ≥ 0 do? The larger the value of β, the higher the relative

probabity given to the all black and to the all white image. In other words, higher β values

encourage ‘smoother’ images. Because of this, sometimes β is referred to as a smoothing

parameter.

Finally, we remark that our probability model for images x,

Pr(x)∝ exp (−β × number of discrepant neighbouring pairs in x) ,

is essentially the same as the Ising model used a lot in physics and statistical mechanics

to represent substances such as magnets. To quote Federico Camia: ‘The Ising model,

introduced by Lenz in 1920 to describe ferromagnetism, is one of the most studied models

of statistical mechanics. Its two dimensional version has played a special role in rigorous

statistical mechanics since Peierls’ proof of a phase transition in 1936 and Onsager’s derivation

of the free energy in 1944.’

19

It can be shown that the

number of discrepant neighbouring pairs in x = 12 n∑i=1 n∑j=1aij (xi − xj)2 ,

in which x = (x1, . . . , xn) is an image and aij is the (i, j)-th element of the adjacence matrix

A such that

aij = {1 if pixels i and j are neighbours0 if pixels i and j are not neighbours.

A consequence of this result is that

Pr(x) ∝ exp (−β × number of discrepant neighbouring pairs in x)

= exp(−β × 12 n∑i=1 n∑j=1aij (xi − xj)2)

= exp(−β2 n∑i=1 n∑j=1aij (xi − xj)2) . (3)

Another irrelevant fact is that more general probability mass functions can be defined by

replacing (xi − xj)2 with φ (xi − xj) where the function φ satisfies certain properties.

20

2.8 A Statistical Model for x given y

It makes sense to consider Pr(x ∣ y), as we are given the observed image y. Pr(x ∣ y) can be

found using Bayes’ theorem as follows:

Pr(x ∣ y) ∝ L(y ∣x) ×Pr(x)

∝ n∏

i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}

xi ×Pr(x) from (1)

∝ n∏

i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}

xi × exp(−β2 n∑i=1 n∑j=1aij (xi − xj)2) from (3) (4)

In the language of Bayesian statistical inference, Pr(x) is the prior probability mass

function of all images x before observing the degraded image y. Pr(x) expresses our belief

about x before seeing the data in terms of probabilities. Pr(x ∣ y) is our belief about x after

seeing the data y and is referred to as our posterior probability mass function.

Please note that, as mentioned, it does not matter if you do not follow all these details,

although they are not hard. The important part is that image restoration leads to an

optimization problem, as we shall see.

21

2.9 The Restoration of a Binary Image as an Optimization Problem

The restored image is defined as the most probable image under the posterior Pr(x ∣ y).

That is, it is the image x that maximizes Pr(x ∣ y). For this reason it is referred to as the

maximum a posteriori or MAP image.

From (4) we see that the restored image maximizes

Pr(x ∣ y)∝ n∏

i=1 f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}

xi × exp(−β2 n∑i=1 n∑j=1aij (xi − xj)2)

over all images x.

We know that it’s very often easier to work with the logarithm of a function, rather than

the function itself, when maximizing. This is because the logarithm of a function and the

function itself are maximized at the same point, and also because taking logarithms turns

products into sums, which are easier to deal with. So, the restored image maximizes

log {Pr(x ∣ y)} = additive constant D + n∑

i=1 log [f(yi ∣0) × {f(yi ∣1)f(yi ∣0)}

xi] − β2 n∑i=1 n∑j=1aij (xi − xj)2

= D + n∑

i=1 log {f(yi ∣0)} + n∑i=1 log [{f(yi ∣1)f(yi ∣0)}

xi] − β2 n∑i=1 n∑j=1aij (xi − xj)2

= D + n∑

i=1 log {f(yi ∣0)} + n∑i=1 xi log{f(yi ∣1)f(yi ∣0)} − β2 n∑i=1 n∑j=1aij (xi − xj)2

= D + n∑

i=1 log {f(yi ∣0)} + n∑i=1 xiλi − β2 n∑i=1 n∑j=1aij (xi − xj)2 (5)

by the definition of λi in (2).

Please note that ∑ni=1 log {f(yi ∣0)} does not depend on x, and so can be considered a constant.

It therefore follows that we need to maximize

n∑

i=1 xiλi − β2 n∑i=1 n∑j=1aij (xi − xj)2

over x = (x1, . . . , xn), or, equivalently, to minimize

− n∑

i=1 xiλi + β2 n∑i=1 n∑j=1aij (xi − xj)2

over all images x.

22

So, we have finally arrived at our image restoration optimization problem. Given a degraded

image y = (y1, . . . , yn), we can find our restored image by solving the problem:

minimize: − n∑

i=1 xiλi + β2 n∑i=1 n∑j=1aij (xi − xj)2 over images x = (x1, . . . , xn). (6)

It can be shown that the minimiation problem given by (6) is equivalent to the problem

minimize:

n∑

i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2 (7)

over images x, in which max(a, b) means the maximum of a and b, so that, for example,

max(0,1) = 1 and max(0,−1) = 0.

This quantity is never negative.

23

2.10 Solving the Minimization Problem using the Ford-Fulkerson

Algorithm

We will now see how Greig, Porteous and Seheult (1989) solved the minimization problem

(7), namely,

minimize:

n∑

i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2

over images x, using the Ford-Fulkerson algorithm. To do this, Greig, Porteous and Seheult

(1989) defined a network with capacities as follows:

2.10.1 Definition of the Network with Capacities

• There are n + 2 vertices:

– A source s;

– A sink t;

– A vertex for each of the n pixels, indexed by i or j below.

• If λi > 0, there is an arc (directed edge) (s, i) from the source s to i with capactiy

csi = λi; otherwise

• If λi < 0, there is an arc (i, t) from i to the sink t with capactiy cit = −λi;

• There is an undirected edge {i, j} between two internal vertices i and j (that corre-

spond to pixels i and j) with capacities cij = β aij. This can be thought of as two arcs,

one between i and j, and the other between j and i, both of which have capacities β aij .

• There are no other arcs.

24

2.10.2 Definition of a Cut

Let us say that x = (x1, . . . , xn) is an image that is being proposed as a reconstruction of the

observed y = (y1, . . . , yn).

• Partitioning the Network

We can use x to partition the network into two sets:

W (x) = {s} ∪ {i ∶ xi = 1}, source s and all white pixels

B(x) = {i ∶ xi = 0} ∪ {t}, all black pixels and the sink t.

This partition of the vertices in the network into the sets W (x) and B(x) defines a cut.

• The Cut Capacity

The associated cut capacity is

C(x) ∶= C(∂+(W (x))) = ∑

kl∈∂+(W (x)) ckl. (8)

Notice that the arcs of ∂+(W (x)) are of three types: arcs starting at s; arcs ending at t; arcs

between intermediate vertices.

The cut capacity (8) can also be written as

C(x) = ∑

k∈W (x) ∑l∈B(x) ckl.

It turns out that

C(x) = n∑

i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2 ,

which is exactly the function that we wanted to minimize in (7) to restore the image from y.,

So, the image restoration problem becomes a max-flow/min-cut problem: we find the cut

that minimizes the cut capacity C(x) and this determines x. We can do this using the

Ford-Fulkerson algorithm (or modern improvements of it). ,

We will expand on and explore the above in Coursework 1 and Coursework 2.

25

3 Coursework 1

Marks will be assigned for Coursework 1 according to this marking grid.

MATH3609 Optimization, Networks and Graphs:

Coursework 1 Marking Grid

Assessment Area Maximum Awarded Feedback

Mark Mark

Task 1: Correctness and va-

lidity of argument. Read-

ability of explanation.

15

Task 2: Correctness and va-

lidity of argument. Read-

ability of explanation.

20

Task 3: Correctness and va-

lidity of argument. Read-

ability of explanation.

10

Task4: Correctness and va-

lidity of argument. Read-

ability of explanation.

25

Task 5: Correctness and va-

lidity of argument. Read-

ability of explanation.

15

Task 6: Correctness and va-

lidity of argument. Read-

ability of explanation.

15

Total1 100

1This mark is provisional. Like all marks it is considered by an External Examiner and an Assessment Panel.

26

3.1 Coursework 1 Tasks

Coursework 1 comprises a number of tasks. You need to produce a document containing

your solutions to these tasks. These may be hand-written and then scanned or typed up

using, for example, LATEX. It is your responsibility to make sure that they are readable.

The page limit is fifteen pages. Please note that this is the limit and should not be

considered as the target. Well written, clear and concise mathematics will receive high credit.

Please arrange your report under subheadings referring to the tasks, that is, Task 1,

Task 2, . . .. This will help you to respond systematically to all parts of Coursework 1.

The marks for each task are shown in the grid on page 26.

Task 1: Let G = (V,E) be a graph with vertex set V = {1,2, . . . , n} and let A be the

adjacency matrix of G with (i, j)-th entry aij. Suppose x∶V → {0,1} is a function which

assigns to each vertex i a number xi, which is either 0 or 1. A pair of neighbouring vertices i

and j are called discrepant if xi ≠ xj.

Explain why the number of discrepant neighbouring (unordered) pairs is given by the formula:

1

2

n∑

i=1

n∑

j=1aij (xi − xj)2 .

Task 2: Explain why the minimization problem given by (6) on page 23 is equivalent to the

problem

minimize:

n∑

i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2

over images x, in which max(a, b) means the maximum of a and b, so that, for example,

max(0,1) = 1 and max(0,−1) = 0.

Task 3: Prove that

n∑

i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2

is never negative.

Task 4: Prove that

C(∂+(W (x))) = n∑

i=1(1 − xi)max(0, λi) + n∑i=1 ximax(0,−λi) + β2 n∑i=1 n∑j=1aij (xi − xj)2 .

Hint: The arcs of ∂+(W (x)) can be split into three types, namely, arcs that start at s, arcs

that end at t and arcs between intermediate vertices. Consider the sums of the capacities

arcs of each type separately.

27

Task 5: Show that every cut K of the network can be written as

K = ∂+(W (x))

for some image x, and deduce using this and the previous tasks that the Ford-Fulkerson

algorithm can be used to find a restored image.

Task 6: Is there always a unique solution x to the minimization problem (i.e. a unique

restored image)? Illustrate your answer with an example.

28

3.2 Coursework 1: What You Need to Submit

Your group needs to submit the following file electronically using the DLE:

• A Portable Document Format file MATH3609_Coursework_1_2020_2021.pdf containing

your report for Coursework 1.

All members of the group should take responsibility for seeing that your work

is submitted before the deadline.

It is also the responsibility of all members to check that the correct file has

been uploaded and that it is readable.

If anything is unclear, you should ask without delay.

29

4 Coursework 2

Marks will be assigned for Coursework 2 according to this marking grid.

MATH3609 Optimization, Networks and Graphs:

Coursework 2 Marking Grid

Assessment Area Maximum Awarded Feedback

Mark Mark

Image Restoration

Tasks: Defining the set-

up. Correctness, clarity

and presentation of code.

25

Image Restoration

Tasks: Optimizing,

presenting the results,

choosing β. Depth of

understanding and quality

of description of insights

gained.

50

Report: display results

clearly, concisely, system-

atically and professionally,

using correct spelling and

grammar. This is important

for employability.

25

Total2 100

2This mark is provisional. Like all marks it is considered by an External Examiner and an Assessment Panel.

30

4.1 Coursework 2 Tasks

Coursework 2 comprises a number of tasks. You need to produce a report of your work

following the instructions below.

Your report should contain well presented and annotated R code.

The page limit is twenty pages, including R code and figures. Please see page 39

for more details. Please do not submit an additional appendix as it will not be considered.

Reports that contain irrelevant or uninteresting discussion or R code will be penalized. Please

note that twenty pages is the limit and should not be considered as the target. Well written,

concise reports will receive high credit.

You may also omit lengthy, technical R code, if you feel that it does not add to

the argument that you are making. Essential code should, however, be shown.

Please arrange your report under subheadings referring to the tasks, that is, Task 1,

Task 2, . . .. This will help you to respond systematically to all parts of Coursework 2.

31

4.1.1 Implementing the Ford-Fulkerson Algorithm to Restore the

Image in R: Getting Started

Following on from the information in the section ‘Image Restoration: Preparation Work’ we

will now produce an image restoration from observed data y using the Ford-Fulkerson

algorithm implemented in R.

First, we need to read y into R. The data that we will use are supplied in the file

R_degraded_image_0.2.csv. We can read in and display these data in the same was

as we did at the beginning of this document. Here is the code again to make sure that you

get started properly:

y <- as.matrix(read.table("R_degraded_image_0.2.csv",

sep = ","))

#

image(y,

col = c("black", "white"))

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

#

n_row <- nrow(y)

n_col <- ncol(y)

#

c(n_row, n_col)

## [1] 23 28

#

y_vec <- as.vector(y)

#

length(y_vec)

## [1] 644

32

We can label the axes in a possibly nicer way:

image(1:n_row,

1:n_col,

y,

xlab = "Row number",

ylab = "Column number",

main = "Data: degraded image y",

col = c("black", "white"))

5 10 15 20

5

10

15

20

25

Data: degraded image y

Row number

Co

lu

m

n

nu

m

be

r

• It is not hard to imagine that y is a degraded version of the image in R_logo seen above.

• It is known that p = 0.2. (In practice, p and any other unknown parameter or parameters

can also be estimated form the data by assigning prior distributions to them and sampling

from the resulting posterior using Markov chain Monte Carlo methods as discussed in

MATH3613 Data Modelling. This is beyond the scope of Coursework 2.)

• We will assume for now that β = 1.

33

4.1.2 Defining the Vector of λ Values

Task 1:

Using the function ifelse or otherwise, define the vector of λ values (λ1, . . . , λn) in R and

show a few λ values. Identify and record in R which λi values are positive and which λi values

are negative using the which function, storing the results in vectors which_lambda_positive

and which_lambda_negative respectively.

4.1.3 Defining the Network with Capacities

We need to define a data frame comprising three variables or columns called from, to and

capacity. This data frame defines the capacities in the network. If, for example, we wish to

specify ckl, then the variable from would take the value k, the variable to would take the

value l and the variable capacity would take the value ckl.

One way to define the required data frame of capacities is to construct three data

frames df_source, df_sink and df_internal, and to combine df_source, df_sink and

df_internal by binding their rows into one data frame df using rbind.

The data frame df_source records the capacities of arcs coming from the source s.

Task 2: Define df_source in R. The square brackets [] may help you here. In order to

check this, the first few rows of the result are given

## from to capacity

## 1 source 2 1.386294

## 2 source 3 1.386294

## 3 source 4 1.386294

## 4 source 6 1.386294

## 5 source 7 1.386294

## 6 source 8 1.386294

The data frame df_sink records the capacities of arcs going to the sink t.

Task 3: Define df_sink in R and show a part of it.

34

Before we define df_internal, we need to obtain the adjacency matrix A. This can be

done by using the function get_neighbourhood_matrix with which you are supplied in the

file of helper functions helper_functions.R; you do not need to understand the details of

these helper functions. Please note that get_neighbourhood_matrix is defined using other

functions, so that you’ll have to make sure that all the functions in the file of helper functions

helper_functions.R are loaded into R. We also need to extract the rows and columns of A

corresonding to neighbouring pixels:

A <- get_neighbourhood_matrix(y, n_row)

#

A[1:5, 1:5]

## [,1] [,2] [,3] [,4] [,5]

## [1,] 0 1 0 0 0

## [2,] 1 0 1 0 0

## [3,] 0 1 0 1 0

## [4,] 0 0 1 0 1

## [5,] 0 0 0 1 0

#

neighbouring_pixels <- which(A == 1,

arr.ind = TRUE)

#

head(neighbouring_pixels)

## row col

## [1,] 2 1

## [2,] 23 1

## [3,] 24 1

## [4,] 622 1

## [5,] 1 2

## [6,] 3 2

35

The data frame df_internal records cij where i and j correspond to pixels i and j.

Task 4: Define df_internal in R and show a part of it. Define the data frame df using

rbind. Use the function graph_from_data_frame from the igraph package to create an

igraph graph called g from df. Produce a plot of your capacitated network, as shown here;

please note, that it’s very big (but rather nice)! Please note that your plot may be displayed in

a different way from the one here. You may find https://igraph.org/r/doc/plot.common.html

helpful.

source

2

23

24

622

1

3

25

623

4

26

624

5

27

625

6

28

626

7

29

627

8

30

628

9

31

629

10

32

630

11

33

631

12

34

632

13

35633

14 36

634

15

37

635

16

38

636

17

39 637

18

40

638

19

41

639

20

42

640

21

43

64122

44

642

45

643

46

64447

48

49

50

51

52

53

54

55

56

57

58

5960

6162

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

8283

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124 125

126

127

128

129

130

131132

133

134

135

136

137

138

139

140

141

142

143

144

145

146 147

148

149

150

151

152153

154

155

156

157

158

159

160

161 162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195196

197198

199

200

201

202

203

204

205

206

207 208

209

210

211

212

213

214

215

216

217218219220

221

222

223

224

225

226

227

228

229

230 231

232 233

234

235

236

237

238

239

240

241242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264265

266

267

268

269

270

271

272

273

274 275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300 301

302

303

304

305

306

307

308

309

310

311

312313314

315

316

317

318

319

320

321 322

323

324

325

326

327

328

329330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384385

386

387

388

389

390

391

392

393

394

395

396

397

398

399400

401

402

40

404

405

406

407

408409

410

411

412

413

414

415

416

417

418

419

420

421

422423

424

425

426

427

428

429

430

431

432

433

434

435

436 437

438

439

440

441

442

443

444

445446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

47

472

473

474

475

476

477

478

79

480

481

482

483

484

485486487

488

489

490

491 492

493

494

495

496

497498

499

500

501

502

503

504

505

506

507

508

509

510

511

512513

514

515

516

517

518

519

520521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

540

541

542543

544

545

546

547

548

549550

551

552553

554

555

556

557

558

559

560

561

562

563

564

565566567

568

569

570

571

572573574

575

576

577

578

579

580

581

582

583

584

585

586 587

588589590

591

592

593

594

595

596

597

598

599

600601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

616

617

618

619

620

621

sink

36

4.1.4 Producing the Image Restoration by Maximizing Flow or Minimizing Cut Capacity

The function max_flow provided by the igraph package finds the maximum flow in a graph

or equivalently the minimum cut capacity. max_flow effectively implements an enhanced

modern version of the Ford-Fulkerson algorithm. Amongst other output, max_flow

returns a partition of the vertices into two sets that it calls partition1 and partition2.

partition1 is the set W (x) and partition2 is the set B(x) defined above, corresponding to

the ‘maximum flow’/‘minimum cut’ x. partition1 and partition2 are quite complicated

R objects. However, the part of interest to us can be extracted using the names function;

a character vector is given by names. Once names has been applied and any unnecessary

element has been removed, the resulting character vector can be turned into pixels numbers

using the as.numeric function.

Task 5: Use max_flow to maximize the flow in g. Report the value of the maximum flow.

Working with partition1 find the numerical values of the white pixels in the set W (x),

corresponding to the max-flow/min-cut, and report the first six of these. Similarly, find the

numerical values of the black pixels in the set B(x) and report the first six.

Task 6: Use these pixel values to define the reconstructed image. An easy way (but not the

only way) to do this is to use rep to define a vector, x_map_vec say, of the same length as

y_vec comprising NAs, and then to assign the value 0 to pixels that have been identified as

black and the value 1 to pixels that have been identified as white. Then turn the vector into

a matrix, remembering that pixels are labelled column by column; please see the help file for

matrix, if necessary. Display your image. Here is the one that we obtained.

0.0 0.2 0.4 0.6 0.8 1.0

0.

0

0.

2

0.

4

0.

6

0.

8

1.

0

37

4.1.5 Choosing β

This reconstruction x obtained from the observed image y was obtained using β = 1. In fact,

the choice of β is in the hands of the user.

Task 7:

Discuss and illustrate how you could choose β.

There is no right answer here. What we are looking for is sensible suggestions for

the choice of β, discussed clearly and critically. Your discussion should be illus-

trated with useful figures, as appropriate. Concise answers comprising sensible

and well illustrated suggestions will receive high credit.

38

4.2 Coursework 2: Report Production

The ability to report technical results clearly and concisely is an important skill.

You should write a report that, as a minimum:

• discusses in detail the above tasks and

• contains well presented R code developed in responding to the tasks.

You are advised, but not obliged, to use R Markdown. We will certainly help you to get

started with R Markdown, if you ask. You can use LATEX from within R Markdown, but you

are under no obligation to use LATEX.

You should use clearly readable fonts and maintain sensible margins and line spacing through-

out. Your work should not look squashed. The default settings of R Markdown are completely

acceptable.

4.3 Coursework 2: What You Need to Submit

Your group needs to submit the following file electronically using the DLE:

• A Portable Document Format file MATH3609_Coursework_2_2020_2021.pdf containing

your report for Coursework 2.

All members of the group should take responsibility for seeing that your work

is submitted before the deadline.

It is also the responsibility of all members to check that the correct file has

been uploaded.

There have been occasions when R code has been submitted instead of the

required report. This error is easy to make, but will have severe consequences.

If anything is unclear, you should ask without delay.

39

学霸联盟