# High-energy X-ray spectrum reconstruction: solving the inverse problem from
optimized multi-material transmission measurements

Arthur Walker
Commissariat à l’Energie Atomique,
Chemin du Ru,
91680 Bruyères-le-Châtel
<EMAIL_ADDRESS>
& Alexandre Friou
Commissariat à l’Energie Atomique,
Chemin du Ru,
91680 Bruyères-le-Châtel
<EMAIL_ADDRESS>
& Kévin Ginsburger
Commissariat à l’Energie Atomique,
Chemin du Ru,
91680 Bruyères-le-Châtel
<EMAIL_ADDRESS>

###### Abstract

Reconstructing the unknown spectrum of a given X-ray source is a common
problem in a wide range of X-ray imaging tasks. For high-energy sources,
transmission measurements are mostly used to recover the X-ray spectrum, as a
solution to an inverse problem. While this inverse problem is usually under-
determined, ill-posedness can be reduced by improving the choice of
transmission measurements. A recently proposed approach optimizes custom
thicknesses of calibration materials used to generate transmission
measurements, employing a genetic algorithm to minimize the condition number
of the system matrix before inversion.

In this paper, we generalize the proposed approach to multiple calibration
materials and show a much larger decrease of the condition number of the
system matrix than thickness-only optimization. Additionally, the spectrum
reconstruction pipeline is tested in a simulation study with a challenging
high-energy Bremsstrahlung X-ray source encountered in Linear Induction
Accelerators (LIA), with strong scatter noise. Using this approach, a
realistic noise level is obtained on measurements. A generic anti-scatter grid
is designed to reduce noise to an acceptable -yet still high- noise range. A
novel noise-robust reconstruction method is then presented, which shows much
less sensitive to initialization than common expectation-maximization
approaches, enables a precise choice of spectrum resolution and a controlled
injection of prior knowledge of the X-ray spectrum.

_Keywords_ spectrum estimation $\cdot$ flash X-ray $\cdot$ transmission
measurement $\cdot$ genetic algorithm $\cdot$ ill-posed inverse problem

##  1 INTRODUCTION

The reconstruction of the unknown spectrum of a given X-ray source is a common
problem in a wide range of X-ray imaging tasks (Sidky et al., 2004; Duan et
al., 2011). If the source flux is low, spectrometers are usually preferred to
estimate the spectrum. If a very precise modeling of the source is available,
a good knowledge of the X-ray spectrum can be obtained, provided that the
model input parameters, such as voltage, are measured precisely enough during
the pulse (Wood, 2018). When none of the two previous conditions are met,
transmission measurements are most frequently used to recover the X-ray
spectrum, as a solution to an inverse problem (Waggener et al., 1999;
Armbruster et al., 2004; Paniak and Charland, 2005; Sidky et al., 2005).

This inverse problem is usually under-determined, because a high resolution of
the reconstructed spectrum is required. It is also ill-conditioned, making the
spectral estimation unstable and very sensitive to noise. The ill-posedness of
this inverse problem can be reduced using a parametric model of the
reconstructed spectrum (Zhao et al., 2015; FitzGerald et al., 2021). However,
these model-based methods restrict, by design, the space of possible
solutions, thus requiring a fine and general enough physical modelling of the
spectrum prior to reconstruction.

Another way to reduce ill-posedness is to improve the quality of the set of
transmission measurements. In particular, the recent approach described in (Li
et al., 2021) proposed to compute custom thicknesses of calibration materials
used to generate transmission measurements, by optimizing on the condition
number of the system matrix used for inversion. Using a genetic algorithm, the
interest of optimized measurements to reconstruct spectra was demonstrated, in
comparison with common linear slab phantoms.

While the approach proved to be very efficient for the configurations tested
in (Li et al., 2021), their simulation study was restricted to unrealistically
small amounts of Poisson noise, and relatively low-energy spectra. In this
work, a challenging high-energy Schiff spectrum (Schiff, 1951) is used to
simulate transmission measurements using the Monte-Carlo N-Particle code
(MCNP4C) with realistic noise levels. The Schiff spectrum is a typical model
for thin-target Bremsstrahlung spectra encountered in radiographic sources
based on Linear Induction Accelerators (LIA), used to perform high-energy
flash X-ray imaging (Estre et al., 2013). We show that the method presented in
(Li et al., 2021) is not readily applicable to this real-world reconstruction
problem. Three improvements are thus proposed to obtain a robust spectrum
reconstruction. Firstly, the transmission measurement optimization, reduced to
variations of material thicknesses in (Li et al., 2021), is extended to
multiple materials by modifying the genetic algorithm. We show that using
multiple materials yields a much larger decrease of the condition number of
the system matrix than thickness-only optimization. Secondly, instead of a
generic Poisson noise, a realistic simulated scatter noise is used to evaluate
the ability to recover the spectrum. Using this approach, we observe a much
higher noise level on measurements, which makes the design of a noise
reduction setup mandatory. As such, an anti-scatter grid is proposed, reducing
noise to an acceptable range for spectrum reconstruction, but still much
higher than noise levels encountered in (Li et al., 2021). Thirdly, once
measurements are optimized and the experimental setup is fixed, a novel
reconstruction method is presented, which shows much less sensitive to
initialization than expectation-maximization approaches (Li et al., 2021),
enables a precise choice of spectrum resolution and a controlled injection of
prior knowledge of the X-ray spectrum.

##  2 METHODS

###  2.1 Transmission measurement model

As illustrated in figure 1, two types of transmission measurements are
considered, which correspond to the two experimental configurations tested in
this study.

Figure 1: Illustration of the two configuration considered for 4 measurement
slabs

In what follows, the detector spectral response $D(E)$ is not accounted for
and taken as unity. For each measurement, the forward transmission model is
given by the Beer-Lambert law. Given an X-ray source with spectrum $S(E)$, the
transmitted intensity $I$ of a measurement writes

$I=\int_{E}S(E)D(E)e^{-l\mu(E)}\,dE$ (1)

where $l$ is the thickness and $\mu(E)$ the linear attenuation coefficient of
the material used for this measurement.

####  2.1.1 One material per measurement

This configuration corresponds to the right side of figure 1. For $M$
measurements $y_{i}$ and with an even discretization of the spectrum into $N$
bins, the forward problem can be cast into a set of linear equations

$y_{i}=\sum_{j=1}^{N}e^{-l_{i}\mu_{i,j}}s_{j},\;i\in[\\![1,M]\\!],\;j\in[\\![1,N]\\!]$
(2)

where $l_{i}$ is the thickness of the $i$-th measurement and $\mu_{i,j}$ is
the linear attenuation coefficient value for the $i$-th measurement at the
discrete energy level $j$.

Equation 2 is conveniently put in the matrix form as

$\textbf{y}=A\times\textbf{s}$ (3)

with $\textbf{y}=(y_{i})\in\mathbb{R}^{M}$ the vector of transmission
measurements, $\textbf{s}=(s_{j})\in\mathbb{R}^{N}$ the vector of the
discretized spectrum and $A=(e^{-l_{i}\mu_{i,j}})\in\mathbb{R}^{M\times N}$
the forward system matrix.

####  2.1.2 Multiple materials per measurement

This configuration is illustrated at the left side of figure 1. Each
measurement contains $K$ layers of distinct materials. Without any loss of
generality, we consider that each measurement contains the same number of
layers with the same material order. Only the thickness of each material layer
is changed, and can be set to zero. With this assumption, the attenuation
coefficient no longer depends on the measurement number $i$. The transmitted
intensity $y_{i}$ for the $i$-th measurement thus writes

$y_{i}=\sum_{j=1}^{N}\prod_{k=1}^{K}e^{-l_{i,k}\mu_{k,j}}s_{j},\;i\in[\\![1,M]\\!],\;j\in[\\![1,N]\\!],\;k\in[\\![1,K]\\!]$
(4)

where $l_{i,k}$ is the thickness of the $k$-th layer of the $i$-th measurement
and $\mu_{k,j}(E)$ is the linear attenuation coefficient of the material $k$
at energy level $j$. The forward system matrix in equation 3 is modified as
$A=(e^{-l_{i,k}\mu_{k,j}})\in\mathbb{R}^{M\times N}$.

###  2.2 Multi-material thickness optimization

Using multiple materials in transmission measurements implies some
modifications to the thickness optimization genetic algorithm described in (Li
et al., 2021). Two optimizations with slightly different implementations,
corresponding to the two types of transmission measurements discussed above,
are here considered and tested.

####  2.2.1 Multiple materials per measurement

When multiple piled up materials are allowed for each measurement, the genetic
algorithm optimizes on two-dimensional matrices of size $M\times K$ instead of
one-dimensional vectors in the single material case. Each cell of the matrices
contains the current thickness of the column’s corresponding material for the
line’s corresponding measurement. In comparison with (Li et al., 2021), the
main steps of the genetic optimization are modified as follows:

  * •

Initialization : $N_{pop}$ matrices are initialized with a randomly chosen
value in each cell, such that the sum on each row (i.e. the total thickness of
the corresponding measurement slab) remains in a chosen interval
$[l_{min},l_{max}]$

  * •

Crossover : the crossover formula used in (Li et al., 2021) is applied
similarly to the rows of the system matrix

  * •

Mutation : instead of simply modifying the value of the mutated cell as in (Li
et al., 2021), the sum of values on the row, corresponding to the total
thickness, is modified as well as the proportion of materials

####  2.2.2 One material per measurement

