#

Chamin Hewa Koneputugodage,1 Rhys Healy,1 Sean Lamont,1 Ian Mallett,1
Matt Brown,1 Matt Walters,1 Ushini Attanayake,1 Libo Zhang,1,2
Roger T. Dean,3 Alexander Hunter,1 Charles Gretton,1 Christian Walder,1,4
1Australian National University
2Nanjing University
3The MARCS Institute for Brain, Behaviour and Development, Western Sydney
University
4Data61, CSIRO
<EMAIL_ADDRESS>westernsydney.edu.au, smail.nju.edu.cn}

#

Chamin Hewa Koneputugodage,1 Rhys Healy,1 Sean Lamont,1 Ian Mallett,1
Matt Brown,1 Matt Walters,1 Ushini Attanayake,1 Libo Zhang,1,2
Roger T. Dean,3 Alexander Hunter,1 Charles Gretton,1 Christian Walder,1,4
1Australian National University
2Nanjing University
3The MARCS Institute for Brain, Behaviour and Development, Western Sydney
University
4Data61, CSIRO
<EMAIL_ADDRESS>westernsydney.edu.au, smail.nju.edu.cn}

###### Abstract

We address the problem of combining sequence models of symbolic music with
user defined constraints. For typical models this is non-trivial as only the
conditional distribution of each symbol given the earlier symbols is
available, while the constraints correspond to arbitrary times. Previously
this has been addressed by assuming a discrete time model of fixed rhythm. We
generalise to continuous time and arbitrary rhythm by introducing a simple,
novel, and efficient particle filter scheme, applicable to general continuous
time point processes. Extensive experimental evaluations demonstrate that in
comparison with a more traditional beam search baseline, the particle filter
exhibits superior statistical properties and yields more agreeable results in
an extensive human listening test experiment.

##  1 Introduction

The recent popularity of deep learning methods has given rise to progress in
computational models and algorithms for music generation (?; ?; ?). While
fully automated systems are interesting in their own right, many researchers
strive to combine machine learning models with human input (?; ?; ?; ?).

Of the various approaches to computer assisted music composition, a number of
broad categories are relevant to understanding the nature of the present
contribution. Perhaps the simplest scheme is the “continuator”, which allows
the user to provide a musical fragment that is then extended forward in time
by the algorithm (?; ?). Assuming that we have the usual causally factorised
model (that is, we have access to the conditional
$p(x_{t}|x_{1},x_{2},\dots,x_{t-1})$), these approaches sample the remainder
of a piece given its beginning. However, the continuator approach is
restrictive in that it may only condition on the start of the piece. This can
be overcome by two distinct families of alternative approaches.

The first family of approaches treats the completion of the musical piece as
an infilling problem (?; ?). Here the input is replaced with random masks, and
a model is learnt to fill in the masked segments. While interesting and
powerful, this approach conflates the modelling problem with the subsequent
composition step. As a result, it is unclear how best to train such a model.

The second general approach (to which this paper belongs) combines a typical
causal model of music with user defined constraints by introducing a
subsequent optimisation or conditioned sampling scheme. Pieces are
subsequently generated to satisfy these constraints while being representative
of the learned music model (?; ?). This has the advantage of permitting
standard models of (unconstrained) music, while placing the burden of
constraining the model on the composition algorithm. The challenging nature of
this subsequent step has meant that only relatively simple models, typically
with discrete time and fixed rhythmic structures, may be considered. For
example, limiting the Markov blanket of a given musical event to a small
temporal window surrounding it permits simple Gibbs sampling techniques (?),
but precludes the modelling of long range structure using advanced neural
sequence models.

In this paper we present the first generic scheme for combining user
constraints with a general continuous time point process model (?) of music
which enjoys the usual convenient causal factorisation structure, and which
has no limitation on the temporal range of the dependencies modeled. We
formulate the problem in Section 2, before introducing our main contribution,
a novel particle filter scheme (?), in Section 3. Section 4 explains how we
apply the generic particle filter scheme to the problem of constrained music
generation, Section 5 details extensive quantitative experiments, and Section
6 concludes the paper.

##  2 Problem Formulation

In the following sub-sections we first introduce some necessary notation and
definitions, and then formally define the conditional sampling problem
designed for generating computer assisted music in continuous time.

###  2.1 Notation and Basic Definitions

#### Notation

Bold such as $\bm{x}$ denotes a vector (or sequence) whose $i$-th element we
denote $x_{i}$. $\mathbb{R}^{*}=\bigcup_{i=0}^{\infty}\mathbb{R}^{i}$ is the
set of arbitrary length sequences, including the empty sequence. $\ell$
returns the length of its argument, so that
$\bm{x}\in\mathbb{R}^{d}\Rightarrow\ell(\bm{x})=d$. We also use $x_{\ell}$ to
denote the last element of $\bm{x}$, and $x_{\ell-1}$ the second last. We
combine this with slices such as $\bm{x}_{<i}$, $\bm{x}_{1:i-1}$ and
$\bm{x}_{<\ell}$, which have the usual meaning (first two being equivalent,
and the third being a sequence with all but the last element of $\bm{x}$).
Finally, in the algorithm pseudo-code we use e.g. $1:S$ to indicate iteration
over $1,2,\dots,S$.

#### Point Process Formulation

Suppose we have a general one dimensional point process (?) defined by

$\displaystyle x_{0}$ $\displaystyle=0$ (1) $\displaystyle x_{i}$
$\displaystyle\sim x_{i-1}+D(\bm{x}_{<i}),~{}~{}~{}i\in\mathbb{N}_{+},$ (2)

where $\bm{x}_{<i}=(x_{0},x_{1},\dots,x_{i-1})^{\top}\in\mathbb{R}^{i}$ and

$\displaystyle D:\mathbb{R}^{*}\rightarrow\Delta,$ (3)

where $\Delta$ is the set of random variables with samples in $(0,+\infty]$.
Hence, $D$ is a function which maps vectors (of the history of the sequence)
to inter-arrival distributions (for the next sample in the sequence). Let the
random variables $\bm{x}_{<i}$ generated in this way be distributed as

$\displaystyle\bm{x}_{<i}\sim X(D,i).$ (4)

This means that $X(D,i)$ is a random variable parametrised by $D$ and $i$,
with samples lying in $\mathbb{R}^{i}$. These samples are non-decreasing
sequences with first element equal to zero. The probability of a sequence
$\bm{x}_{1:i}$ can be written in terms of the conditional probability density
function (p.d.f.) as