When only one material is allowed per measurement, two-dimensional $M\times K$
matrices are also used but with only one non-zero value per row, at the column
corresponding to the employed material. The main steps of the genetic
optimization are modified as follows:

  * •

Initialization : $N_{pop}$ matrices are initialized with a randomly chosen
value in only one cell per row, and the other cells are initialized to zero

  * •

Crossover : If the two parents are using the same material, the situation
refers to the case of (Li et al., 2021). Otherwise, the two children receive
one material each, with thickness values crossed with the same formula

  * •

Mutation : Each mutation on a row also has a probability of changing the
material used for the corresponding measurement

###  2.3 Geometry design and simulation

A basic MCNP4C geometry was built for each set of measurements. Transmission
measurement phantoms were represented as cylinders of known thicknesses and
radii. Cylinders axes are aiming at the photon source, which is modeled as a
point source located 180 cm away from the phantoms. A Schiff spectrum,
corresponding to a maximum electron energy of 20 MeV, incident on a 1.2 mm
thick tantalum target, is used. The photons are emitted within a cone with a
constant angular distribution, chosen in order to cover all the measurement
devices. Both photon and electron interactions are modeled, so as to
accurately account for scattering. The filling medium is air. All these
parameters were chosen to be as close as possible to reality.

The detectors in the simulation are modeled as simple fluence tallies, placed
behind each measurement slab. Consistent modeling of detectors in the
simulation and accounting for their spectral responses in the optimization
process, is left to future work.

As mentioned earlier, the MCNP4C simulation accounts for the presence of
scatter noise in the measurements, produced during the passage of X-rays
through materials. A specific simulation setup was designed to reduce this
scatter noise and obtain clean enough measurements for the spectrum
reconstruction.

The first approach employed to decrease scatter noise was to optimize the
placement of the measurement slabs within the experimental setup and space
them out in order to reduce the interaction between particles coming through
the materials. The differential evolution algorithm Storn and Price (1997) was
used to compute the optimal placement of $M$ points constrained in a circle by
minimizing the electrostatic potential between them. With this method, the
scatter noise was reduced by half compared to slabs placed on a circle. While
significant, this reduction was not sufficient to exploit transmission
measurements for spectrum reconstruction. To further reduce this noise, a
specific anti-scatter grid was designed. This grid consists of a 20 cm deep
lead layer between the measurement slabs and the dose sensors. Cylinders of
small radius are carved in the lead behind each measurement slab in order to
absorb all particles except source photons, yielding the expected signal.
Using this anti-scatter grid reduces the scatter noise to an exploitable
amount of around $3\%$, allowing the measurements to be used for the spectrum
reconstruction.

###  2.4 Reconstruction algorithm

####  2.4.1 Adaptive resampling based on a typical spectrum

Prior to the reconstruction algorithm itself, a preliminary step of dimension
reduction is realized on a typical spectrum presenting roughly the same
characteristics as the unknown spectrum. As shown in figure 2, by sampling
uniformly from the cumulative integral of this typical spectrum’s derivative,
an approximately optimal choice of the $N$ energy bins is obtained, allowing
an efficient spectrum representation, later employed during the optimization
process of the reconstruction.

Figure 2: Typical Schiff spectrum optimally discretized to $N=100$ energy
points

####  2.4.2 Spectrum reconstruction as an interpolation

The spectrum is reconstructed on the optimal energy sampling presented above,
using experimental or simulated measurements y. A candidate spectrum
$\textbf{s}_{cand}$ is initialized and then modified during the optimization
process. The measurements computed through the forward model
$\textbf{y}_{cand}=A\times\textbf{s}_{cand}$ are expected to be as close to y
as possible.

The optimization problem for the spectrum reconstruction thus writes :

$\underset{\textbf{s}_{cand}}{\mathrm{argmin}}\quad\lVert\textbf{y}-A\times\textbf{s}_{cand}\rVert^{2}$
(5)

A straightforward method would consist in performing a spectrum discretization
over the optimal energy sampling values $E_{1},...,E_{N}$ where $N$ is usually
large to obtain a good resolution of the spectrum. However, a large $N$
increases the degrees of freedom and makes the optimization process harder. It
is therefore necessary to reduce the number of parameters used to describe the
spectrum. The method chosen here is to sample the spectrum with a small number
$P$ of interpolation points. This requires the spectrum to be continuous and
deprived of high variation peaks, which is a general feature of high-energy
Bremsstrahlung spectra.

As shown in figure 3, the spectrum is thus fully described by $P$
interpolation points, which coordinates can evolve both in energy and
magnitude (e.g. $2P$ degrees of freedom). At each optimization step, the
points are interpolated by a piecewise cubic polynomial function Akima (1970),
then projected over the $N$ energy intervals. $\textbf{y}_{cand}$ is then
computed and new values of coordinates for the interpolation points can be
determined.

Figure 3: Schiff spectrum interpolated by P = 12 interpolation points
distributed in the energy interval

####  2.4.3 Reconstruction algorithm

A trust-region algorithm for constrained optimization Conn et al. (2000) is
used to compute the interpolation points corresponding to the approximated
spectrum. Minimum and maximum energy levels are enforced, based on a known low
energy cut-off and the maximum energy of electrons. The spectrum magnitude is
constrained using a simple step function. In addition to this restriction, the
spectrum is normalized by its integral during each step of the optimization
process.

For this algorithm, the $P$ abscissa of the interpolation points are fixed,
which improves optimization performances and decreases execution time.
However, other algorithms might be more efficient with the abscissa taken as
additional degrees of freedom. In the case of this study, these abscissa
values are taken as evenly spaced inside the energy interval, in logarithmic
scale. The optimization algorithm is then applied to the $P$ magnitudes of the
interpolation points, with a fixed number of iterations.

##  3 RESULTS

###  3.1 Impact of multiple materials

A comparison between our multi-material genetic algorithm and the version
implemented in (Li et al., 2021) was performed using the same setup. A number
of $M=8$ measurements was fixed, $N=100$ energy bins and $K=4$ materials
(iron, copper, tantalum and lead) were used. For a fair comparison, all other
hyper-parameters, including the number of generations, the population size,
the crossover and mutation probabilities, the target distribution mean and the
best fitness boosting factor, were kept identical to (Li et al., 2021).

####  3.1.1 Algorithms performance

The work done in (Li et al., 2021) led to a great reduction in orders of
magnitude of the forward matrix condition number, and highlighted an
exponential distribution of thicknesses for the optimal measurement set using
one material. Using this one-material genetic algorithm, the optimal
arrangement of thicknesses plotted on figure 4(a) was obtained, with a
condition number of $1.8\times 10^{6}$.

Figure 4(b) displays the results obtained using the first type of transmission
measurements with several materials, described in section 2.2.1. When multiple
piled up materials are allowed per measurement, a significant improvement of
the previous results is observed, with a further reduced condition number of
$1.61\times 10^{4}$.

Figure 4(c) shows the results obtained using a single material per
measurement. With this second version of the algorithm presented in 2.2.2, an
even greater improvement is observed with a highly reduced condition number of
$2.11\times 10^{3}$.

(a) Original single-material version of the genetic algorithm

(b) Version with multiple materials allowed per measurement, graphic
representation of the optimized slab for measurement 8

(c) Version with only one material allowed per measurement, ordered by
material

Figure 4: Optimal arrangement of thicknesses and achieved condition number for
the different versions of the genetic algorithm

####  3.1.2 Stability of the optimization

Regardless of the performance of the different algorithms, another main goal
of the proposed algorithms is to ensure a good stability of the optimized
measurements. To compare the stability of the multi-material algorithms
developed in this study, both were executed several times with the same
parameters. The set of materials obtained for each optimization was then
plotted.

(a) Mean optimal thicknesses and standard deviation for 20 optimizations with
the same parameters when multiple materials are allowed per measurement

(b) Optimal arrangement of thicknesses for 4 optimizations with similar
parameters when only one material is allowed per measurement

Figure 5: Stability study for both versions of the multi-material genetic
algorithm

Results are shown in figure 5. In figure 5(a) where multiple materials per
measurement are allowed, the mean thickness value and the standard deviation
of each piled up material for each measurement are plotted. However, in figure
5(b) where only one material is allowed per measurement, the plots have to be
separated because a given measurement can correspond to different materials
for different optimizations.

The result of the stability comparison is clear. As illustrated on figure
5(a), the algorithm version with piled up materials is quite unstable, with a
significant standard deviation on material thicknesses between experiments and
a strong variance on the condition number. Conversely, for the algorithm with
only one material per measurement, the stability is satisfying. As shown in
figure 5(b), each measurement is attributed almost always the same material
with the same thickness across experiments. The reduction of the research
space size thus plays a key role in the stability improvement.

###  3.2 Reduction of scatter noise

Different geometric configurations have been tested in MCNP4C to reduce
scatter noise while keeping an easy-to-design setup. For each of them, a
simulation was run using the measurement slabs returned by the one-material-
per-measurement version of the genetic algorithm for $M=12$ measurements. In
this study we consider that scattered rays represent the entire measurement
noise, which is a reasonable approximation. The Monte-Carlo simulation code
allows the calculation of the theoretical unscattered rays along with the
measurement of total (unscattered and scattered) rays. The objective is to
design a configuration in which the measured total rays are as close to the
theoretical direct rays as possible. Total and unscattered rays have been
measured for each simulation and are plotted in figure 6.