$\displaystyle p(\bm{x}_{1:i})$
$\displaystyle\propto\prod_{k=1}^{i}f(x_{k}|\bm{x}_{<k}).$ (5)

#### Restriction Operator

Define the operator $\Pi_{\mathcal{I}}$ which restricts sequences to the set
$\mathcal{I}\subset\mathbb{R}$:

$\displaystyle\Pi_{\mathcal{I}}:\,$
$\displaystyle\mathbb{R}^{*}\rightarrow\mathbb{R}^{*}$ (6)
$\displaystyle\bm{x}\mapsto\Pi_{\mathcal{I}}(\bm{x})=\mathcal{I}\cap\bm{x}$
(7)

where $\mathcal{I}\cap\bm{x}$ means the subsequence containing all elements of
$\bm{x}$ lying in $\mathcal{I}$, in their original order. For example,
$\Pi_{[0,10]}((-1,0,1,\sqrt{117}))=(0,1)$.

#### Restricted Point Process

We define the restriction of our original process to the unit interval as
follows:

$\displaystyle\hat{X}(D)=\lim_{i\rightarrow\infty}\Pi_{[0,1]}\left(X(D,i)\right).$
(8)

Where we have abused the notation slightly, by applying $\Pi_{[0,1]}$ to a
random variable. In other words, by repeatedly sampling $x_{i}$ from (2) until
the first $i$ such that $x_{i}>1$, we obtain $\bm{x}_{<i}$ as a sample from
$\hat{X}(D)$.

###  2.2 Conditional Sampling Problem

Roughly speaking, we wish to sample from the restricted process $\hat{X}(D)$,
which is conditional on some fixed countable set of values
$Z=\\{z_{1}<z_{2}<\dots<z_{r}\\}\subset(0,1]$ being contained within the
sample. $Z$ represents the input provided by the user — e.g. $Z$ may represent
a melodic musical fragment which we wish to harmonise using our algorithm. We
introduce an additional degree of control by way of a binary indicator
sequence
$B=\\{b_{1},b_{2},\dots,b_{r}\\}\in\\{\operatorname*{True},\operatorname*{False}\\}^{r}$.
$b_{i}$ denotes whether the elements are allowed to be added between $z_{i}$
and $z_{i+1}$, which allows sections of the musical piece (contiguous ranges
of time) to be completely specified by the human. Also note that if the last
element $b_{r}=\operatorname*{False}$, then $z_{r}$ will be the final element
of the generated sequence. Define the indicator function

$\displaystyle\mathbbm{I}_{Z,B}(\bm{x})=$ (9) $\displaystyle\big{(}Z\subseteq
X\big{)}\text{ and }\big{(}\forall
i:b_{i}=\operatorname*{True},\,(z_{i},z_{i+1})\cap X=\emptyset\big{)},$ (10)

where $X\equiv\\{x_{i}\\}_{i=1}^{\ell(\bm{x})}$ is the set of elements in the
sequence $\bm{x}$. In words, $\mathbbm{I}_{Z,B}(\bm{x})$ is true iff $\bm{x}$
contains all of the elements of $Z$ while obeying the constraints indicated by
$B$. Then we may denote the random variable of interest by

$\displaystyle\hat{X}(Z,D),$ (11)

and define it implictly by way of its p.d.f., namely

$\displaystyle p_{\bm{x}\sim\hat{X}(Z,D)}(\bm{x}=\bm{x}^{\prime})\propto
p_{\bm{x}_{Z}\sim\hat{X}(D)}\left(\bm{x}_{Z}=\bm{x}^{\prime}\,|\,\mathbbm{I}_{Z}(\bm{x^{\prime}})\right),\quad$
(12)

where the right hand side involves conditioning on a zero measure set, and as
such should be interpreted as a regular conditional probability (?).

##  3 General Continuous Time Particle Filter

Algorithm 1 is a sequential Monte Carlo scheme for approximately sampling
$\hat{X}(Z,D)$. We now demonstrate the correctness of the algorithm.

#### Partitioned Formulation

Define the edge conditions $z_{0}=-\infty$ and $z_{R+1}=\infty$, and let

$\displaystyle r_{i}=\Pi_{(z_{i-1},z_{i}]}\,,$ (13)

and

$\displaystyle\mathcal{Z}:\,$
$\displaystyle\mathbb{R}^{*}\rightarrow(\mathbb{R}^{*})^{R+1}$ (14)
$\displaystyle\bm{x}\mapsto\mathcal{Z}(\bm{x})=\big{(}r_{1}(\bm{x}),r_{2}(\bm{x}),\dots,r_{R+1}(\bm{x})\big{)}.$
(15)

$\mathcal{Z}$ maps a sequence onto a sequence of sequences partitioned by the
$z_{i}$. Note that the inverse of $\mathcal{Z}$ is the concatenation of the
sequences. Now define

$\displaystyle Y(Z,D)=\mathcal{Z}(\hat{X}(Z,D)),$ (16)

where we have again abused the notation slightly by applying $\mathcal{Z}$ to
a distribution. By this we mean that applying $\mathcal{Z}$ to a sample from
$\hat{X}(Z,D)$ yields a sample from $Y(Z,D)$.

Algorithm 1 Our sequential Monte Carlo sampler for $\hat{X}_{Z}$. See
Algorithm 3 in the supplementary material for the SystematicResample function.

1:procedure ConditionalSample
input: $S$, $D$, $0<z_{1}<\dots<z_{R}\leq
1,b_{i}\in\\{\operatorname*{T},\operatorname*{F}\\};\,\forall i$
output: approximate samples $\bm{x}^{(1:S)}$ from $\hat{X}_{Z}$

2: for $i$ in $1,2,\dots,R+1$ do

3: for $s$ in $1,2,\dots,S$ do

4: if $b_{i-1}$ then notes allowed in $(z_{i-1},z_{i})$

5: while not $\big{(}x_{\ell}^{(s)}=z_{i}$ or $x_{\ell}^{(s)}\geq 1\big{)}$ do

6: sample $d\sim D(\bm{x}^{(s)})$ time increment

7: $\bm{x}^{(s)}\leftarrow(\bm{x}^{(s)},\min(z_{i},x_{\ell}^{(s)}+d))$

8: $d\leftarrow x_{\ell}^{(s)}-x_{\ell-1}^{(s)}$

9:
$w_{s}\leftarrow\frac{f_{D(\bm{x}_{<\ell})}(d^{(s)})}{\left(1-F_{D(\bm{x}_{<\ell})}(d^{(s)})\right)}$
from (32)

10: else

11: if $i\neq R+1$ thenif not the last element

12: $\bm{x}^{(s)}\leftarrow(\bm{x}^{(s)},z_{i})$ append

13: $d\leftarrow x_{\ell}^{(s)}-x_{\ell-1}^{(s)}$

14: $w_{s}\leftarrow f_{D(\bm{x}_{<\ell})}(d^{(s)})$

15: if $z_{i}\neq\infty$ thenif we are not finished

16: $k_{1:S}=\textsc{SystematicResample}(w_{1:S})$

17:
$\left(\bm{x}^{(s)}\right)_{s=1}^{S}\leftarrow\left(\bm{x}^{(k_{s})}\right)_{s=1}^{S}$
duplicate / delete

18: return $\bm{x}_{<\ell}^{(1:S)}$

Since samples from $Y(Z,D)$ are fixed length sequences (albeit with elements
that are arbitrary length sequences), we can appeal to sequential Monte Carlo
logic to approximately sample from it by approximating the density of $Y(Z,D)$
by a weighted empirical distribution of $S$ particles, namely

$\displaystyle\tilde{p}(\bm{y}_{1:R+1})=\sum_{s=1}^{S}\omega(\bm{y}_{1:R+1}^{(s)})\,\delta_{\bm{y}_{1:R+1}^{(s)}}(\bm{y}_{1:R+1}),$
(17)

where $\delta$ is the Dirac delta function. Note that $\bm{y}^{(s)}_{1:R+1}$
is the $s$-th particle, which consists of $R+1$ sequences
$\bm{y}^{(s)}_{i}\in\mathbb{R}^{*}$, for $i=1,2,\dots,R+1$. The normalised
weights $\omega$ are

$\displaystyle\omega(\bm{y}_{1:R+1}^{(s)})=\frac{w(\bm{y}_{1:R+1}^{(s)})}{\sum_{s^{\prime}=1}^{S}w(\bm{y}_{1:R+1}^{(s^{\prime})})}.$
(18)

#### Proposal

We take a factorised proposal in the partitioned sequence representation,

$\displaystyle
q(\bm{y}_{1:R+1})=\prod_{i=1}^{R+1}q(\bm{y}_{i}|\bm{y}_{<i},Z).$ (19)

The concrete proposal we employ, which is designed to yield a tractable
particle filter scheme, is most directly understood by its sampling procedure
which we provide as Algorithm 2 of the supplementary material (see Section
A.1). The basic idea of the proposal is to sample forward using the recursion
(2) until such time as the relevant partition boundary $z_{i}$ has been
crossed, at which point the element $z_{i}$ is included rather than the sample
which crossed that value. Mathematically, the proposal may be written

$\displaystyle q(\bm{y}_{i}^{(s)}|\bm{y}_{<i}^{(s)},Z)$
$\displaystyle=q(\bm{x}_{k+1:k+m}^{(s)}|\bm{x}_{<k+1}^{(s)},Z)$ (20)
$\displaystyle=\prod_{i=1}^{m}q(\bm{x}_{k+i}^{(s)}|\bm{x}_{<k+i}^{(s)},Z),$
(21)

where we assume there are $k$ elements $\bm{x}_{1:k}^{(s)}$ in
$\bm{y}_{1:i-1}^{(s)}$, i.e. $\sum_{j=1}^{i-1}|\bm{y}_{j}^{(s)}|=k$ and $m$
elements ($\bm{x}_{k+1:k+m}^{(s)}$) in $\bm{y}_{i}^{(s)}$. To match the
proposal sampler (Algorithm 2 of the appendix),

$\displaystyle
q(x_{k+1}^{(s)}|\bm{x}_{1:k}^{(s)},Z)=q(x_{k+1}^{(s)}|\bm{x}_{1:k}^{(s)},z_{i})^{b_{i}}\delta_{z_{i}}(x_{k+1}^{(s)})^{1-b_{i}},$
(22)

(where we have let the $b_{i}$ act as $\\{0,1\\}$ indicators), and

$\displaystyle q(x_{k+1}^{(s)}|\bm{x}_{1:k}^{(s)},z_{i})=$ (23) $\displaystyle
f(x_{k+1}^{(s)}|\bm{x}_{1:k}^{(s)})\mathbbm{I}(x_{k+1}^{(s)}<z_{i})+\delta_{z_{i}}(x_{k+1}^{(s)})(1-F(z_{i}|\bm{x}_{1:k}^{(s)})),$
(24)

where $F$ is the c.d.f. of $x_{k+m}^{(s)}$ and $\mathbbm{I}$ is the
$\\{0,1\\}$-valued indicator function. These expressions are most easily
understood by the sketches in Figure 1.

(a) $f(x_{k+m}^{(s)}|\bm{x}_{<k+m}^{(s)})$ (of the model)

(b) $q(x_{k+m}^{(s)}|\bm{x}_{<k+m}^{(s)})$ (of the proposal)

Figure 1: The model and proposal conditionals. For clarity we have omitted to
notate the conditioning in the annotations of the figures.

#### Particle Filter Resampling Weight

By sampling in the partitioned domain (of the $\bm{y}_{i}$), we have a fixed
length sequence (albeit of arbitrary length sequences), and may employ
standard importance sampling arguments (see e.g. ?; ? (?; ?)) to derive the
particle filtering weight of each sample,

$\displaystyle w(\bm{y}_{1:R+1}^{(s)})$
$\displaystyle=\frac{p_{\bm{y}_{1:R+1}^{(s)}\sim
Y(Z,D)}(\bm{y}_{1:R+1}^{(s)})}{q(\bm{y}_{1:R+1}^{(s)})}$ (25)
$\displaystyle=\frac{\prod_{i=1}^{R+1}p(\bm{y}_{i}^{(s)}|\bm{y}_{<i}^{(s)})}{\prod_{i=1}^{R+1}q(\bm{y}_{i}^{(s)}|\bm{y}_{<i}^{(s)})}.$
(26)

Note that we have omitted the subscript to $p$ in the second line. This can be
decomposed recursively as