Three configurations were tested : a straightforward design where measurement
slabs were placed on the edge of a circle of given radius (config. 1, figure
6(a)), a reworked configuration in which slabs were optimally distributed in a
circle (config. 2, figure 6(b)), and a last design where a thick lead anti-
scatter grid was added between the measurement slabs and the sensors of
configuration 2 (config. 3, figure 6(c)).

(a) Case of naive configuration (config. 1)

(b) Case of optimally distributed slabs (config. 2)

(c) Case of optimally distributed slabs with a lead layer (config. 3)

Figure 6: Total ($\textrm{y}_{tot}$) and direct ($\textrm{y}_{dir}$) measured
rays for the 3 tested experimental configurations

###  3.3 Schiff spectrum reconstruction

####  3.3.1 Nominal configuration

Figure 7: Spectrum reconstruction in the nominal configuration

In this study, the nominal configuration for the spectrum reconstruction
consists in : $M=12$ measurements, $N=100$ energy intervals, $K=4$ different
materials (iron, copper, tantalum and lead), only one material allowed per
measurement in the genetic measurement set optimization (hyperparameters : 500
generations, 1000 in population size), config. 3 for the simulation
configuration, and $P=10$ interpolation points for the reconstruction
algorithm.

For this nominal configuration, the reconstruction algorithm has been applied
multiple times. Reconstructed spectrums (dotted lines) are plotted on figure
7, along with the objective theoretical spectrum (solid black line).

####  3.3.2 Ablation studies

Finally, ablation studies were performed to evaluate the influence of every
part of the spectrum estimation pipeline on the reconstruction accuracy.
Figure 8 highlights the relative improvements obtained with each major step of
the reconstruction method.

More precisely, reconstructions have been performed independently in the
nominal configuration in the cases where :

  * •

No anti-scatter grid was used for the measurements (figure 8(a))

  * •

Only one material (iron) was used for the experimental slabs (figure 8(b))

  * •

The expectation-maximization (EM) algorithm was used for the reconstruction
instead of the algorithm proposed in this study (figure 8(c))

(a) Reconstruction with no anti-scatter grid used

(b) Reconstruction with only one material

(c) Reconstrucion with expectation-maximization algorithm (EM reconstruction)
compared to the algorithm described above (Reconstruction)

Figure 8: Ablation study of the spectrum reconstruction in nominal
configuration

##  4 DISCUSSION

###  4.1 Impact of multiple materials

The work done in (Li et al., 2021) led to a great reduction in orders of
magnitude of the forward matrix condition number, and highlighted an
exponential distribution of thicknesses for the optimal measurement set using
one material. Results shown in section 3.1 illustrate the impact of using
multiple materials on the condition number, with two distinct cases.

Using multiple piled up materials is beneficial in comparison to the previous
setup from (Li et al., 2021). This can be explained by the much larger size of
the research space when multiple materials are allowed, leading to a better
optimum than the single material algorithm. However, as shown in figure 5(a),
the expected drawback of this larger research space is a worse stability. The
reduction of the research space size thus plays a key role in the stability
improvement.

With the second version of our algorithm presented in 1, where only one
material is allowed per measurement, an even greater improvement of results is
observed with a highly reduced condition number and a good stability, shown in
figures 4(c), 5(b). Even though the research space is smaller than in the
previous version of the algorithm (which theoretically includes this one’s
particular case), the relevant constraints put on the optimizer allow a more
thorough exploration, leading to the discovery of a better optimum.

###  4.2 Reduction of scatter noise

Reduction of the measurement scatter noise has proven to be necessary when a
challenging high-energy spectrum is reconstructed. Indeed, even though the
ill-posedness of the inverse problem is reduced by the search of optimal
measurement sets, the condition number remains significant enough to disturb
the reconstruction when the measurement noise is too high.

As displayed in figure 6(a), when a naive simulation configuration is used
(config. 1) the scatter noise tends to become as large as $100\%$, leading to
measurements unusable for reconstruction.

When the measurement slabs are optimally distributed within the experimental
circle (config. 2), figure 6(b) shows a substantial reduction of scattering,
with unscattered rays noised by an amount of $50\%$. This scatter noise is
however still way too high for the inversion.

Lastly, when a thick lead anti-scatter grid is added behind the measurement
slabs (config. 3), the scatter noise is reduced to a negligible amount of
around $3\%$ as shown in figure 6(c). The lead layer allows the absorption of
almost every scattered ray and the detection of the unscattered rays that pass
through material slabs parallel to its axis. The scatter noise obtained with
this last configuration appears to be small enough for the measurements to be
used to solve the inverse problem and reconstruct the spectrum.

###  4.3 Schiff spectrum reconstruction

####  4.3.1 Nominal configuration

As shown in figure 7, the overall reconstruction accuracy is satisfying in the
nominal configuration, but two areas can be distinguished. For high energy
($E>100$ keV) there is no restriction as the constraint step function is set
to unity, and the reconstruction shows great accuracy even for as few as 12
interpolation points. However, in the low energy range, the points are
constrained to 0 and are thus not optimized. In practice, this is necessary
because of the much lower contribution of low energies in the transmission
measurements. Indeed, when dense materials are subjected to an X-ray beam of
given spectrum, most of the low-energy photons are absorbed and are not
detected at the end. The difficulty to reconstruct the low energy part of
spectra is thus an intrinsic issue for the inverse problem at hand.

####  4.3.2 Ablation studies

Figure 8(a) emphasizes the substantial impact of the presence of an anti-
scatter grid in the experimental setup on the reconstruction accuracy. When
nothing is done to reduce the scattered rays, the measurements are very noisy,
leading inevitably to an inaccurate reconstruction.

Similarly, the influence of using multiple materials is illustrated in figure
8(b). Because only iron is allowed in the first optimization problem, the
genetic algorithm cannot converge to a satisfying minimum of the condition
number of the system. Thus, even with low noise, the reconstruction is
unstable and inaccurate.

Finally, the interest of using the "trust-constr" algorithm for the second
optimization problem formulated above is clearly highlighted in figure 8(c).
When the expectation-maximization algorithm is used instead, as it has usually
been done in other studies Li et al. (2021), the quality of the reconstruction
is very poor for low and medium energy ranges, and the spectrum peak is not
retrieved at the expected energy level.

##  5 CONCLUSION

In this article, a full spectrum reconstruction pipeline for high-energy X-ray
sources was presented. Building on prior work which introduced the
optimization of transmission measurements, the present work generalizes this
approach to multiple calibration materials, enabling to reach better
performance than thickness-only optimization. This work also demonstrated the
importance of noise reduction to perform spectrum reconstruction in realistic
experimental setups. As such, it showed that the design of an adapted anti-
scatter grid is a precious asset to solve the inverse problem and obtain
faithful spectrum estimations. Finally, a novel noise-robust reconstruction
method was shown to outperform common expectation-maximization approaches,
enabling a precise choice of spectrum resolution and a controlled injection of
prior knowledge of the X-ray spectrum.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

KG and AF lead the research. AW wrote the code and conducted numerical
experiments. AW, KG and AF wrote the article.

## Funding

No funding information applicable.