$\displaystyle w(\bm{y}_{1:i}^{(s)})$
$\displaystyle=\underbrace{\frac{p(\bm{y}_{i}^{(s)}|\bm{y}_{<i}^{(s)})}{q(\bm{y}_{i}^{(s)}|\bm{y}_{<i}^{(s)})}}_{\doteq\,w(\bm{y}_{i}^{(s)}|\bm{y}_{<i}^{(s)})}\,\times\,\underbrace{\frac{\prod_{j=1}^{i-1}p(\bm{y}_{j}^{(s)}|\bm{y}_{<j}^{(s)})}{\prod_{j=1}^{i-1}q(\bm{y}_{j}^{(s)}|\bm{y}_{<j}^{(s)})}}_{\doteq\,w(\bm{y}_{<i}^{(s)})}.$
(27)

For $i\leq R$, assuming there are $k$ elements in $\bm{y}_{<i}$ and $m$
elements in $\bm{y}_{i}$ then the weight may be written

$\displaystyle w(\bm{y}_{i}^{(s)}|\bm{y}_{<i}^{(s)})$
$\displaystyle=\frac{p(\bm{y}_{i}^{(s)}|\bm{y}_{<i}^{(s)})}{q(\bm{y}_{i}^{(s)}|\bm{y}_{<i}^{(s)})}=\frac{\prod_{j\in\mathcal{J}}p(x_{j}^{(s)}|\bm{x}_{<j}^{(s)})}{\prod_{j\in\mathcal{J}}q(x_{j}^{(s)}|\bm{x}_{<j}^{(s)})}$
(28)
$\displaystyle=\frac{\prod_{j=k+1}^{k+m}p(x_{j}^{(s)}|\bm{x}_{<j}^{(s)})}{\prod_{j=k+1}^{k+m}q(x_{j}^{(s)}|\bm{x}_{<j}^{(s)})}$
(29)
$\displaystyle=\frac{p(x_{k+m}^{(s)}|\bm{x}_{<k+m}^{(s)})\prod_{j=k+1}^{k+m-1}p(x_{j}^{(s)}|\bm{x}_{<j})^{(s)}}{q(x_{k+m}^{(s)}|\bm{x}_{<k+m}^{(s)})\prod_{j=k+1}^{k+m-1}p(x_{j}^{(s)}|\bm{x}_{<j}^{(s)})}$
(30)
$\displaystyle=\frac{p(x_{k+m}^{(s)}|x_{<k+m}^{(s)})}{q(x_{k+m}^{(s)}|x_{<k+m}^{(s)})}$
(31)
$\displaystyle\propto\frac{f_{D(\bm{x}_{<\ell}^{(s)})}(d^{(s)})}{1-F_{D(\bm{x}_{<\ell}^{(s)})}(d^{(s)})},$
(32)

as per line 9 of Algorithm 1. Here, $f_{D}$ and $F_{D}$ are the p.d.f. and
c.d.f. of $x_{k+m}^{(s)}$ with the $d$ denoting the distance between the
variable $x_{k+m}^{(s)}$ and the former note $x_{k+m-1}^{(s)}$ respectively.
We have let $\bm{x}^{(s)}=\mathcal{Z}^{-1}\bm{y}_{1:i}^{(s)}$, recalling that
$\mathcal{Z}^{-1}$ is merely the concatenation operator, and let
$d^{(s)}=x_{\ell}^{(s)}-x_{\ell-1}^{(s)}$. $\mathcal{J}$ are the indices into
$\bm{x}^{(s)}$ corresponding to $\bm{y}_{i}^{(s)}$ (in increasing order). The
final expression (32) follows from the preceding line because all but the
final factor in the numerator and denominator cancel, due (loosely speaking)
to the proposal sampling all but the last element of $\bm{y}_{i}^{(s)}$ from
$p$ — a key simplification which arises from our appropriate choice of
proposal distribution. The denominator of (32) is due to the $\min$ operation
on line 12 of Algorithm 2 along with the definition of the c.d.f. Finally, we
denote (32) as a proportionality by neglecting an irrelevant (and generally
intractable) factor. This factor would merely renormalise the appropriate
joint, to obtain the conditional on the r.h.s. of (12).

For the final step $i=R+1$ note that the $\min$ operation on line 12 of
Algorithm 2 compares with $z_{i}=\infty$, thereby not altering the simple
recursion of (2). Hence, all factors in (27) cancel, so that all particle
weights equal one. In this case SystematicResample (Algorithm 3) returns the
sequence $1,2,\dots,S$, leading to an identity mapping as the re-sampling step
of line 15 of Algorithm 1. Hence, for this (final) step of Algorithm 1 we skip
the re-sampling step using the condition on line 15.

###  3.1 Alternative Derivation

We now sketch an alternative derivation of Algorithm 1 (for the simplified
case $b_{i}=\operatorname*{True},\forall i$) which offers an illuminating
connection with the point process literature (?; ?). Assuming the process is
simple (generates no duplicate points), then we may equivalently represent
samples of $\hat{X}(D)$ by way of an infinitesimally discretized indicator
vector, as follows. Recall that samples of $\hat{X}(D)$ are sequences with
$0=x_{0}<x_{1}<\dots<x_{\ell}\leq 1$. Let $N$ denote the resolution of the
discretization, and for $i=0,1,\dots,N-1$ define

$\displaystyle v_{i}=\begin{cases}1&\text{if
$x_{j}\in\big{(}i/N,(i+1)/N\big{]}$ for any $j$}\\\
0&\text{otherwise.}\end{cases}$ (33)

Consider the distribution on $\bm{v}$ implied by that on $\bm{x}$ and the
above (deterministic) mapping from $\bm{x}$ to $\bm{v}$. For any $D$ (of (9))
there is a corresponding $g_{D}:\left\\{0,1\right\\}^{*}\rightarrow[0,1]$ with

$\displaystyle v_{i}\sim\operatorname*{Bernoulli}(g_{D}(\bm{v}_{<i})).$ (34)