## References

  * Sidky et al. [2004] Emil Y. Sidky, Lifeng Yu, and Xiaochuan Pan.  Application of expectation maximization to x-ray spectrum estimation for medical accelerators from transmission data.  page 856, San Diego, CA, May 2004.  doi:10.1117/12.535989.  URL http://proceedings.spiedigitallibrary.org/proceeding.aspx?doi=10.1117/12.535989.
  * Duan et al. [2011] Xinhui Duan, Jia Wang, Lifeng Yu, Shuai Leng, and Cynthia H. McCollough.  CT scanner x-ray spectrum estimation from transmission measurements: CT scanner x-ray spectrum estimation from transmission measurements.  _Medical Physics_ , 38(2):993–997, January 2011\.  ISSN 00942405.  doi:10.1118/1.3547718.  URL http://doi.wiley.com/10.1118/1.3547718.
  * Wood [2018] Wm M. Wood.  Shot-by-shot spectrum model for rod-pinch, pulsed radiography machines.  _AIP Advances_ , 8(2):025105, February 2018.  ISSN 2158-3226.  doi:10.1063/1.5016299.  URL http://aip.scitation.org/doi/10.1063/1.5016299.
  * Waggener et al. [1999] Robert G. Waggener, Melissa M. Blough, James A. Terry, Di Chen, Nina E. Lee, Sean Zhang, and William D. McDavid.  X-ray spectra estimation using attenuation measurements from 25 kVp to 18 MV.  _Medical Physics_ , 26(7):1269–1278, July 1999\.  ISSN 00942405.  doi:10.1118/1.598622.  URL http://doi.wiley.com/10.1118/1.598622.
  * Armbruster et al. [2004] Benjamin Armbruster, Russell J Hamilton, and Arthur K Kuehl.  Spectrum reconstruction from dose measurements as a linear inverse problem.  _Physics in Medicine and Biology_ , 49(22):5087–5099, November 2004.  ISSN 0031-9155, 1361-6560.  doi:10.1088/0031-9155/49/22/005.  URL https://iopscience.iop.org/article/10.1088/0031-9155/49/22/005.
  * Paniak and Charland [2005] L D Paniak and P M Charland.  Enhanced bremsstrahlung spectrum reconstruction from depth–dose gradients.  _Physics in Medicine and Biology_ , 50(14):3245–3261, July 2005.  ISSN 0031-9155, 1361-6560.  doi:10.1088/0031-9155/50/14/004.  URL https://iopscience.iop.org/article/10.1088/0031-9155/50/14/004.
  * Sidky et al. [2005] Emil Y. Sidky, Lifeng Yu, Xiaochuan Pan, Yu Zou, and Michael Vannier.  A robust method of x-ray source spectrum estimation from transmission measurements: Demonstrated on computer simulated, scatter-free transmission data.  _Journal of Applied Physics_ , 97(12), 2005.  ISSN 0021-8979.  doi:10.1063/1.1928312.  URL http://www.scopus.com/inward/record.url?scp=21644449154&partnerID=8YFLogxK.
  * Zhao et al. [2015] Wei Zhao, Kai Niu, Sebastian Schafer, and Kevin Royalty.  An indirect transmission measurement-based spectrum estimation method for computed tomography.  _Physics in Medicine and Biology_ , 60(1):339–357, January 2015.  ISSN 0031-9155, 1361-6560.  doi:10.1088/0031-9155/60/1/339.  URL https://iopscience.iop.org/article/10.1088/0031-9155/60/1/339.
  * FitzGerald et al. [2021] Paul FitzGerald, Stephen Araujo, Mingye Wu, and Bruno De Man.  Semiempirical, parameterized spectrum estimation for x-ray computed tomography.  _Medical Physics_ , 48(5):2199–2213, May 2021\.  ISSN 0094-2405, 2473-4209.  doi:10.1002/mp.14715.  URL https://onlinelibrary.wiley.com/doi/10.1002/mp.14715.
  * Li et al. [2021] Mengzhou Li, Feng-Lei Fan, Wenxiang Cong, and Ge Wang.  EM Estimation of the X-Ray Spectrum With a Genetically Optimized Step-Wedge Phantom.  _Frontiers in Physics_ , 9:678171, May 2021.  ISSN 2296-424X.  doi:10.3389/fphy.2021.678171.  URL https://www.frontiersin.org/articles/10.3389/fphy.2021.678171/full.
  * Schiff [1951] L. I. Schiff.  Energy-Angle Distribution of Thin Target Bremsstrahlung.  _Physical Review_ , 83(2):252–253, July 1951\.  ISSN 0031-899X.  doi:10.1103/PhysRev.83.252.  URL https://link.aps.org/doi/10.1103/PhysRev.83.252.
  * Estre et al. [2013] Nicolas Estre, Daniel Eck, Jean-Luc Pettier, Emmanuel Payan, Christophe Roure, and Eric Simon.  High-energy X-ray imaging applied to non destructive characterization of large nuclear waste drums.  In _2013 3rd International Conference on Advancements in Nuclear Instrumentation, Measurement Methods and their Applications (ANIMMA)_ , pages 1–6, June 2013.  doi:10.1109/ANIMMA.2013.6727987.
  * Storn and Price [1997] Rainer Storn and Kenneth Price.  Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces.  _Journal of global optimization_ , 11:341–359, 1997.
  * Akima [1970] Hiroshi Akima.  A new method of interpolation and smooth curve fitting based on local procedures.  _Journal of the ACM (JACM)_ , 17(4):589–602, 1970\.
  * Conn et al. [2000] Andrew R Conn, Nicholas IM Gould, and Philippe L Toint.  _Trust region methods_.  SIAM, 2000.


212019275374

# Uniquely-Wilf classes

Michael Albert1  Jinge Li1  Department of Computer Science, University of
Otago

(2019-4-12; 2019-9-25; 2019-10-15)

###### Abstract

Two permutations in a class are Wilf-equivalent if, for every size, $n$, the
number of permutations in the class of size $n$ containing each of them is the
same. Those infinite classes that have only one equivalence class in each size
for this relation are characterised provided either that they avoid at least
one permutation of size 3, or at least three permutations of size 4.

##  1 Introduction

In the classical study of permutation patterns, two permutations $\pi$ and
$\sigma$ are said to be _Wilf-equivalent_ if the class,
$\operatorname{Av}(\pi)$, of permutations avoiding $\pi$ and the class,
$\operatorname{Av}(\sigma)$, of permutations avoiding $\sigma$ have the same
enumeration sequence. Notational definitions are given in the next section.
Arguably, one of the main results that initiated this area of study was that
any two permutations of size 3 are Wilf-equivalent. Due to underlying
symmetries this relation could have had at most two equivalence classes on
permutations of size 3. It has since been discovered that among the
permutations of size 4 there are only three Wilf-equivalence classes while,
taking only symmetry into account, there could have been up to seven.

In this paper we are concerned with a relativised version of Wilf-equivalence.
Rather than considering all the permutations that avoid $\pi$ we consider only
those that belong to some given class ${\mathcal{C}}$ – so the original notion
just corresponds to taking ${\mathcal{C}}$ to be the class of all
permutations.

######  Definition 1.

Let ${\mathcal{C}}$ be a permutation class. Two permutations
$\pi,\sigma\in{\mathcal{C}}$ are _Wilf-equivalent relative to ${\mathcal{C}}$_
if the two classes ${\mathcal{C}}\cap\operatorname{Av}(\pi)$ and
${\mathcal{C}}\cap\operatorname{Av}(\sigma)$ have the same enumeration
sequence. In this case we write $\pi\equiv_{{\mathcal{C}}}\sigma$.

If ${\mathcal{C}}$ is clear from context we sometimes omit mentioning it.

Let $(c_{n})_{n\in\naturals}$ be the enumeration sequence of a class
${\mathcal{C}}$. If $\pi\in{\mathcal{C}}$ has size $k$ then the enumeration
sequence $(p_{n})_{n\in\naturals}$ of
${\mathcal{C}}\cap\operatorname{Av}(\pi)$ satisfies $p_{n}=c_{n}$ for $n<k$
and $p_{k}=c_{k}-1$. Therefore, if $\pi$ and $\sigma$ are Wilf-equivalent
relative to ${\mathcal{C}}$ they must have the same size.

######  Definition 2.

Let ${\mathcal{C}}$ be a permutation class. The _Wilf-sequence_ of
${\mathcal{C}}$ is the sequence $(w_{n})_{n\in\naturals}$ where $w_{n}$ is the
number of $\equiv_{\mathcal{C}}$-equivalence classes among the permutations of
${\mathcal{C}}$ of size $n$.

Obviously $w_{n}\leqslant c_{n}$ in general but it has been observed in
multiple contexts that it is frequently the case that $w_{n}=o(c_{n})$. We
refer to this phenomenon as a _Wilf-collapse_. Since it is known that the
partially ordered sets corresponding to permutation classes do not have
infinite automorphism groups (Albert et al. (2015)) this cannot occur solely
due to global symmetries of the containment relation within a pattern class.
Our particular concern in this paper is with the ultimate Wilf-collapse:

######  Definition 3.

A permutation class ${\mathcal{C}}$ is _uniquely-Wilf_ if the equivalence
relation $\equiv_{{\mathcal{C}}}$ has only one equivalence class on
permutations of size $n$ in ${\mathcal{C}}$ for all $n$.

In Albert and Li (2019) it was shown that the class
$\operatorname{Av}(312,213)$ is uniquely-Wilf. It is the aim of this paper to
explore the landscape of uniquely-Wilf classes with the ultimate (but
unrealised) goal of classifying them completely.

##  2 Preliminaries

In this section we introduce our notation, definitions and make some basic
observations. Notations chosen are somewhat standard in the field of classical
permutation patterns so we shall be quite brief. A primer on this notation can
be found in Bevan (2015) and for introductions to and/or surveys of the field
we suggest the books Bóna (2012) and Kitaev (2011) and a survey article Vatter
(2015).

We denote by ${\mathcal{S}}$ the class of all finite permutations. For a set
of permutations $X$ we denote the set of permutations in $X$ of size $k$ by
$X_{k}$ and similarly use $X_{\leq k}$ etc. The _enumeration sequence_
$(x_{n})_{n\in\mathbb{N}}$ of $X$ is defined by $x_{k}=|X_{k}|$.

The set ${\mathcal{S}}$ carries a partial order $\leq$ called _classical
pattern containment_. Specifically, $\pi\leq\sigma$ if $\sigma$, considered as
a sequence in one line notation, contains a subsequence order-isomorphic to
$\pi$. There are eight automorphisms of ${\mathcal{S}}$ as a partially ordered
set generated by the operations on sequences of: reverse, complement (subtract
each term for a permutation of size $n$ from $n+1$) and inverse. All the
questions we consider respect these automorphisms.

If $\pi\leq\sigma$ we say that $\sigma$ _involves_ $\pi$, while if
$\pi\not\leq\sigma$ we say that $\sigma$ _avoids_ $\pi$. Further we define:

$\displaystyle\operatorname{Inv}(\pi)$
$\displaystyle=\\{\sigma\in{\mathcal{S}}\,:\,\pi\leq\sigma\\}$
$\displaystyle\operatorname{Av}(\pi)$
$\displaystyle=\\{\sigma\in{\mathcal{S}}\,:\,\pi\not\leq\sigma\\}.$

When $\pi\leq\sigma$ and $|\sigma|=|\pi|+1$ we say that $\sigma$ _covers_
$\pi$.

The $\operatorname{Av}(\pi)$ notation is extended to sets of permutations with
the requirement that all the permutations in the set must be avoided. A set of
permutations is a _class_ if it is downwards-closed with respect to $\leq$ or
equivalently if it is of the form $\operatorname{Av}(X)$ for some set of
permutations $X$. As already seen we generally use Greek letters to denote
individual permutations, upper case Latin letters to denote sets of
permutations, and calligraphic upper case Latin letters to denote permutation
classes. Additionally we set ${\mathcal{I}}=\operatorname{Av}(21)$ to be the
class of all monotone-increasing permutations and
${\mathcal{D}}=\operatorname{Av}(12)$ to be the class of all monotone-
decreasing permutations.

######  Definition 4.

Let $k$ and $n$ be positive integers. A permutation class ${\mathcal{C}}$ is
_$(k,n)$ -balanced_ if, for all $\pi,\sigma\in{\mathcal{C}}_{k}$,
$|{\mathcal{C}}_{n}\cap\operatorname{Inv}(\pi)|=|{\mathcal{C}}_{n}\cap\operatorname{Inv}(\sigma)|$.

Any permutation class is $(k,n)$-balanced for any $k\geq n$ (since the two
sets in question are singletons if $k=n$ and empty otherwise) so the issue of
whether or not a class is $(k,n)$-balanced is really of interest only if
$k<n$.

######  Observation 5.

A permutation class ${\mathcal{C}}$ is uniquely-Wilf if and only if it is
$(k,n)$-balanced for all pairs of positive integers $(k,n)$.

###### Proof.

Suppose that ${\mathcal{C}}$ is uniquely-Wilf and that
$\pi,\sigma\in{\mathcal{C}}_{k}$. Then $\pi\equiv_{{\mathcal{C}}}\tau$ so the
two classes ${\mathcal{C}}\cap\operatorname{Av}(\pi)$ and
${\mathcal{C}}\cap\operatorname{Av}(\sigma)$ have the same enumeration
sequence. In particular
$|{\mathcal{C}}_{n}\cap\operatorname{Av}(\pi)|=|{\mathcal{C}}_{n}\cap\operatorname{Av}(\sigma)|$.
But
${\mathcal{C}}_{n}\cap\operatorname{Inv}(\pi)={\mathcal{C}}_{n}\setminus({\mathcal{C}}_{n}\cap\operatorname{Av}(\pi))$
and similarly for $\sigma$ so this is equivalent to
$|{\mathcal{C}}_{n}\cap\operatorname{Inv}(\pi)|=|{\mathcal{C}}_{n}\cap\operatorname{Inv}(\sigma)|$,
i.e., ${\mathcal{C}}$ is $(k,n)$-balanced.

Conversely, suppose that ${\mathcal{C}}$ is $(k,n)$-balanced for all pairs of
positive integers $(k,n)$. Let $\pi,\sigma\in{\mathcal{C}}_{k}$ be given.
Then, for any $n$,

$|\operatorname{Av}(\pi)\cap{\mathcal{C}}_{n}|=c_{n}-|\operatorname{Inv}(\pi)\cap{\mathcal{C}}_{n}|=c_{n}-|\operatorname{Inv}(\sigma)\cap{\mathcal{C}}_{n}|=|\operatorname{Av}(\sigma)\cap{\mathcal{C}}_{n}|,$

i.e., $\pi\equiv_{{\mathcal{C}}}\sigma$, hence ${\mathcal{C}}$ is uniquely-
Wilf. ∎

######  Observation 6.

Let ${\mathcal{C}}$ be a uniquely-Wilf class that contains both $12$ and $21$.
Then, for each $k$, either ${\mathcal{C}}$ contains both $12\cdots k$ and
$k(k-1)\cdots 1$ or it contains neither of those permutations. In the latter
case ${\mathcal{C}}$ is finite.

###### Proof.

If, for some $k$, ${\mathcal{C}}$ contained $12\cdots k$ but not $k(k-1)\cdots
1$ then $12\not\equiv_{{\mathcal{C}}}21$ since ${\mathcal{C}}$ contains a
permutation of size $k$ avoiding 21 but none avoiding $12$. This contradiction
establishes the first claimed result. The second part follows immediately
since, by the Erdős-Szekeres theorem (Erdős and Szekeres (1935)), a class
containing neither $12\cdots k$ nor $k(k-1)\cdots 1$ cannot contain any
permutation of size greater than $(k-1)^{2}$. ∎

Observation 6 establishes a trichotomy for uniquely-Wilf classes: they might
be finite, consist only of monotone permutations of a single type, or contain
both monotone permutations of every size. Somewhat disingenuously we consider
the first two cases as trivial ones and focus on the last. This is
disingenuous since any uniquely-Wilf class ${\mathcal{C}}$ is the union of all
the finite uniquely-Wilf classes ${\mathcal{C}}_{\leq k}$ for each positive
integer $k$. This observation is central to the experimental results described
in the next section.

##  3 Experimental results

Let ${\mathcal{F}}$ be a finite uniquely-Wilf class whose permutations are all
of size at most $n$. We are interested in determining those uniquely-Wilf
classes ${\mathcal{C}}$ with the property that
${\mathcal{C}}\cap{\mathcal{S}}_{\leq n}={\mathcal{F}}$. This condition
already puts some restrictions on ${\mathcal{C}}$ – specifically that any
permutation of size at most $n$ that does not belong to ${\mathcal{F}}$ is
avoided by all the elements of ${\mathcal{C}}$, i.e.,
${\mathcal{C}}\subseteq\operatorname{Av}({\mathcal{S}}_{\leq
n}\setminus{\mathcal{F}})$. For convenience we denote
$\operatorname{Av}({\mathcal{S}}_{\leq n}\setminus{\mathcal{F}})$ by
${\mathcal{F}}^{\uparrow}$ – this is the maximum class (in terms of set
containment) whose restriction to ${\mathcal{S}}_{\leq n}$ is ${\mathcal{F}}$.

######  Observation 7.

Let ${\mathcal{C}}$ be a permutation class and let $k$ be a positive integer.
Then ${\mathcal{C}}=\left({\mathcal{C}}_{\leq k}\right)^{\uparrow}$ if and
only if ${\mathcal{C}}=\operatorname{Av}(X)$ for some
$X\subseteq{\mathcal{S}}_{\leq k}$.

###### Proof.

This is essentially a restatement of the definition. If
${\mathcal{C}}=\left({\mathcal{C}}_{\leq k}\right)^{\uparrow}$ then we can
take $X$ to be $\operatorname{Av}({\mathcal{S}}_{\leq
k}\setminus{\mathcal{C}}_{\leq k})$. Conversely if
${\mathcal{C}}=\operatorname{Av}(X)$ for some $X\subseteq{\mathcal{S}}_{\leq
k}$ then $X\subseteq{\mathcal{S}}_{\leq k}\setminus{\mathcal{C}}_{\leq k}$ and
so $\operatorname{Av}({\mathcal{S}}_{\leq k}\setminus{\mathcal{C}}_{\leq
k})\subseteq{\mathcal{C}}$. But the reverse containment always holds so
${\mathcal{C}}=\operatorname{Av}({\mathcal{S}}_{\leq
k}\setminus{\mathcal{C}}_{\leq k})$. ∎

If ${\mathcal{F}}\subseteq{\mathcal{S}}_{\leq n}$ is a finite uniquely-Wilf
class, then an obvious question to address is:

> _Which subsets $X\subseteq{\mathcal{S}}_{n+1}$ have the property that
> ${\mathcal{F}}\cup X$ is a uniquely-Wilf class?_

Call such subsets _potential extensions_ of ${\mathcal{F}}$ to size $n+1$.
It’s quite easy to determine whether $X$ is a potential extension of
${\mathcal{F}}$ of size $n+1$. First of all $X$ must not contain any
permutations that cannot lie in any extension of ${\mathcal{F}}$, i.e.,
$X\subseteq{\mathcal{F}}^{\uparrow}$. This ensures that ${\mathcal{F}}\cup X$
is a class. Secondly, ${\mathcal{F}}\cup X$ must be $(k,n+1)$-balanced for all
$k\leq n$. Since ${\mathcal{F}}$ was already $(k,m)$ balanced for all $m\leq
n$ and ${\mathcal{F}}\cup X$ contains no permutations of size greater than
$n+1$ these two conditions are both necessary and sufficient for $X$ to be a
potential extension of ${\mathcal{F}}$.

This suggests an experimental technique to build uniquely-Wilf classes from
the bottom up. Start with a finite permutation class ${\mathcal{F}}$ that is
uniquely-Wilf and contains permutations up to size $n$. Consider binary
vectors (i.e., indicator functions)
$\mathbf{x}:\left({\mathcal{F}}^{\uparrow}\right)_{n+1}\to\\{0,1\\}$ where, as
usual we denote $\mathbf{x}(\pi)$ by $\mathbf{x}_{\pi}$ and associate
$\mathbf{x}$ with the set $X=\mathbf{x}^{-1}(1)$. Then par abus de langage we
call $\mathbf{x}$ a potential extension of ${\mathcal{F}}$ to size $n+1$ if
the associated set $X$ is. The important thing is that this construction
already captures the first condition above since we only consider permutations
from ${\mathcal{F}}^{\uparrow}$. Each of the balancing conditions that any two
permutations of the same size, say $\sigma,\tau\in{\mathcal{F}}_{k}$, have the
same number of extensions in $X$ can be expressed as a constraint:

$\sum_{\sigma\leq\pi}\mathbf{x}_{\pi}-\sum_{\tau\leq\pi}\mathbf{x}_{\pi}=0.$

Here, and in the following few equations, $\pi$ ranges over
$\left({\mathcal{F}}^{\uparrow}\right)$ (i.e., the set of permutations on
which the notation $\mathbf{x}_{\pi}$ is defined).

From a practical standpoint we might wish to change the form in which we
express these constraints. For instance we could use instead:

$\sum_{\sigma\leq\pi,\tau\not\leq\pi}\mathbf{x}_{\pi}-\sum_{\tau\leq\pi,\sigma\not\leq\pi}\mathbf{x}_{\pi}=0.$