Our conditional sampling problem may now be viewed as one of sampling
sequences $\bm{v}\in\\{0,1\\}^{N}$ conditional on $v_{j}=1$ for all $j$
corresponding to the observed values $Z$. We apply sequential Monte Carlo to
this sampling problem. Let the proposal $q(v_{i}|\bm{v}_{{<i}})$
deterministically generate the value $1$ for the observed indices, and be
distributed according to (34) for the unobserved indices. By the arguments in
? (?), the associated particle weights are uniform (for the unobserved
indices), and equal to $g_{D}(\bm{v}_{<i})$ (for the observed indices). Now,
letting $k=\sum_{j=0}^{i-1}v_{j}$ be the index into $\bm{x}$ of the last point
accounted for by $\bm{v}_{<i}$, we have

$\displaystyle g_{D}(v_{i}|\bm{v}_{<i})$
$\displaystyle=\mathbb{P}[v_{i}=1|\bm{v}_{<i}]$ (35)
$\displaystyle=\mathbb{P}[x_{k+1}\in\big{(}i/N,(i+1)/N\big{]}\big{|}\bm{x}_{\leq
k}].$ (36)

In the limit $N\rightarrow\infty$ with $t=i/N$ this is by definition
$\lambda(t)\text{d}t$, where $\lambda$ is the conditional intensity function.
It is well known that (?)

$\displaystyle\lambda(t)=\frac{f^{*}(t)}{1-F^{*}(t)},$ (37)

where $f^{*}$ and $F^{*}$ are respectively the p.d.f. and c.d.f. of $x_{k+1}$
given $\bm{x}_{\leq k}$. We see that (37) is in agreement with (32) (and line
9 of Algorithm 1), and thereby provides an alternative intuition for Algorithm
1, namely that the resampling weight is proportional to the conditional
intensity function of the original point process.

##  4 Sequential Monte Carlo for Music

For simplicity, the previous section assumed a one-dimensional point process.
However, multiple notes can occur at any timestep in music (e.g. chords), and
as they may have arbitrary duration we need events to model their start and
end. As a result, we can view music as having two dimensions, time and
actions. The action space is an on or off action for each pitch. Although a
one-dimensional sampler was constructed in the last section, an appropriate
conversion can adapt it to sampling music encodings in this two-dimensional
fashion. The conversion is illustrated in Figure 2 and explained in the next
section.

###  4.1 Modelling an Event Sequence

We define our representation for music in two dimensions. Let a music event
(which is either a note on or off event) be denoted as
$e=(t(e),a(e)),t(a)\in[0,+\infty),a(e)\in[1,a_{max}]$. $t(a)$ represents the
time of the event, and $a(e)$ represents the type of event, which is an
action. We model music with fine temporal discretization for this conversion
(each corresponding to a $2400^{\text{th}}$ of a musical quarter note) and 128
different MIDI pitches. We set $a_{max}=256$ with $a(e)=m,m\in[1,128]$
representing turning on note $m$, and $a(e)=m,m\in[129,256]$ representing
turning off note $m-128$. To reduce to one dimension we map each event to a
unique integer value by setting $e=t(e)a_{max}+a(e)$. Given a one-dimensional
event $e$, we can then derive its time and action by
$t(e)=\lfloor(e/a_{max})\rfloor$, $a(e)=e\ \%\ a_{max}$, where $\%$ is the
modulus function.

(a) Two-dimensional events

(b) Unrolled one-dimensional events

Figure 2: A sketch of the reduction from two- to one-dimensional event
domains.

###  4.2 Modelling with an RNN

We model the distribution of music using a Recurrent Neural Network (RNN)
based model with a music representation inspired by the Performance RNN (?),
expanding the action space to include different parts (MIDI channels). Each
symbol modeled by the RNN is either an action (a musical note being turned on
or off in a specific part) or a shift in time between music events. Notes
turned on prior to the shift are kept on during the time shift, and there is
be silence otherwise.

#### Resolving Ambiguity

We improve slightly on ? (?) by ensuring that our representation for a given
musical sequence is unique — or roughly speaking that there is a bijection
between MIDI data and what the RNN sees. This is achieved by simply enforcing
that shifts not follow shifts, and at a given time, actions can only occur in
ascending order of enumeration.

Variable $t$ | $t(z)=t(x_{l-1})$ | $t(z)>t(x_{l-1})$
---|---|---
$x_{l}\sim\bm{x}_{<l}$ | sample to $a(z)$ | shift to time $t(z)$, then sample to action $a(z)$
$f_{\bm{x}_{<l}}(x_{l})$ | $f_{\bm{x}_{<l}}(a(z))$ | $f_{\bm{x}_{<l}}(t(z)-t_{x})\cdot f_{z}(a(z))$
$F_{\bm{x}_{<l}}(x_{l})$ | $\sum_{a(x)}^{a(z)}{f_{\bm{x}_{<l}}(a)}$ |  | $\sum_{a_{x}}^{a_{z}}{f_{\bm{x}_{<l}}(a)}+\sum_{1}^{t(z)-t(x)-1}{f_{\bm{x}_{<l}}(\Delta t)}+$
---
${f_{\bm{x}_{<l}}(t(z))\cdot\sum_{1}^{a(z)}{f_{z}(a)}}$
Table 1: Summary of the key quantities required by Algorithm 1. Here, $f_{z}$
(respectively $F_{z}$) denotes the conditional p.d.f. (c.d.f.) of the next
step, given that the particle has reached $t(z)$.

To this end, we define a constant mask matrix $M$ which has the following
structure:

$M=\begin{gathered}\end{gathered}{}$

and, denoting the output layer before the softmax layer as
$y_{l}=f_{\bm{x}_{<l}}(x_{l})$, we then replace $y_{l}$ with

$\displaystyle z_{l}=y_{l}+M(a(x_{l-1})),$ (38)

where $M(a(x_{l-1}))$ denotes the $a(x_{l-1})$-th column of $M$.

#### Unrolled Probabilities in terms of Event Probabilities

To apply Algorithm 1, we must write the required quantities as they pertain to
the unrolled music representation (that is, the bottom of Figure 2) in terms
of the probabilities available in the original, unrolled representation (the
top of Figure 2). Suppose the current event step in the sequence is $x$, then
the probability mass function (PMF) vector of the next potential step $f_{x}$
is obtained through a softmax layer that is divided into two parts:
$f_{x}=\big{[}f_{x}(a):a\in[1,a_{max}],f_{x}(\Delta t):\Delta
t\in[1,s_{max}]\big{]}$. $f_{x}(a)$ gives the probabilities for different
actions at the current time step $t(x)$, and $f_{x}(\Delta t)$ gives the
probabilities for different shift times. Each element of $f_{x}(\Delta t)$ is
the sum of the probabilities of all the actions at $t(x)+\Delta t$. The key
variables are summarized in Table 1, which provides the basis for computing
the weights of the particles in (32) as required in Algorithm 1.