Or, for each $2\leq k\leq n$ we could assign a target, $t_{k}$, with

$0\leq t_{k}\leq\left|\left({\mathcal{F}}^{\uparrow}\right)_{n+1}\right|$

and, for $\sigma\in{\mathcal{F}}_{k}$, use the constraint

$\sum_{\sigma\leq\pi}\mathbf{x}_{\pi}=t_{k}.$

Or we could simply look for all solutions that give $(n,n+1)$-balance
(possibly with a target) and then filter out those that do not also give
$(k,n+1)$-balance for $k<n$.

Whatever our chosen system of constraints might be it can then be passed to a
constraint solver. We used the system Minion (Gent et al. (2006)). All
solutions can be found and then (if necessary) these can be filtered to
provide all the potential extensions. Then we can proceed in the same manner
on each of these (trying to further extend to size $n+2$) until either some
resolution is achieved or we exhaust our patience or machine resources. What
sorts of resolution can we hope for?

  * •

We might find that none of the potential extensions of ${\mathcal{F}}$ contain
the monotone permutations of size $n+1$. In this case we know that there is no
infinite uniquely-Wilf class whose restrictions to permutations of size at
most $n$ is ${\mathcal{F}}$.

  * •

We might find that the only potential extension of ${\mathcal{F}}$ that
contains both monotone permutations of size $n+1$ is
${\mathcal{F}}^{\uparrow}_{n+1}$. If this continues to be the case when we
extend to $n+2$, $n+3$ etc. then we might form the opinion that
${\mathcal{F}}^{\uparrow}$ is uniquely-Wilf and is the only such class whose
restriction to permutations of size at most $n$ is ${\mathcal{F}}$. However,
this does not constitute a proof of that opinion – but it is something we can
attempt to prove.

Since we are interested only in infinite uniquely-Wilf classes containing both
12 and 21 our starting point is the finite class ${\mathcal{S}}_{\leq
2}=\\{1,12,21\\}$. The $(1,n)$-balancing condition is always trivially
satisfied and $(2,n)$-balance is satisfied whenever we include both monotone
permutations of size $n$. So our first sets of possible extensions are the
union of $\\{123,321\\}$ with each of the sixteen subsets of
$\\{132,213,231,312\\}$. Since all questions about Wilf-equivalence are
invariant under the action of the automorphism group of $({\mathcal{S}},\leq)$
we can further restrict our experiments to one representative of each orbit of
that group on sets of permutations. Our findings are as follows:

######  Experimental results.

Let ${\mathcal{F}}\subseteq{\mathcal{S}}_{\leq 3}$ be a finite class with
${\mathcal{F}}_{2}=\\{12,21\\}$. Then:

  1. 1.

If $|{\mathcal{F}}_{3}|=2$ the only infinite uniquely-Wilf class containing
${\mathcal{F}}$ is ${\mathcal{F}}^{\uparrow}={\mathcal{I}}\cup{\mathcal{D}}$.

  2. 2.

If $|{\mathcal{F}}_{3}|=3$ the only infinite uniquely-Wilf class containing
${\mathcal{F}}$ is ${\mathcal{F}}^{\uparrow}$.

  3. 3.

If $|{\mathcal{F}}_{3}|=4$ and ${\mathcal{F}}_{3}=\\{123,213,312,321\\}$ or
one of its symmetries then the only infinite uniquely-Wilf class containing
${\mathcal{F}}$ is ${\mathcal{F}}^{\uparrow}$ while, for all other choices of
${\mathcal{F}}_{3}$ with $|{\mathcal{F}}_{3}|=4$, there are no such classes.

  4. 4.

If $|{\mathcal{F}}_{3}|=5$ then there is no infinite uniquely-Wilf class
containing ${\mathcal{F}}$.

  5. 5.

If ${\mathcal{F}}_{3}={\mathcal{S}}_{3}$ and $|{\mathcal{F}}_{4}|<22$ then
there are no infinite uniquely-Wilf classes containing ${\mathcal{F}}$.

In the next section we prove the first and second results above along with the
first half of the third result. Provided that we trust our constraint solver
(which we do) and our translation of the constraints into code111This code is
freely available from the authors on request. Owing to its specialised nature
we do not feel it appropriate, at this time, to make it part of a more widely
available repository. for the constraint solver (which we have to) the second
half of the third result, the fourth result, and the fifth result require no
further proof.

Of note is that we have not completely classified infinite uniquely-Wilf
classes since the cases where ${\mathcal{F}}_{3}={\mathcal{S}}_{3}$ and
$|{\mathcal{F}}_{4}|\in\\{22,23,24\\}$ remain open – this is where exhaustion
set in as there were simply too many possible extensions to deal with using
our present ideas. While the case $|{\mathcal{F}}_{4}|=22$ could probably be
handled by our current techniques using some additional computational power it
does seem that the remaining cases are definitely out of reach using the
straight forward constraints presented here. In order to complete the
classification would seem to require either determining some new constraints
that have been overlooked, some clever ordering of the CSP that facilitates
its (non-)solution, or a new structural idea.

However, we are confident enough to present:

######  Conjecture 8.

There are no infinite uniquely-Wilf classes, ${\mathcal{C}}$, with
${\mathcal{C}}_{3}={\mathcal{S}}_{3}$.

##  4 Proofs

In this section we provide the required proofs that justify the resolutions of
the experiment above.

######  Proposition 9.

The only infinite uniquely-Wilf class, ${\mathcal{C}}$ with
${\mathcal{C}}_{3}=\\{123,321\\}$ is ${\mathcal{I}}\cup{\mathcal{D}}$.

###### Proof.

Certainly ${\mathcal{I}}\cup{\mathcal{D}}$ is uniquely-Wilf. If
${\mathcal{C}}_{3}=\\{123,321\\}$ then ${\mathcal{C}}_{\leq
3}^{\uparrow}={\mathcal{I}}\cup{\mathcal{D}}$ as any non-monotone permutation
of size 4 contains at least one non-monotone permutation of size 3. But any
infinite uniquely-Wilf class containing 12 and 21 must contain
${\mathcal{I}}\cup{\mathcal{D}}$. ∎

Up to symmetry, there is only one set of three permutations of size three
containing $\\{123,321\\}$ so the second result we need is:

######  Proposition 10.

The only infinite uniquely-Wilf class ${\mathcal{C}}$ with
${\mathcal{C}}_{3}=\\{123,132,321\\}$ is
${\mathcal{C}}=\operatorname{Av}(213,231,312)$.

###### Proof.

Let ${\mathcal{C}}=\operatorname{Av}(213,231,312)$. Since all the permutations
avoided in ${\mathcal{C}}$ are of size three,
${\mathcal{C}}={\mathcal{C}}_{\leq 3}^{\uparrow}$. To verify that
${\mathcal{C}}$ is uniquely-Wilf is easy. The permutations in ${\mathcal{C}}$
have the structure shown in the picture below:

That is, they can be represented by an increasing sequence followed by a
decreasing one. There is an obvious one to one correspondence between such
permutations and pairs $(a,b)$ of integers with $a\geq 0$ and $b\geq 1$. Here
$a$ represents the size of the increasing sequence and $b$ the size of the
decreasing one – and we always choose the maximum element to be in the
decreasing sequence. Thus, the representation of the increasing permutation of
size $n$ in this form is as the pair $(n-1,1)$. With respect to this
representation $(a,b)\leq(c,d)$ as permutations if and only if $a\leq c$ and
$b\leq d$ i.e., the containment order of permutations and the usual
coordinatewise order coincide. In particular each permutation, $(a,b)$, of
size $n$ has exactly $k+1$ extensions of size $n+k$ namely those corresponding
to pairs $(a+x,b+k-x)$ for $0\leq x\leq k$. Thus, ${\mathcal{C}}$ is uniquely-
Wilf.

Note that, in ${\mathcal{C}}$ every permutation has exactly two covers since
$(a,b)$ is covered by $(a+1,b)$ and by $(a,b+1)$.

Suppose that ${\mathcal{E}}$ were a proper uniquely-Wilf subclass of
${\mathcal{C}}$ containing all monotone permutations and with
${\mathcal{E}}_{3}={\mathcal{C}}_{3}$. Let $n+1$ be the least size of any
permutation in ${\mathcal{C}}\setminus{\mathcal{E}}$ and note that $n\geq 3$.
Since we omit one element of ${\mathcal{C}}$ from ${\mathcal{E}}$ in size
$n+1$ there is an element of ${\mathcal{E}}$ of size $n$ with only one cover
in ${\mathcal{E}}$. Therefore, as ${\mathcal{E}}$ is uniquely-Wilf, every
element of ${\mathcal{E}}$ of size $n$ must have exactly one cover of size
$n+1$. In particular, the monotone permutation $(n,1)$ must be the unique
cover of $(n-1,1)$ belonging to ${\mathcal{E}}$ and hence its other cover
$(n-1,2)\not\in{\mathcal{E}}$.

As $(n-2,2)\in{\mathcal{E}}$, but $(n-1,2)\not\in{\mathcal{E}}$, it must be
the case that $(n-2,3)\in{\mathcal{E}}$ since that is the only remaining cover
of $(n-2,2)$ in ${\mathcal{C}}$. Now if $n=3$ we already have a contradiction
since $(1,3),(0,4)\in{\mathcal{E}}$ are the two covers of
$(0,3)\in{\mathcal{E}}$. But if $n\geq 4$ then since $(n-3,4)$ also covers
$(n-3,3)$ we have $(n-3,4)\not\in{\mathcal{E}}$. Since
$(n-1,2),(n-3,4)\not\in{\mathcal{E}}$ none of their covers can lie in
${\mathcal{E}}$. But these include $(n-1,3)$ and $(n-2,4)$ which are the two
covers of $(n-2,2)$. So $(n-2,2)$ has no covers in ${\mathcal{E}}$ and we have
our contradiction (since ${\mathcal{E}}$ is infinite all monotone permutations
have one cover in ${\mathcal{E}}$). The latter part of the argument is perhaps
better understood through the following picture where filled circles represent
elements of ${\mathcal{E}}$ and open ones elements of
${\mathcal{C}}\setminus{\mathcal{E}}$:

$n$$n+1$$n+2$$(n-1,1)$$(n-2,2)$$(n-3,3)$$(n-4,4)$$(n-1,3)$$(n-2,4)$

∎

The second half of the proof above can also be adapted to characterise the
finite uniquely-Wilf classes containing 123, 321 and one other permutation of
size three. They agree with the maximum such class to some size, then include
none or half (rounded up if necessary) of the permutations in that class at
the next size, and then no further permutations.

######  Proposition 11.

Let ${\mathcal{F}}\subseteq{\mathcal{S}}_{3}$ be the finite class with
${\mathcal{F}}_{3}=\\{123,132,231,321\\}$. The only infinite uniquely-Wilf
class containing ${\mathcal{F}}$ is
${\mathcal{F}}^{\uparrow}=\operatorname{Av}(213,312)$.

###### Proof.

That $\operatorname{Av}(213,312)$ is uniquely-Wilf is proven in Albert and Li
(2019) so we provide only a short sketch of the argument here for the reader’s
convenience. Permutations in this class have the wedge-shaped structure shown
below:

The permutations of size $n$ in the class are naturally encoded by binary
sequences of length $n-1$ over the alphabet $\\{L,R\\}$ obtained by reading
the permutation from bottom to top and noting which “side” of the picture the
corresponding element belongs. The containment relation in the class is
exactly the subsequence containment relation for these words. Consider two
words $\alpha L$ and $\beta R$ where $\alpha$ and $\beta$ are words of the
same length. A word contains $\alpha L$ but none of its prefixes do if and
only if it is of the form $AR^{k}L$ where $A$ contains $\alpha$ but none of
its prefixes do and $k\geq 0$. But similarly a word contains $\beta R$ but
none of its prefixes do if and only if it is of the form $BL^{k}R$ where $B$
contains $\beta$ but none of its prefixes do. This sets up a length-preserving
bijection between the minimal words containing $\alpha L$ and those containing
$\beta R$ whenever we have such a bijection with respect to $\alpha$ and
$\beta$. But of course such a bijection exists when $\alpha$ and $\beta$ are
both the empty word, so induction shows that this bijection exists in general
for any two words of the same length. This bijection then extends to a length-
preserving bijection between all words containing $\alpha L$ and all words
containing $\beta R$ since these are obtained from words of the first type by
concatenation with an arbitrary word. This establishes that
$\operatorname{Av}(213,312)$ is uniquely-Wilf.

Now suppose that ${\mathcal{E}}$ is a proper subclass of
$\operatorname{Av}(213,312)$ that contains ${\mathcal{F}}_{3}$. The
permutations in ${\mathcal{E}}$ can be partitioned into their maximal monotone
intervals which must be alternately increasing and decreasing when considered
from bottom to top. There must be a bound on the number of such intervals
(since ${\mathcal{E}}$ is proper). Therefore the structure of a permutation in
${\mathcal{E}}$ is determined (essentially) by a finite tuple of non-negative
integers summing to its size and the enumeration of ${\mathcal{E}}$ is
polynomially bounded and hence ultimately polynomial (Kaiser and Klazar
(2002/03)).

The structure of polynomial classes is well-understood and we will make use of
a form introduced in Homberger and Vatter (2016) – we sketch only the
necessary details here and refer to that paper for the full results. Define a
_peg permutation_ to be a permutation written in one line notation where each
of its elements is decorated with a superscript from the set
$\\{+,-,\bullet\\}$. For example, $2^{-}3^{-}1^{\bullet}$ is a peg
permutation.

Given a pegged permutation $\tilde{\rho}$ the set of permutations
$\operatorname{Grid}(\tilde{\rho})$ consists of all those permutation that can
be obtained from $\rho$ by inflating each point into a monotone (possibly
empty) interval according to the associated symbol: for $+$ an arbitrary
increasing interval, for $-$ an arbitrary decreasing interval, and for
$\bullet$ either a single point or nothing at all. For instance $4321$, $2431$
and $432651$ all belong to $\operatorname{Grid}(2^{-}3^{-}1^{\bullet})$ but
$123$ does not.

Departing from the terminology of Homberger and Vatter (2016) we say that a
pegged permutation is _properly pegged_ if it contains no two element monotone
interval of any of the following forms:

$1^{+}2^{+},\,1^{\bullet}2^{+},\>1^{+}2^{\bullet},\,2^{-}1^{-},\,2^{\bullet}1^{-},\,2^{-}1^{\bullet}.$
(1)

To be explicit, the first condition above means that there are no two
consecutive elements $a$ and $a+1$ in the pegged permutation both of which are
decorated with $+$. We define $\operatorname{Grid}^{f}(\tilde{\rho})$ to be
the subset of $\operatorname{Grid}(\tilde{\rho})$ consisting of all inflations
of $\rho$ in which any symbol decorated with $+$ or $-$ is inflated into a
monotone interval of size at least two, and any symbol decorated with
$\bullet$ is replaced by a single point (the superscript stands for “filled”).
Restating the conclusion of Theorem 1.4 of Homberger and Vatter (2016)
slightly:

_There exists a finite set $\tilde{H}$ of properly pegged permutations such
that ${\mathcal{E}}$ is the disjoint union of the sets
$\operatorname{Grid}^{f}(\tilde{\rho})$ for $\tilde{\rho}\in\tilde{H}$._

We remark that since ${\mathcal{E}}$ is a class that, if
$\tilde{\rho}\in\tilde{H}$ then in fact
$\operatorname{Grid}(\tilde{\rho})\subseteq{\mathcal{E}}$ since
$\operatorname{Grid}(\tilde{\rho})$ is just the downward closure of
$\operatorname{Grid}^{f}(\tilde{\rho})$.

Choose $\tilde{\rho}\in\tilde{H}$ in such a way that the total number, $m$ of
$+$ and $-$ decoration on $\tilde{\rho}$ is maximal (among all the elements of
$\tilde{H}$) and, among all such, the underlying permutation $\rho$ is maximal
in the pattern containment order. Equality is allowed in either case. For
instance it might be the case that
$2^{+}4^{+}1^{\bullet}3^{-},2^{-}4^{\bullet}1^{-}3^{+}\in\tilde{H}$ . If no
other permutation in $\tilde{H}$ contains more than three non-$\bullet$
decorations and no permutation with three non-$\bullet$ decorations contains
2413 then either of these would be a suitable choice for $\tilde{\rho}$. Form
a permutation $\pi\in\operatorname{Grid}^{f}(\tilde{\rho})$ by inflating each
such non-$\bullet$ symbol to a monotone interval having at least $K$ elements
where $K$ is greater than the square of the largest size of any permutation
underlying an element of $\tilde{H}$. This choice of $K$ ensures (easily) that
the only covers of $\pi$ in ${\mathcal{E}}$ are obtained by extending the
inflated intervals by one point. In particular $\pi$ has exactly $m$ covers in
${\mathcal{E}}$. Therefore, as ${\mathcal{E}}$ is uniquely-Wilf, _all_
sufficiently large permutations in ${\mathcal{E}}$ must have exactly $m$
covers.

Now suppose that
$\operatorname{Grid}(1^{+}3^{+}2^{\bullet})\subseteq{\mathcal{E}}$. Then any
permutation that consists of a monotone increasing sequence followed by one
more element belongs to ${\mathcal{E}}$ as these all occur as elements of
$\operatorname{Grid}(1^{+}3^{+}2^{\bullet})$. In particular, for any $n$, the
monotone increasing permutation of size $n$ has at least $n+1$ covers in
${\mathcal{E}}$ which contradicts the observation above that all sufficiently
large permutations in ${\mathcal{E}}$ must have exactly $m$ covers. So it
follows by contradiction that
$\operatorname{Grid}(1^{+}3^{+}2^{\bullet})\not\subseteq{\mathcal{E}}$. This,
together with the fact that ${\mathcal{E}}$ is a class, implies that no
$\tilde{\rho}\in\tilde{H}$ contains a 132 pattern with both the 1 and the 3
carrying $+$ decorations – we’ll abbreviate this as “no $1^{+}3^{+}2$
pattern”. Similarly, there can be no $\tilde{\rho}\in\tilde{H}$ containing a
$23^{-}1^{-}$ pattern.

This technique of showing that certain restrictions on the labelings of the
permutations in $\tilde{H}$ must be satisfied because otherwise some
permutation has “too many” covers contradicting the assumption that
${\mathcal{E}}$ is uniquely-Wilf is used throughout the remainder of the
argument.

The underlying permutations of $\tilde{H}$ belong to
$\operatorname{Av}(312,213)$ so if there were three or more non-$\bullet$
decorations then two of them would occur on the same side and would therefore
have the same sign. Since the intervals described in (1) cannot occur they
could not be consecutive, and this would create one of the forbidden
situations above. It follows that there can be at most two non-$\bullet$
decorations on each permutation of $\tilde{H}$ and, if there are two, they
must have differing signs.