##  5 Results

Algorithm 1 was used to generate music samples using the conversion for music
from Section 4. An RNN trained on the Symbolic Music Data Version 1.0 (SMD-V1)
(?) was applied for the sampling distribution $D$ in Algorithm 1. SMD-V1
consists of 5 different music datasets, and the validation sets of two
datasets (PMD and NOT) were used for testing. The pieces in these two datasets
are separated into two parts (MIDI channels). These are labelled as Part 0 and
Part 1. Using this combined validation set, the following 6 categories of
samples were generated:

GroundTruthSamples:

20s samples taken at random directly from samples in the validation sets with
both parts intact. Labelled as GT in the figures.

ConditionedPF-$\\{s\\}$:

40s samples were taken at random from the test set of the dataset. The RNN was
conditioned on he first 20s, and then Algorithm 1 was run with this
conditioned model as distribution $D$ and the notes in Part 0 of the second
20s as the constraints. $s$ particles were used in the algorithm. Labelled as
PF1, PF2, PF3 in order of ascending $s\in\\{30,100,300\\}$ in the figures.

ConditionedBeamSearch-${b}$-${f}$:

Used the same setup as ConditionedPF-$\\{s\\}$ but with a beam search based
algorithm in place of Algorithm 1. The beam search utilizes Algorithm 2 while
retaining multiple beams. In the loop (lines 5-9) of Algorithm 2 it samples
$b\times f$ trajectories at random (specifically $b$ samples from each the
previously kept $f$ trajectories), and then keeps the top $f$. The ranking is
based on maximising log probabilities of the current (incomplete)
trajectories. Labelled as BS1, BS2 in order of ascending $f$ in the figures.

#### Parameter choice

The parameters ($s,b$ and $f$) were chosen to give an interesting comparison
in the number of trajectories, amount of memory and time taken. In particular,
BS1 has as many trajectories as PF3, while BS2 has vastly more than the other
algorithms, giving it an unfair advantage. These performance metrics are
detailed in the Supplementary Material, Section A.3.

To measure the quality of the resulting samples, in the following sections we
investigate two aspects of the output: the probability of our samples under
the distribution of our underlying RNN, and the results of human listening
tests.

###  5.1 Probability of the Generated Sequences

From (5), we evaluate the log probability of a sample by

$\displaystyle\log p(\bm{x}_{1:i})$ $\displaystyle=\sum_{k=1}^{i}\log
f(x_{k}|\bm{x}_{<k}).$

where $f$ is our model distribution given by the RNN. A histogram of the log
probabilities of the samples from each category is shown in Figure 3.
BeamSearch’s distribution of log probabilities is much higher than
ConditionedPF’s, with larger proportions at $-100$ and smaller range (-400 to
-50 compared to -800 to -50 for ConditionedPF). Furthermore, as the number of
maintained trajectories increases, the distributions for both algorithms
become thinner and sharper, overall being more weighted towards higher log
probabilities. This is the difference between the beam search based approach,
which is actively maximising the whole subsequence between constraints, and
Algorithm 1 which is resampling based on the probability of constrained events
only. On the other hand, the GroundTruthSamples have a range from $-1000$ to
$-200$, and a mode at $-300$, and which is most similar to the ConditionedPF.
The distribution of GroundTruthSamples being very different to the BeamSearch
is the expected result of actively maximising log probabilities rather than
generating samples from the conditional distribution of interest as in the
particle filter method. This also does not allow for the diversity typical of
human music (?) and makes the samples rather bland and repetitive as we
demonstrate in the next section.

Figure 3: Histograms of log-probabilities.

###  5.2 Listening Tests

| BS1 | BS2 | PF1 | PF2 | PF3 | GT
---|---|---|---|---|---|---
BS1 | - | 0.51 | 0.51 | 0.54 | 0.49 | 0.39
BS2 | 0.49 | - | 0.52 | 0.44 | 0.42 | 0.44
PF1 | 0.49 | 0.48 | - | 0.54 | 0.52 | 0.46
PF2 | 0.46 | 0.56 | 0.46 | - | 0.59 | 0.40
PF3 | 0.51 | 0.58 | 0.48 | 0.41 | - | 0.48
GT | 0.61 | 0.56 | 0.54 | 0.60 | 0.52 | -
Table 2: Row beats column ratio for the “Better” test.  | BS1 | BS2 | PF1 | PF2 | PF3 | GT
---|---|---|---|---|---|---
BS1 | - | 0.42 | 0.33 | 0.41 | 0.33 | 0.28
BS2 | 0.58 | - | 0.35 | 0.37 | 0.37 | 0.30
PF1 | 0.67 | 0.65 | - | 0.46 | 0.40 | 0.40
PF2 | 0.59 | 0.63 | 0.54 | - | 0.45 | 0.34
PF3 | 0.67 | 0.63 | 0.60 | 0.55 | - | 0.43
GT | 0.72 | 0.70 | 0.60 | 0.66 | 0.57 | -
Table 3: Row beats column ratio for the “Interesting” test.

#### Setup

A listening test similar to ? (?) was used for the six categories. Sixteen 20s
samples from each category were taken for comparative listening tests between
each pair from different categories. After listening to the two samples
(labelled ‘A’ and ‘B’), the listeners were asked “Which sounds better
overall?” and “Which sounds more interesting?”. For both questions, the
listeners were given the choice ‘A’ or ‘B’. Both the order of the tests and
the order of the two samples within each test were randomised. Every
individual test was also given to 3 different listeners, making a total of
$3\times(^{6}_{2})\times 16^{2}=11,520$ listening tests.

#### Overall Wins

Overall win proportions can be seen in Figure 4. For the ‘Better’ listening
test, all categories did well (around 50%) while the GroundTruth category
achieved 56%. This shows that while people clearly percceived the ground truth
samples as superior, samples from each category were convincing enough to win
almost half the time. On the other hand, for the ‘Interesting’ test,
BeamSearch performed significantly lower than ConditionedPF (at almost 40%
wins compares to 50-60% wins), with ConditionedPF-300 performing very well
(57% wins) and GroundTruth clearly peforming the best (65% wins). This
suggests BeamSearch acts as expected of a greedy algorithm: it seeks to
optimise its likelihood (which intuitively would make it more likely to win
the ‘Better’ test), and therefore, it plays very monotonous but ‘safe’ pieces,
leading to uninteresting music samples. However, ConditionedPF is able to
sample more ‘risky’ and interesting pieces while still maintaining high
likelihoods, thus doing well on both tests, especially with higher particle
sizes.