Therefore, every sufficiently large permutation in ${\mathcal{E}}$ has either
exactly one or exactly two covers. The former case is easily dismissed – since
the increasing permutation has only the increasing cover the only permutation
in $\tilde{H}$ carrying a $+$ decoration would have to be $1^{+}$ and
similarly $1^{-}$ would be the only one containing the $-$ decoration (and all
the monotone permutations of ${\mathcal{E}}$ would be represented by these
along with $1^{\bullet}$). The rest of the permutations (if any) would carry
only $\bullet$s so define a finite set – but if any existed there would a non-
monotone permutation with no covers. So
${\mathcal{E}}={\mathcal{I}}\cup{\mathcal{D}}$ but supposedly
${\mathcal{E}}_{3}={\mathcal{F}}_{3}$ which contains non-monotone
permutations.

Supposing that each sufficiently large permutation in ${\mathcal{E}}$ has
exactly two covers also implies that

$\operatorname{Grid}(2^{+}3^{-}1^{\bullet})\not\subseteq{\mathcal{E}}$

since otherwise the monotone increasing permutation of size $n$ would then
have all of:

$123\cdots(n+1),\,23\cdots(n+1)1,\,\mbox{and}\,123\cdots(n-1)(n+1)n$

as covers. Similarly,
$\operatorname{Grid}(1^{\bullet}3^{+}2^{-})\not\subseteq{\mathcal{E}}$. Also,

$\operatorname{Grid}(1^{+}3^{\bullet}4^{\bullet}2^{-})\not\subseteq{\mathcal{E}}\,\text{and}\,\operatorname{Grid}(2^{+}4^{\bullet}3^{\bullet}1^{-})\not\subseteq{\mathcal{E}}$

as otherwise the monotone increasing (or decreasing in the second case)
permutation would have at least three extensions.

Putting all these observations together shows that the only possible elements
of $\tilde{H}$ containing both a $+$ and a $-$ decoration are $1^{+}2^{-}$ and
$2^{+}1^{-}$. Only one can occur as, if they both did, increasing permutations
would have at least three covers in ${\mathcal{E}}$.

Suppose without loss of generality that $1^{+}2^{-}\in\tilde{H}$. We claim
that, other than $1^{+}$, $1^{-}$, and $1^{\bullet}2^{-}$ all other patterns
in $\tilde{H}$ contain only $\bullet$s as decorations.

To rule out other patterns containing symbols $x^{+}$ observe first that
$\operatorname{Grid}(2^{+}1^{\bullet})\not\subseteq{\mathcal{E}}$ since if it
were, there would be an extra extension to the monotone increasing permutation
beyond the two already available from $\operatorname{Grid}(1^{+}2^{-})$. It
follows, that the only possible $+$ decoration on an element of $\tilde{H}$
can occur on the first (and least) element of a pattern. No decorations of
$132$ can belong to $\tilde{H}$ since we can’t place $\bullet$s on both the 3
and the 2 as that case is already covered by $1^{+}2^{-}$ and other labelings
are forbidden by (1) or previous observations. By (1) we cannot begin
$1^{+}2^{\bullet}$. Finally,
$\operatorname{Grid}(1^{+}3^{\bullet}4^{\bullet}2^{\bullet})$ also gives extra
extensions to the monotone increasing permutation. Any larger labeled pattern
includes one of these forbidden cases, and so we’ve ruled out the presence of
a $+$ except for $1^{+}2^{-}$ and $1^{+}$.

Now we want to rule out patterns in $\tilde{H}$ other than $1^{-}$,
$1^{\bullet}2^{-}$ and $1^{+}2^{-}$ that contain symbols $x^{-}$. We cannot
appeal to symmetry because we’ve broken that by assuming
$\operatorname{Grid}(1^{+}2^{-})$ is contained in ${\mathcal{E}}$. In
$\operatorname{Grid}(2^{\bullet}3^{-}1^{\bullet})$ we gain an extra extension
to the monotone decreasing permutation so no pattern in $\tilde{H}$ can
contain a $23^{-}1$ pattern. Similarly no $231^{-}$ patterns can occur.
Together these rule out all possibilities of having a symbol $x^{-}$ other
than the ones which are explicitly allowed.

Now $231\in{\mathcal{E}}$ but is not in $\operatorname{Grid}(1^{+})$,
$\operatorname{Grid}(1^{-})$, $\operatorname{Grid}(1^{\bullet}2^{-})$, or
$\operatorname{Grid}(1^{+}2^{-})$ so its representation must be as
$\operatorname{Grid}(2^{\bullet}3^{\bullet}1^{\bullet})$. The same applies for
its supposed two covers in ${\mathcal{E}}$, each of their covers, and so on.
However, $\tilde{H}$ is finite, so we have our final contradiction. ∎

The proof of Proposition 11 has obviously been the most intricate of the
arguments so far. Fundamentally though it is not so different in character to
the previous ones (indeed it would have been possible to absorb those into its
proof). Working in a polynomial class is a powerful restriction because of the
limit it imposes on the number of extensions of any sufficiently large
permutation. This then enables the “surplus extensions” technique to rule out
many configurations. The second author, as part of his current PhD research
(to be submitted in 2020), has proven that the only uniquely-Wilf polynomial
classes are the ones that already arise in Proposition 11.

##  5 Conclusion

Our goal of characterising the uniquely-Wilf permutation classes has not been
fully realised as the possibility still exists that there may be some which we
have not discovered. If so, they will contain all of ${\mathcal{S}}_{3}$ and
at least 22 elements of ${\mathcal{S}}_{4}$.

It is notable that most of the arguments by contradiction in Proposition 11
focused on the existence of surplus extensions to a monotone permutation. The
monotone permutations are already somewhat special in the order of pattern
containment since they are the only permutations that cover only one other
permutation. Therefore it does not seem unreasonable to ask:

######  Question 12.

What can be determined about classes in which all non-monotone permutations
are Wilf-equivalent?

For polynomial classes, the second author has answered this question
completely (forthcoming as part of his PhD thesis in 2020) and, other than the
examples already discovered above there are only a few new examples for
instance the class $\operatorname{Av}(123,312)$ which was considered in Albert
and Li (2019).

The authors would like to acknowledge the helpful comments of two referees,
which have improved the presentation of this work significantly.

## References

  * Albert and Li (2019) M. Albert and J. Li.  Wilf-collapse in permutation classes having two basis elements of size three.  _Austral. J. Comb._ , 73(3):440–455, 2019.  ISSN 2202-3518.
  * Albert et al. (2015) M. H. Albert, M. D. Atkinson, and A. Claesson.  Isomorphisms between pattern classes.  _J. Comb._ , 6(1-2):3–17, 2015.  ISSN 2156-3527.  10.4310/JOC.2015.v6.n1.a1.  URL https://doi.org/10.4310/JOC.2015.v6.n1.a1.
  * Bevan (2015) D. Bevan.  Permutation patterns: basic definitions and notation.  _arXiv e-prints_ , art. arXiv:1506.06673, Jun 2015.
  * Bóna (2012) M. Bóna.  _Combinatorics of permutations_.  Discrete Mathematics and its Applications (Boca Raton). CRC Press, Boca Raton, FL, second edition, 2012.  ISBN 978-1-4398-5051-0.  10.1201/b12210.  URL https://doi.org/10.1201/b12210.  With a foreword by Richard Stanley.
  * Erdős and Szekeres (1935) P. Erdős and G. Szekeres.  A combinatorial problem in geometry.  _Compos. Math._ , 2:463–470, 1935.
  * Gent et al. (2006) I. P. Gent, C. Jefferson, and I. Miguel.  Minion: A fast, scalable, constraint solver.  In _Proceedings of the 2006 Conference on ECAI 2006: 17th European Conference on Artificial Intelligence August 29 – September 1, 2006, Riva Del Garda, Italy_ , pages 98–102, Amsterdam, The Netherlands, The Netherlands, 2006. IOS Press.  ISBN 1-58603-642-4.  URL http://dl.acm.org/citation.cfm?id=1567016.1567043.
  * Homberger and Vatter (2016) C. Homberger and V. Vatter.  On the effective and automatic enumeration of polynomial permutation classes.  _J. Symbolic Comput._ , 76:84–96, 2016.  ISSN 0747-7171.  10.1016/j.jsc.2015.11.019.  URL https://doi.org/10.1016/j.jsc.2015.11.019.
  * Kaiser and Klazar (2002/03) T. Kaiser and M. Klazar.  On growth rates of closed permutation classes.  _Electron. J. Combin._ , 9(2):Research paper 10, 20, 2002/03.  ISSN 1077-8926.  URL http://www.combinatorics.org/Volume_9/Abstracts/v9i2r10.html.  Permutation patterns (Otago, 2003).
  * Kitaev (2011) S. Kitaev.  _Patterns in permutations and words_.  Monographs in Theoretical Computer Science. An EATCS Series. Springer, Heidelberg, 2011.  ISBN 978-3-642-17332-5; 978-3-642-17333-2.  10.1007/978-3-642-17333-2.  URL https://doi.org/10.1007/978-3-642-17333-2.  With a foreword by Jeffrey B. Remmel.
  * Vatter (2015) V. Vatter.  Permutation classes.  In _Handbook of enumerative combinatorics_ , Discrete Math. Appl. (Boca Raton), pages 753–833. CRC Press, Boca Raton, FL, 2015.