#### Comparisons with the Ground Truth

Direct comparisons for each category versus the GroundTruth samples can be
seen in Figure 5. There is a larger difference between the BeamSearch and the
ConditionedPF on the ‘Better’ test, while the BeamSearch samples compare
poorly with the GroundTruth. The ConditionedPF with 300 particles performs the
best, nearly winning half of the tests (47.5%). Likewise, the BeamSearch rates
poorly on the ‘Interesting’ test, and again ConditionedPF-300 performs the
best (43%). Detailed comparisons between each pair of categories can be seen
in Tables 2 and 3.

Figure 4: Listening test wins against all other categories.  Figure 5:
Listening test wins vs. GroundTruth.

##  6 Conclusion

We have presented a novel sequential Monte Carlo algorithm for continuous time
point processes that are conditioned on a subset of events. By a suitable
choice of proposal distribution, we obtained Monte Carlo resampling weights
that are computationally tractable and simple to compute. We demonstrated how
this general approximate sampling algorithm may be adapted to symbolic music
generation, by letting the user specify a subset of musical events to be
included in the piece, along with controls on which parts of the piece to hold
fixed. Our results show that the proposed scheme produces more interesting
samples compared to a beam search baseline, without sacrificing musical
quality as perceived by human test subjects. The samples produced by our novel
approach sound more like human music, perform better under two listening test
metrics, and have a distribution which more closely resembles human
compositions when evaluated using our generative model for music.

## References

  * [2010] Baddeley, A.  2010\.  Multivariate and marked point processes.  Handbook of Spatial Statistics 371–402.
  * [2017] Briot, J.-P.; Hadjeres, G.; and Pachet, F.  2017\.  Deep learning techniques for music generation - a survey.  ArXiv abs/1709.01620.
  * [1999] Carpenter, J.; Clifford, P.; and Fearnhead, P.  1999\.  Improved particle filter for nonlinear problems.  IEE Proc., Radar Sonar Navig. 146(1):2.
  * [2009] Doucet, A., and Johansen, A. M.  2009\.  A tutorial on particle filtering and smoothing: Fifteen years later.  Handbook of nonlinear filtering 12(656-704):3.
  * [1985] Faden, A.  1985\.  The existence of regular conditional probabilities: Necessary and sufficient conditions.  In Annals of Probability.
  * [2013] Fernandez, J. D., and Vico, F. J.  2013\.  Ai methods in algorithmic composition: A comprehensive survey.  J. Artif. Intell. Res. (JAIR) 48:513–582.
  * [1993] Gordon, N.; Salmond, D.; and Smith, A.  1993\.  Novel approach to nonlinear/non-Gaussian Bayesian state estimation.  IEEE Proceedings F, Radar and Signal Processing 140(2):107–113.
  * [2017] Hadjeres, G.; Pachet, F.; and Nielsen, F.  2017\.  Deepbach: a steerable model for bach chorales generation.  In Precup, D., and Teh, Y. W., eds., Proceedings of the 34th International Conference on Machine Learning, ICML 2017, Sydney, NSW, Australia, 6-11 August 2017, volume 70 of Proceedings of Machine Learning Research, 1362–1371.  PMLR.
  * [2019] Hawthorne, C.; Stasyuk, A.; Roberts, A.; Simon, I.; Huang, C.-Z. A.; Dieleman, S.; Elsen, E.; Engel, J.; and Eck, D.  2019\.  Enabling factorized piano music modeling and generation with the MAESTRO dataset.  In International Conference on Learning Representations.
  * [2019] Huang, C.-Z. A.; Vaswani, A.; Uszkoreit, J.; Shazeer, N.; Simon, I.; Hawthorne, C.; Dai, A.; Hoffman, M.; Dinculescu, M.; and Eck, D.  2019\.  Music transformer: Generating music with long-term structure.
  * [2018] Ippolito, D.; Huang, A.; Hawthorne, C.; and Eck, D.  2018\.  Infilling piano performances.  NIPS Workshop on Machine Learning for Creativity and Design.
  * [2008] Nierhaus, G.  2008\.  Algorithmic Composition: Paradigms of Automated Music Generation.  Springer Publishing Company, Incorporated, 1st edition.
  * [2018] Oore, S.; Simon, I.; Dieleman, S.; Eck, D.; and Simonyan, K.  2018\.  This time with feeling: learning expressive musical performance.  Neural Computing and Applications.
  * [2010] Pachet, F.  2010\.  The continuator: Musical interaction with style.  Journal of New Music Research 32:333–341.
  * [2015] Papadopoulos, A.; Pachet, F.; Roy, P.; and Sakellariou, J.  2015\.  Exact sampling for regular and markov constraints with belief propagation.  In Proceedings of the 21st International Conference on Principles and Practice of Constraint Programming, CP’15, 341–350.  Switzerland: Springer.
  * [2016] Papadopoulos, A.; Roy, P.; and Pachet, F.  2016\.  Assisted lead sheet composition using flowcomposer.  In Rueher, M., ed., Principles and Practice of Constraint Programming, 769–785.  Cham: Springer International Publishing.
  * [2011] Rasmussen, J. G.  2011\.  Temporal point processes: the conditional intensity function.  Technical report, Aalborg University.
  * [2019] Roberts, A.; Engel, J.; Mann, Y.; Gillick, J.; Kayacik, C.; Nørly, S.; Dinculescu, M.; Radebaugh, C.; Hawthorne, C.; and Eck, D.  2019\.  Magenta studio: Augmenting creativity with deep learning in ableton live.  In Proceedings of the International Workshop on Musical Metacreation (MUME).
  * [2016] Walder, C., and Kim, D.  2016\.  Computer assisted composition with recurrent neural networks.  Asian Conference on Machine Learning.
  * [2016] Walder, C.  2016\.  Symbolic music data version 1.0.  CoRR abs/1606.02542.
  * [1987] Zicarelli, D.  1987\.  M and jam factory.  Computer Music Journal 11(4):13–29.

##  Appendix A Supplementary Material for “”

###  A.1 Proposal Sampler

We define the proposal $q(\bm{y}_{i}|\bm{y}_{<i})$ procedurally as Algorithm
2. The proposal samples appropriately from the interval $(z_{i-1},z_{i}]$ by
applying the recursion (2), but with a clipping operation on line 12 in order
to include $z_{i}$ in the proposed sequence with probability one. Once $z_{i}$
is included the algorithm terminates. Note that the edge conditions at lines 2
and 17 of Algorithm 2 vanish in the context of the complete Algorithm 1, which
is more elegant than the proposal distribution.

Algorithm 2 A sampler for the proposal $q(\bm{y}_{i}|\bm{y}_{<i})$

1:procedure Proposal
input: $i,\bm{y}_{<i}$, $D$, $0<z_{1}<\dots<z_{R}\leq 1$, $b_{i-1}.$
output: $\bm{y}_{i}\sim q(\bm{y}_{i}|\bm{y}_{<i})$

2: if $i=1$ then

3: $x\leftarrow 0;\quad\bm{y}_{i}^{(s)}\leftarrow(0)$

4: else

5: $x\leftarrow(\bm{y}_{i-1})_{\ell};\quad\bm{y}_{i}^{(s)}\leftarrow()$

6: $x\leftarrow(\bm{y}_{i-1})_{\ell}$ last value of previous subsequence

7: $\bm{y}_{i}^{(s)}\leftarrow()$ length zero sequence

8: if $b_{i-1}$ then notes allowed in $(z_{i-1},z_{i})$

9: while $\big{(}x<z_{i}\big{)}$ do

10: $\bm{x}^{(s)}\leftarrow\mathcal{Z}^{-1}(\bm{y}_{1:i}^{(s)})$
$\mathcal{Z}^{-1}$ concatenates

11: sample $d\sim D(\bm{x}^{(s)})$

12: $\bm{y}_{i}\leftarrow(\bm{y_{i}},\min(z_{i},x+d))$ append

13: $x\leftarrow(\bm{y}_{i})_{\ell}$

14: else

15: $\bm{y}_{i}\leftarrow(\bm{y_{i}},z_{i})$ append

16: return $\bm{y}_{i}$

17: if $x>1$ then

18: return $(\bm{y}_{i})_{<\ell}$

19: else

20: return $\bm{y}_{i}$

Algorithm 3 Systematic Resampling
(?).

1:procedure SystematicResample
input: $w_{1:S}$
output: $k_{1:S}$

2: $\omega_{1:S}\leftarrow w_{1:S}/\sum_{s=1}^{S}w_{s}$

3: sample $u\sim\textsc{Uniform}([0,1])$

4: $\bar{u}\leftarrow u/S$

5: $j\leftarrow 1$

6: $S_{\omega}\leftarrow\omega_{1}$

7: for $l$ in $1,2,\dots,S$ do

8: while $S_{\omega}<\bar{u}$ do

9: $j\leftarrow j+1$

10: $S_{\omega}\leftarrow S_{\omega}+\omega_{j}$

11: $k_{l}\leftarrow j$

12: $\bar{u}\leftarrow\bar{u}+1/S$

13: return $k_{1:S}$

###  A.2 Music Samples and Listening Tests

The samples used for the listening test can be found at

    http://bit.ly/GenMusicSamples

It is separated into 6 folders, one for each category, and within each there
are 16 music samples rendered into mp3 format.

Furthermore, in our supplementary material we provide four selected music
samples from each category, which correspond to the four most interesting
pieces as determined by the “interesting” listening tests.

The output music events from our algorithm were converted into MIDI (our
representation gives a unique representation to MIDI). These were then
synthesized using a custom SoundFont based on piano, and then converted to
mp3. No other adjustments were made.

The listening tests were conducted on Amazon Turk. Each listening test was set
as a different task, and there were 3 assignments per task (this means 3
different workers were required to complete each listening test). 767 workers
completed our 11,520 listening tests, taking a median of 69 seconds to
complete.

###  A.3 Comparison of Algorithms

Algorithm | Trajectories | Memory | Time Taken (s) | Survived (%)
---|---|---|---|---
ConditionedPF-30 (PF1) | 30 | 30 | 157 | 0.44
ConditionedPF-100 (PF2) | 100 | 100 | 543 | 0.83
ConditionedPF-300 (PF3) | 300 | 300 | 1829 | 0.84
BeamSearch-30-10 (BS1) | 300 | 30 | 1057 | 0.73
BeamSearch-30-30 (BS2) | 900 | 30 | 3344 | 0.74
Table 4: Results for each generated category

Table 4 shows a summary of the performance of the algorithms mentioned in
Section 5. Trajectories indicates how many different subsequences between each
two constraints are sampled, $b\times f$ for BeamSeach and $s$ for
ConditionedPF. It can be considered as a measure of the search space of the
algorithm. Memory is the number of trajectories kept when a constraint is
reached, $f$ for BeamSearch and $s$ for ConditionedPF. Thus it is the number
promising sequences the algorithm keeps. Time taken is the average time for
each algorithm to make a sample. The parameters $b$, $f$ and $s$ for the last
5 categories were chosen to give a good comparison between BeamSearch and
ConditionedPF with similar Trajectory and Memory values, within a reasonable
amount of time (an hour). Experiments were run on an i7-6500 Intel CPU.

Each algorithm also has a survival percentage. This is the proportion of times
the algorithm returned a complete conditioned music sample rather than
failing. A failure happens when there is no choice but to sample a zero
probability event on all trajectories. This is caused by the forced event of
the constraints (the $\min$ operator returning $z_{i}$ on line 8 in Algorithm
2 or line 7 in Algorithm 1) being too unnatural to make any other
continuations likely. Thus the survival rate can be considered as how likely
the algorithm will avoid impossible sequences and find probable sequences
subject to the constraints.

Time taken is most dependent on the number of trajectories and then on memory.
The BeamSearch algorithms have the most trajectories and take $\approx$ 20mins
and $\approx$ 1 hour respectively, while having a survival rate of at most
74%. However, ConditionedPF-100 takes just under 10mins with a much higher
rate of survival. It is interesting to note that changing from 30 particles to
100 particles in ConditionedPF increases the survival rate greatly, from 44%
to 83%. Having 300 particles only increased the survival rate marginally, to
84%.