%\VignetteIndexEntry{eulerr under the hood}
%\VignetteEngine{knitr::knitr}

\PassOptionsToPackage{
  usenames,
  fixpdftex,
  dvipsnames,
  svgnames,
  x11names
}{xcolor}
\PassOptionsToPackage{
  font=small,
  labelfont=bf,
  labelsep=period,
  justification=RaggedRight,
  format=plain,
  singlelinecheck=off
}{caption}
\PassOptionsToPackage{hyphens}{url}

\documentclass[
  numbers=noendperiod,
  bibliography=totoc,
  parskip=half
]{scrartcl}

% Page layout
\usepackage[
  a4paper,
  left=23mm,
  top=27.4mm,
  bottom=27.4mm,
  textwidth=107mm,
  marginparsep=8mm,
  marginparwidth=49mm,
]{geometry}
\usepackage{multicol}

% Margin paragraph
\usepackage{marginfix}

% Headers and footers
\usepackage[
  automark,
  headsepline=false,
  headwidth=textwithmarginpar,
  footwidth=head
]{scrlayer-scrpage}

\clearpairofpagestyles
\rofoot[\pagemark]{\pagemark}
\lohead{\headmark}

% Title page
\title{eulerr under the hood}
\author{Johan larsson}
\date{\today}

\usepackage{titling}
\setlength{\droptitle}{-5em}
\pretitle{\begin{flushleft}\LARGE\bfseries}
\posttitle{\end{flushleft}}
\preauthor{\begin{flushleft}\Large}
\postauthor{\end{flushleft}}
\predate{\begin{flushleft}\large\itshape}
\postdate{\end{flushleft}}

\usepackage[font=footnotesize]{subcaption}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage[document]{ragged2e}
\usepackage[shortlabels]{enumitem}
\usepackage[noend]{algpseudocode}
\usepackage{etoolbox}
%\usepackage[swedish]{babel}
\usepackage[square,numbers,sort&compress]{natbib}
\usepackage{floatrow}
\usepackage{environ}
\usepackage{xspace}

% Fonts
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage[lining]{libertine}
\usepackage{textcomp}
\usepackage[varqu,varl,scaled=0.93]{inconsolata}
\usepackage{mathtools}
\usepackage{amsthm}
\usepackage[libertine,vvarbb,libaltvw,liby]{newtxmath}
\usepackage[scr=rsfso]{mathalfa}
\usepackage{bm}
\useosf
\usepackage{microtype}

% Two-column bibliography
\patchcmd{\thebibliography}{\list}{\begin{multicols}{2}\small\list}{}{}
\appto{\endthebibliography}{\end{multicols}}

% Block margins for caption around figure and table environments
\BeforeBeginEnvironment{figure}{\blockmargin}
\AfterEndEnvironment{figure}{\unblockmargin}
\BeforeBeginEnvironment{table}{\blockmargin}
\AfterEndEnvironment{table}{\unblockmargin}

% Floats
\DeclareNewFloatType{algorithm}{fileext=alg, name=Algorithm}
\DeclareFloatSeparators{marginparsep}{\hskip\marginparsep}
\floatsetup[figure]{%
  margins=hangright,
  capposition=beside,
  capbesidesep=marginparsep,
  capbesideposition={top,right},
  floatwidth=\textwidth}
\floatsetup[table]{%
  margins=hangright,
  capposition=beside,
  capbesidesep=marginparsep,
  capbesideposition={top,right},
  floatwidth=\textwidth}
\floatsetup[widefigure]{%
  margins=hangright,
  capposition=bottom}
\floatsetup[widetable]{%
  margins=hangright,
  capposition=top}
\floatsetup[algorithm]{%
  margins=hangright,
  capposition=beside,
  capbesidesep=marginparsep,
  capbesideposition={top,right},
  floatwidth=\textwidth}
\captionsetup[subfigure]{position=b,justification=RaggedRight}
\captionsetup[sub]{justification=RaggedRight}

% Provide marginfigure
\NewEnviron{marginfigure}{%
  \marginpar{%
    \captionsetup{type=figure}
    \BODY
  }%
}

% Graphics options
\setkeys{Gin}{%
  width=\linewidth,
  totalheight=\textheight,
  keepaspectratio
}
\graphicspath{{graphics/}}

% Extended verbatim environments
\usepackage{fancyvrb}
\fvset{fontsize=\small}% default font size for fancy-verbatim environments

% Provide fullwidth environment
\newlength{\overhang}
\setlength{\overhang}{\marginparwidth}
\addtolength{\overhang}{\marginparsep}

\newenvironment{fullwidth}{%
  \blockmargin
  \begin{addmargin*}[0em]{-\overhang}%
}{%
  \end{addmargin*}%
  \unblockmargin
}

\newcommand{\proglang}[1]{\textsf{#1}}
\newcommand{\pkg}[1]{{\fontseries{b}\selectfont #1}}
\newcommand{\code}[1]{\texttt{#1}}

% Custom operators
\DeclareMathOperator{\E}{E}
\DeclareMathOperator{\Pois}{Pois}
\DeclareMathOperator{\B}{Bin}
\DeclareMathOperator{\V}{Var}
\DeclareMathOperator{\Exp}{Exp}

% Cross-referencing and colors
\usepackage{hyperref}
\usepackage[noabbrev,capitalize,nameinlink]{cleveref}
\usepackage{xcolor}
\hypersetup{%
  linkcolor=SteelBlue4,
  citecolor=SteelBlue4,
  urlcolor=SteelBlue4,
  colorlinks=true,
  breaklinks=true
}

% Theorems
\newtheorem{mydef}{Definition}

% Redefine footnotes to be marginnotes
\newcounter{snmark}
\setcounter{snmark}{0}

\newcommand*{\makesidenotemark}{\/\textsuperscript{\thesnmark}}

\renewcommand{\footnote}[1]{%
  \refstepcounter{snmark}%
  \makesidenotemark{}%
  \marginpar{\RaggedRight\small\textsuperscript{\thesnmark}\,#1}%
}

% Add cref
\newcommand{\secref}[1]{\hyperref[#1]{\ref*{#1}~\nameref*{#1}}}

% Fix vdots and ddots for smallmatrix
\newcommand{\svdots}{\raisebox{3pt}{$\scalebox{.75}{\vdots}$}}
\newcommand{\sddots}{\raisebox{3pt}{$\scalebox{.75}{$\ddots$}$}}

\let\ts\textsuperscript

\begin{document}

<<setup, include = FALSE>>=
knitr::knit_hooks$set(crop = knitr::hook_pdfcrop)
knitr::knit_hooks$set(
  crop = knitr::hook_pdfcrop,
  document = function(x) {
    sub('\\usepackage[]{color}', '\\usepackage{xcolor}', x, fixed = TRUE)
  }
)

knitr::opts_chunk$set(
  echo = FALSE,
  warning = FALSE,
  message = FALSE,
  fig.width = 2.1,
  fig.height = 2.1,
  crop = TRUE,
  fig.align = "center",
  dev = if (capabilities("cairo"))
    "cairo_pdf"
  else "pdf"
)

library(grid)
library(RConics)
library(eulerr)
library(lattice)

lattice::lattice.options(default.theme = list(fontsize = list(text = 8,
                                                              points = 6)))

pal <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442",
         "#0072B2", "#D55E00", "#CC79A7")

ellipse_arc <- function(saxes = c(1, 1),
                        loc = c(0, 0),
                        theta = 0,
                        n = 200,
                        rng = c(0, 2*pi)) {
  b <- min(saxes[1], saxes[2])
  a <- max(saxes[1], saxes[2])
  d2 <- (a - b) * (a + b)
  if (length(rng) == 1)
    phi <- rng - theta
  else
    phi <- seq(rng[1], rng[2], len = n) - theta
  sp <- sin(phi)
  cp <- cos(phi)
  r <- a * b/sqrt((saxes[2] * cp)^2 + (saxes[1] * sp)^2)
  P <- matrix(nrow = n, ncol = 2)
  P[, 1] <- r * cp
  P[, 2] <- r * sp

  if (theta != 0) {
    P <- P %*% matrix(c(cos(theta), sin(theta), -sin(theta), cos(theta)),
                      byrow = TRUE, nrow = 2, ncol = 2)
  }
  P <- P + matrix(loc[1:2], nrow = nrow(P), ncol = 2, byrow = TRUE)
  P
}

set.seed(1)
@

\maketitle

% Always avoid colored links in TOC
% {\hypersetup{linkcolor=black}
% \tableofcontents
% }

\section{Introduction}
\label{sec:introduction}

\textbf{eulerr} relies on an extensive machinery to turn user input into
a pretty Euler diagram. Little of this requires any
tinkering from the user. To make that happen, however, \textbf{eulerr}
needs to make several well-formed decisions about the design of the
diagram on behalf of the user, which is not a trivial task.

This document outlines the implementation of \textbf{eulerr} from input
to output. It is designed to be an evolving documentation on the innards of the
program.

\section{Input}
\label{sec:input}

Euler diagrams present relationships between sets, wherefore the data must
describe these relationships, either directly or indirectly.
\pkg{eulerr} allows several alternatives for this data, namely,
\begin{itemize}
\item intersections and relative
complements\footnote{\(A \setminus B = 3 \quad B \setminus A = 2 \quad A \cap B=1\)},
\item unions and identities\footnote{\(A=4 \quad B=3 \quad A \cap B=1\)},
\item a matrix of binary (or boolean) indices\footnote{%
  $ \begin{bmatrix}
      \bm{A} & \bm{B}\\
      1      & 0     \\
      1      & 0     \\
      1      & 0     \\
      1      & 1     \\
      0      & 1     \\
      0      & 1
    \end{bmatrix}$},
\item a list of sample spaces\footnote{%
  $ \begin{array}{l}
      A = \{a,\,b,\,c,\,d\}\\
      B = \{a,\,e,\,f\}
    \end{array}$}, or
\item a two- or three-way table\footnote{%
  \begin{tabular}{lcc}
    \toprule
    & $A$ & $A^\mathsf{c}$ \\
    \midrule
    $B$ & 1 & 2 \\
    $B^\mathsf{c}$ & 3 & 0 \\
    \bottomrule
  \end{tabular}}.
\end{itemize}

As an additional feature for the matrix form, the user may supply a factor
variable with which to split the data set before fitting the diagram,
which sometimes improves diagrams where the set relationships vary
across categories.

Whichever type of input is provided, \pkg{eulerr} will translate it to the first
and second types, \emph{intersections and relative complements} and
\emph{unions and identities}, which will be used in the steps to come.

The Euler diagram is then fit in two
steps: first, an initial layout is formed with circles using only the
sets' pairwise relationships. Second, this layout is fine-tuned taking
all $2^N-1$ intersections into consideration.

\section{Initial layout}
\label{sec:initConfig}

For our initial layout, we adopt a constrained version of
multi-dimensional scaling~(MDS) that has been adapted from
\pkg{venn.js}~\citep{Frederickson_2016}, which in turn is a modification of
an algorithm used in \pkg{venneuler}~\citep{Wilkinson_2012}.
In it, we consider only the pairwise intersections between sets, attempting to
position their respective shapes so as to minimize the difference between the
separation between their centers
required to obtain an optimal overlap
and the actual overlap of the shapes in the diagram.

This problem is unfortunately intractable for ellipses, being that there is an
infinite number of ways by which we can position two ellipses to obtain a given
overlap. Thus, we restrict ourselves to circles in our initial layout, for which we can use the
circle--circle overlap formula~\eqref{eq:circleOverlap} to numerically find the
required distance, $d$, for each pairwise relationship.
\begin{fullwidth}
\begin{multline}
O_{ij} = r_i^2\arccos\left(\frac{d_{ij}^2 + r_i^2 - r_j^2}{2d_{ij}r_i}\right) +
r_j^2\arccos\left(\frac{d_{ij}^2 + r_j^2 - r_i^2}{2d_{ij}r_j}\right) - \\
\frac{1}{2}\sqrt{(-d_{ij} + r_i + r_j)(d_{ij} + r_i - r_j)(d_{ij} - r_i + r_j)(d_{ij} + r_i + r_j)},
\label{eq:circleOverlap}
\end{multline}
\end{fullwidth}
where $r_i$ and $r_j$ are the radii of the circles representing the $i$\ts{th} and
$j$\ts{th} sets respectively, $O_{ij}$ their overlap, and $d_{ij}$ their separation.

\begin{marginfigure}
<<circleOverlap, fig.width=2.2, fig.height = 2.4>>=
c0 <- ellipseToConicMatrix(c(1, 1), c(0, 0), 0)
c1 <- ellipseToConicMatrix(c(0.7, 0.7), c(1.2, 0), 0)
pp <- intersectConicConic(c0, c1)

theta0 <- atan2(pp[2, 1], pp[1, 1])
theta1 <- atan2(pp[2, 2], pp[1, 2])

phi0 <- atan2(pp[2, 1], pp[1, 1] - 1.2)
phi1 <- atan2(pp[2, 2], pp[1, 2] - 1.2)

seg <- rbind(ellipse_arc(c(0.7, 0.7), c(1.2, 0), 0, n = 50, c(-pi, phi0)),
             ellipse_arc(c(1, 1), c(0, 0), 0, n = 100, c(theta0, theta1)),
             ellipse_arc(c(0.7, 0.7), c(1.2, 0), 0, n = 50, c(phi1, pi)))

ospot <- c(mean(seg[, 1]), mean(seg[, 2]))

xyplot(1~1, xlim = c(-1.1, 2), ylim = c(-1.5, 1.1), asp = "iso",
       xlab = NULL, ylab = NULL,
       scales = list(draw = FALSE),
       par.settings = list(axis.line = list(col = "transparent")),
       panel = function() {
         panel.polygon(seg, col = "steelblue1", alpha = 0.5,
                       border = "transparent")
         grid::grid.circle(0, 0, r = 1, default.units = "native",
                           gp = gpar(fill = "transparent"))
         grid::grid.circle(1.2, 0, r = 0.7, default.units = "native",
                           gp = gpar(fill = "transparent"))
         pBrackets::grid.brackets(1.2, -1.2, 0, -1.2, h = 0.05)
         pBrackets::grid.brackets(-1, 0, 0, 0, h = 0.05)
         pBrackets::grid.brackets(1.2, 0, 1.9, 0, h = 0.05)
         panel.lines(c(0, 0), c(0, -1.2), lty = 2, col = "grey65")
         panel.lines(c(1.2, 1.2), c(0, -1.2), lty = 2, col = "grey65")
         panel.text(0.6, unit(-1.3, "native"),
                    labels = expression(italic(d[ij])), pos = 1)
         panel.text(ospot[1], ospot[2], labels = expression(italic(O[ij])))
         panel.text(-0.5, 0.1, labels = expression(italic(r[i])), pos = 3,
                    default.units = "native")
         panel.text(1.55, 0.1, labels = expression(italic(r[j])), pos = 3,
                    default.units = "native")
         panel.points(c(0, 1.2), c(0, 0), pch = 21, col = 1, cex = 1,
                      fill = "white")
       })
@
\caption{The circle--circle overlap is computed as a function of the discs'
separation ($d_{ij}$), radii ($r_i,r_j$), and area of overlap ($O_{ij}$).}
\label{fig:circleCircle}
\end{marginfigure}

Setting $r_i = \sqrt{F_i/\pi}$, where $F_i$ is
the size of the $i$\ts{th} set, we are able to obtain $d$ numerically using
the squared difference between $O$ and the desired overlap
as loss function~\eqref{eq:dopt},
\begin{equation}
  \mathcal{L}(d_{ij}) = \left(O_{ij} - (F_i \cap F_j)  \right)^2, \quad \text{for } i <
    j \leq n,
\label{eq:dopt}
\end{equation}
which we optimize using \code{optimize()}\footnote{According to the
documentation, \code{optimize()} consists of a ``combination of golden section
search and successive parabolic interpolation.''} from \pkg{stats}.

For a two-set combination, this is all we need to plot an exact diagram, given
that we now have the two circles' radii and separation and may place the circles
arbitrarily as long
as their separation, $d$, remains the same. This is not, however, the case
with more than two sets.

With three or more sets, the circles need to be arranged so that
they interfere minimally with one another. In some cases, the set
configuration allows this to be accomplished flawlessly,
but often, compromises must me made. As is often the case
in this context, this turns out to be another optimization problem.
It can be
tackled in many ways; \pkg{eulerr}'s approach is based on a method
developed by~\citet{Frederickson_2015}, which the author describes as
constrained multi-dimensional scaling.

The algorithm tries to position the circles so that the separation
between each pair of circles matches the separation required from~\eqref{eq:dopt}.
If the two sets are disjoint, however, the algorithm is indifferent to the relative
locations of those circles as long as they do not intersect. The equivalent
applies to subset sets: as long as the circle representing the smaller set
remains within the larger circle, their locations are free to vary.
In all other cases, the loss function~\eqref{eq:initLoss} is the residual sums
of squares of the optimal separation of circles, $d$, that we found
in~\eqref{eq:circleOverlap}, and the actual distance in the layout we are
currently exploring.

\begin{fullwidth}
\begin{equation}
  \mathcal{L}(h,k) = \displaystyle\smashoperator[l]{\sum_{1\leq i<j\leq N}}
  \begin{cases}
    0 & F_i \cap F_j = \emptyset \text{ and } O_{ij} = 0\\
    0 & (F_i \subseteq F_j \text{ or } F_i \supseteq F_j) \text{ and } O_{ij}=0\\
   \left(\left(h_i-h_j\right)^2+\left(k_i-k_j\right)^2-d_{ij}^2\right)^2  & \text{otherwise} \\
  \end{cases}. \label{eq:initLoss}
\end{equation}
\end{fullwidth}
The analytical gradient~\eqref{eq:initGrad} is retrieved as usual by
taking the derivative of the loss function,
\begin{fullwidth}
\begin{equation}
  \vec{\nabla} f(h_i) = \sum_{j=1}^N
  \begin{cases}
    \vec{0} & F_i \cap F_j = \emptyset \text{ and } O_{ij} = 0\\
    \vec{0} & (F_i \subseteq F_j \text{ or } F_i \supseteq F_j) \text{ and } O_{ij}=0\\
    4\left(h_i-h_j\right)\left(\left(h_i-h_j\right)^2+\left(k_i-k_j\right)^2-d_{ij}^2\right) & \text{otherwise}, \\
  \end{cases} \label{eq:initGrad}
\end{equation}
\end{fullwidth}
where $\vec{\nabla} f(k_i)$ is found as in~\eqref{eq:initGrad} with $h_i$
swapped for $k_i$ (and vice versa).

The Hessian~\eqref{eq:initHess} for our loss function is given next. However,
because the current release of R suffers from a bug\footnote{The current development version
of R features a fix for this bug; \pkg{eulerr} will be updated to
use \eqref{eq:initHess} as soon as it is introduced in a stable version of R.} causing the
analytical Hessian to be updated improperly, the current
release of \pkg{eulerr} instead relies on the numerical approximation of
the Hessian offered by the optimizer.
\begin{fullwidth}
\begin{equation}
  \bm{H}(h,k) = \smashoperator[l]{\sum_{1\leq i<j\leq N}}
  \begin{bsmallmatrix}
    4\left(\left(h_i-h_j\right)^2+\left(k_i-k_j\right)^2-d_{ij}^2\right)+8\left(h_i-h_j\right)^2 &
    \cdots &
    8\left(h_i-h_j\right)\left(k_i-k_j\right)\\
    \svdots & \sddots & \svdots \\
    8\left(k_i-k_j\right)\left(h_i-h_j\right) &
    \cdots &
    4\left(\left(h_i-h_j\right)^2+\left(k_i-k_j\right)^2-d_{ij}^2\right)+8\left(k_i-k_j\right)^2
  \end{bsmallmatrix}.
\label{eq:initHess}
\end{equation}
\end{fullwidth}
Note that the constraints given in~\eqref{eq:initLoss} and \eqref{eq:initGrad}
still apply to each element of~\eqref{eq:initHess} and have been omitted for
convenience only.

We optimize~\eqref{eq:initLoss} using the nonlinear optimizer \code{nlm()} from
the R core package \pkg{stats}. The underlying code for \code{nlm()} was
written by \citet{Schnabel_1985}. It was ported to R by Saikat DebRoy and the R
Core team~\citep{RCT_2017} from a previous \proglang{FORTRAN} to \proglang{C}
translation by Richard H. Jones. \code{nlm()} consists of a system
of Newton-type algorithms and performs well for
difficult problems~\cite{Nash_2014}.

The initial layout outlined above will sometimes turn up perfect diagrams, but
only reliably so when the diagram is completely determined by its pairwise
intersections. More pertinently, we have not yet considered the higher-order
intersections in our algorithm and neither have we approached the problem
of using ellipses---as we set out to do.

\section{Final layout}
\label{sec:final-config}

We now need to account for all the sets' intersections and, consequently,
all the overlaps in the diagram. The goal is to map each area uniquely
to a subset of the data from the input and for this purpose we will
use the sets' intersections and the relative complements of these
intersections, for which we will use the shorthand $\omega$. We introduced this
form in \cref{sec:input}, but now define it rigorously in \cref{def:omega}.

\begin{mydef}
\label{def:omega}
For a family of \emph{N} sets, $F = F_1, F_2, \dots, F_N$, and their $n=2^N-1$
intersections, we define $\omega$ as the intersections of these sets and their
relative complements, such that
\begin{align*}
  \omega_{1} & = F_1 \setminus \bigcap_{i=2}^N F_i  \\
  \omega_{2} & = \bigcap_{i=1}^2 F_i \setminus \bigcap_{i=3}^{N} F_i\\
             & \vdotswithin{=} \\
    \omega_n & = \bigcap_{i=1}^{N}F_i
\end{align*}
with
\[
  \sum_{i = 1}^n \omega_i =  \bigcup_{j=1}^N F_i.
\]
Analogously to $\omega$, we also introduce the $\&$-operator, such that
\[
  F_i \& F_j = (F_i \cap F_j)\setminus (F_i \cap F_j)^\textsf{c}.
\]
\end{mydef}

The fitted diagram's area-equivalents for $\omega$ will be defined as $A$, so that
an exact diagram requires that $\omega_i = A_i$ for $i=1,2,\dots,2^N-1$,
where $N$ is the number of sets in the input.

In \cref{sec:initConfig}, we restricted
ourselves to circles but now extend ourselves also to ellipses. From now on,
we abandon the practice of treating circles separately---they are only
a special case of ellipses, and, hence, everything that applies to an ellipse does so
equally for a circle.

\subsection{Intersecting ellipses}
\label{sec:intersecting-ellipses}

We now need the ellipses'
points of intersections. \pkg{eulerr}'s approach to this is outlined in
\citet{Richter-Gebert_2011} and based in \emph{projective}, as opposed to
\emph{Euclidean}, geometry.

To collect all the intersection points, we naturally need only to consider two
ellipses at a time. The canonical form of an ellipse is given by
\begin{equation*}
\frac{\left[ (x-h)\cos{\phi}+(y-k)\sin{\phi} \right]^2}{a^2}+
  \frac{\left[(x-h) \sin{\phi}-(y-k) \cos{\phi}\right]^2}{b^2} = 1,
\end{equation*}
where $\phi$ is the counter-clockwise angle from the positive x-axis to the
semi-major axis $a$, $b$ is the semi-minor axis, and $h, k$ are the x- and
y-coordinates, respectively, of ellipse's center~(\cref{fig:ellipse}).

<<ellipse, fig.env="marginfigure", fig.cap="A rotated ellipse with semimajor axis $a$, semiminor axis $b$, rotation $\\phi$, and center $h,k$.">>=
h <- 0.7
k <- 0.5
a <- 1
b <- 0.6
phi <- pi/5
n <- 200
ellipse <- RConics::ellipse(c(a, b), c(h, k), phi,  n = n)

xyplot(1 ~ 1, xlim = c(-0.3, 1.7), ylim = c(-0.3, 1.3), asp = "iso",
       xlab = NULL, ylab = NULL,
       par.settings = list(axis.line = list(col = "transparent")),
       scales = list(draw = FALSE),
       panel = function(x, y, ...) {
         # grid
         panel.refline(h = 0)
         panel.refline(v = 0)

         # rotation
         arc <- ellipse_arc(c(0.2, 0.2), c(h, k), rng = c(0, phi))
         theta <- ellipse_arc(c(a, b), c(h, k), phi, n = 1)

         panel.lines(x = c(h, ellipse[1, 1]), y = c(k, ellipse[1, 2]), lty = 2,
                     col = "grey50")
         panel.lines(x = c(h, theta[1, 1]), y = c(k, theta[1, 2]), lty = 2,
                     col = "grey50")
         panel.lines(arc, col = "grey50")
         panel.text(x = arc[nrow(arc)/2, 1] + 0.06, y = arc[nrow(arc)/2, 2] + 0.02,
                    labels = expression(italic(phi)))

         # semiaxes
         pBrackets::grid.brackets(h, k, ellipse[n/2, 1], ellipse[n/2, 2],
                                  col = "grey50")
         pBrackets::grid.brackets(h, k, ellipse[n/4, 1], ellipse[n/4, 2],
                                  h = 0.04, col = "grey50")

         panel.text(0.37, 0.12, labels = "a", font = 3)
         panel.text(0.43, 0.68, labels = "b", font = 3)

         # center
         panel.text(x = h, y = k + 0.025, pos = 3, label = c("h,k"), font = 3)
         panel.points(x = h, y = k, pch = 21, fill = "white", col = "black")

         # ellipse
         panel.lines(ellipse, col = 1)
       })
@

However, because an ellipse
is a conic\footnote{The circle, parabola, and hyperbola are the other types of conics.} it
can be represented in quadric form,
\begin{equation*}
Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
\label{eq:quadric}
\end{equation*}
that in turn can be represented as a matrix,
\begin{equation*}
\begin{bmatrix}
A   & B/2 & D/2 \\
B/2 & C   & E/2 \\
D/2 & E/2 & F
\end{bmatrix},
\end{equation*}
which is the form we need to intersect our ellipses. We now proceed to
\begin{enumerate}
\item form three degenerate conics from a linear combination of the
  two ellipses we wish to intersect,
\item split one of these degenerate conics into two lines, and
\item intersect one of the ellipses with these lines, yielding 0 to 4
  intersection points points~(\cref{fig:intersection}).
\end{enumerate}

\begin{figure*}[thbp]
\begin{subfigure}[t]{0.32\linewidth}
<<intersect-1>>=
C1 <- ellipseToConicMatrix(c(8, 2), c(0, 0), -pi/3)
C2 <- ellipseToConicMatrix(c(5, 2), c(1, -2), pi/5)

ellipses <- data.frame(rbind(ellipse(c(8, 2), c(0, 0), -pi/3),
                             ellipse(c(5, 2), c(1, -2), pi/5)))
colnames(ellipses) <- c("x", "y")
ellipses$fac <- rep(c("A", "B"), each = 201)
pp <- intersectConicConic(C1, C2)
x <- pp[1, ]
y <- pp[2, ]
a <- double(6)
b <- double(6)
k <- 1
for (i in 1:3) {
  for (j in (i + 1):4) {
    b[k] <- (y[j] - y[i])/(x[j] - x[i])
    a[k] <- b[k]*(-x[i]) + y[i]
    k <- k + 1
  }
}
a <- a[c(1, 6, 3, 4, 2, 5)]
b <- b[c(1, 6, 3, 4, 2, 5)]

xyplot(
  y ~ x, data = ellipses, type = "l", groups = fac, asp = 1,
  xlab = NULL, ylab = NULL,
  scales = list(draw = FALSE, axes = FALSE),
  par.settings = list(axis.line = list(col = "transparent")),
  panel = function(x, y, ...) {
    panel.abline(a = a[3], b = b[3], col = "#FEFEFE")
    panel.abline(a = a[4], b = b[4], col = "#FEFEFE")
    panel.xyplot(x, y, ..., col = 1)
})
@
\subcaption{Our objective is two intersect these two ellipses.}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.32\linewidth}
<<intersect-2>>=
xyplot(y ~ x, data = ellipses, type = "l", groups = fac, asp = 1,
  xlab = NULL, ylab = NULL,
  scales = list(draw = FALSE, axes = FALSE),
  par.settings = list(axis.line = list(col = "transparent")),
  panel = function(x, y, ...) {
    panel.xyplot(x, y, ..., col = 1)
    for (i in 1:length(a)) {
      panel.abline(a = a[i], b = b[i], col = pal[ceiling(i/2) + 1])
    }
})
@
\subcaption{Three degenerate conics (\textcolor{orange}{orange},
\textcolor{Teal}{teal}, and \textcolor{SteelBlue4}{blue}) are formed
from the linear combination of our two ellipses.}
\end{subfigure}
\hfill
\begin{subfigure}[t]{0.32\linewidth}
<<intersect-3>>=
xyplot(
  y ~ x, data = ellipses, type = "l", groups = fac, asp = 1,
  xlab = NULL, ylab = NULL,
  scales = list(draw = FALSE, axes = FALSE),
  par.settings = list(axis.line = list(col = "transparent")),
  panel = function(x, y, ...) {
    panel.xyplot(x, y, ..., col = c("transparent", "black"))
    panel.abline(a = a[3], b = b[3], col = pal[3])
    panel.abline(a = a[4], b = b[4], col = pal[3])
    panel.points(t(pp[1:2, ]), col = 1, pch = 21, fill = "white")
})
@
\subcaption{One of the degenerate conics is split it into two lines and
intersected with one of the ellipses to yield four intersection
points.}
\end{subfigure}
\caption{The process used to intersect two ellipses, here
yielding four points. This figure was inspired by an example from
\citet{Richter-Gebert_2011}.}
\label{fig:intersection}
\end{figure*}

\subsection{Overlap areas}
\label{sec:overlap-areas}

Using the intersection points of a set of ellipses that we retrieved in
\cref{sec:intersecting-ellipses}, we can now find the
overlap of these ellipses. We are only interested in the points that are \emph{contained
within all of these ellipses}, which together form a geometric shape consisting of a
convex polygon, the sides of which are made up of straight lines between
consecutive points, and a set of elliptical arcs---one for each
pair of points~(\cref{fig:polyarea}).

\begin{figure}[hbtp]
\caption{The overlap area between three ellipses is the sum of a convex
polygon (in \textcolor{Grey}{grey}) and 2--3 ellipse segments (in
\textcolor{SteelBlue4}{blue}).\label{fig:polyarea}}
<<polyarea, fig.height = 4.2, fig.width = 4.2>>=
x <- c(0, -0.3, 0.2)
y <- c(0, 0.1, 0.3)
ra <- a <- c(0.3, 0.5, 0.4)
rb <- b <- c(0.3, 0.3, 0.6)
phi <- c(-pi/6, 2, -2)

ee <- data.frame(x, y, ra, rb, phi)

tx <- atan2(-b*tan(phi), a)
ty <- atan2(b*tan(pi/2L - phi), a)

xlim <- range(x + a*cos(tx)*cos(phi) - b*sin(tx)*sin(phi),
              x + a*cos(tx + pi)*cos(phi) - b*sin(tx + pi)*sin(phi))
ylim <- range(y + b*sin(ty)*cos(phi) + a*cos(ty)*sin(phi),
              y + b*sin(ty + pi)*cos(phi) + a*cos(ty + pi)*sin(phi))

pp <- matrix(NA, ncol = 4, nrow = 0)

for (i in 1:2) {
  for (j in (i + 1):3) {
    e1 <- ellipseToConicMatrix(c(ra[i], rb[i]), c(x[i], y[i]), phi[i])
    e2 <- ellipseToConicMatrix(c(ra[j], rb[j]), c(x[j], y[j]), phi[j])
    pp <- rbind(pp, cbind(t(intersectConicConic(e1, e2)[1:2, ]), i, j))
  }
}

sel <- logical(nrow(pp))
for (k in 1:nrow(pp)) {
  in_which <- ((pp[k, 1] - x)*cos(phi) + (pp[k, 2] - y)*sin(phi))^2/ra^2 +
    ((pp[k, 1] - x)*sin(phi) - (pp[k, 2] - y)*cos(phi))^2/rb^2 <= 1 + 0.1
  sel[k] <- all(in_which)
}

pp <- pp[sel, ]
mid <- cbind(mean(pp[, 1]), mean(pp[, 2]))
seglines <- matrix(NA, ncol = 2, nrow = 0)

ang <- atan2(pp[, 1] - mid[1], pp[, 2] - mid[2])
ord <- order(ang)

pp <- pp[ord, ]

j <- nrow(pp)
for (i in 1:nrow(pp)) {
  k <- intersect(pp[i, 3:4], pp[j, 3:4])
  start <- atan2(pp[j, 2] - y[k], pp[j, 1] - x[k])
  stop <- atan2(pp[i, 2] - y[k], pp[i, 1] - x[k])
  arc <- ellipse_arc(c(a[k], b[k]), c(x[k], y[k]), theta = phi[k],
                     rng = c(start, stop))
  seglines <- rbind(seglines, arc)
  j <- i
}

els <- eulerr:::ellipse(x, y, ra, rb, phi)
x0 <- c(lapply(els, "[[", "x"), recursive = TRUE)
y0 <- c(lapply(els, "[[", "y"), recursive = TRUE)

xyplot(
  y ~ x, data = ee, asp = "iso",
  xlim = extendrange(xlim, f = 0.01), ylim = extendrange(ylim, f = 0.01),
  scales = list(draw = FALSE), xlab = NULL, ylab = NULL,
  par.settings = list(axis.line = list(col = "transparent")),
  panel = function(x, y, ...) {
    grid.polygon(x0, y0,  id.lengths = rep(200, 3),
                 default.units = "native",
                 gp = grid::gpar(fill = "transparent"))
    panel.polygon(seglines, col = "slategray2")
    panel.polygon(pp[, 1:2], col = "grey90")
    panel.points(mid, pch = 4, col = 1, cex = 1.5)
    panel.points(pp[, 1:2, drop = FALSE], col = 1, pch = 21, fill = "white",
                 cex = 1.5)
  }
)
@
\end{figure}

We continue by ordering the points around their centroid. It is then trivial to
find the area of the
polygon section since it is always convex~\citep{Finley_2006}. Now, because each
elliptical segment is formed from the arcs that connect successive
points, we can establish the segments' areas algorithmically~\citep{Eberly_2016}.
For each ellipse and its related pair of points (located at angles $\theta_0$ and
$\theta_1$ from the semimajor axis), we proceed to find its area by
\begin{enumerate}
\item centering the ellipse at $(0, 0)$,
\item normalizing its rotation, which is not needed to compute the area,
\item integrating the ellipse over [$0,\theta_0$] and [$0,\theta_1$], producing
  elliptical sectors $F(\theta_0)$ and $F(\theta_1)$,
\item subtracting the smaller ($F(\theta_0$)) of these sectors from the larger ($F(\theta_0$), and
\item subtracting the triangle section to finally find the segment area,
  \begin{equation*}
  F(\theta_1) - F(\theta_0) - \frac{1}{2}\left|x_1y_0 - x_0y_1\right|,
  \label{eq:segmentArea}
  \end{equation*}
  \[
  \text{where } F(\theta) = \frac{a}{b}\left[ \theta - \arctan{\left(\frac{(b - a)\sin{2\theta}}{b + a +(b - a )\cos{2\theta}} \right)}\right].
  \]
\end{enumerate}
This procedure is illustrated in~\cref{fig:ellipsesegment}. Note that there are
situations where this algorithm is altered, such that when the sector angle
ranges beyond $\pi$---we refer the interested reader to \citet{Eberly_2016}.

\begin{marginfigure}
<<ellipsesegment>>=
ellipse <- ellipse(c(1, 0.6), c(0, 0), 0)
i0 <- 15
i1 <- 45

tri <- rbind(cbind(0, 0),
             ellipse[i0, , drop = FALSE],
             ellipse[i1, , drop = FALSE])

xyplot(1~1, xlim = c(-1.2, 1.2), ylim = c(-0.8, 0.8), asp = "iso",
       xlab = NULL, ylab = NULL,
       par.settings = list(axis.line = list(col = "transparent")),
       scales = list(draw = FALSE),
       panel = function(x, y, ...) {
         panel.grid(x = 0, y = 0, h = 1,v = 1)
         panel.lines(ellipse, col = 1)
         panel.polygon(ellipse[i0:i1, ], col = "steelblue1", alpha = 0.5)
         panel.polygon(tri, col = "grey95")
         panel.points(ellipse[c(i0, i1), ], col = 1, pch = 21, fill = "white")
         panel.text(ellipse[c(i0, i1), ],
                    adj = c(-0.5, -0.5),
                    labels = c(expression(italic(theta[0])),
                               expression(italic(theta[1]))))
                    pBrackets::grid.brackets(1, 0, 0, 0, h = 0.05,
                                  type = 1)
         pBrackets::grid.brackets(0, 0, 0, 0.6, h = 0.04)
         panel.text(0.5, unit(-0.05, "npc"), labels = "a", font = 3, pos = 1)
         panel.text(-0.05, unit(0.3, "native"), labels = "b", font = 3, pos = 2)
       })
@
\caption{The elliptical segment in \textcolor{SteelBlue4}{blue} is found by
first subtracting the elliptical sector from $(a, 0)$ to $\theta_0$ from the one
from $(a, 0)$ to $\theta_1$ and then subtracting the triangle part
(in \textcolor{Grey}{grey}).}
\label{fig:ellipsesegment}
\end{marginfigure}

Finally, the area of the overlap is then obtained by adding the area of the
polygon and all the elliptical arcs together.

Note that this does not yet give us
the areas that we require, namely $A$: the area-equivalents to the set intersections
and relative complements from~\cref{def:omega}. For this, we must
decompose the overlap areas so that each area maps uniquely
to a subspace of the set configuration. This, however, is simply a matter of transversing down the
hierarchy of overlaps and subtracting the higher-order overlaps from the
lower-order ones. For a three-set relationship of sets $A$, $B$, and $C$,
for instance, this means subtracting the $A\cap B \cap C$ overlap from
the $A \cap B$ one to retrieve the equivalent of $(A \cap B) \setminus C$.

The exact algorithm may in rare instances\footnote{1 out of approximately 7000 in
our simulations.}, break down, the culprit being numerical precision
issues that occur when ellipses are tangent or
completely overlap. In these cases, the algorithm will approximate
the area of the involved overlap by
\begin{enumerate}
\item spreading points across the ellipses using Vogel's
  method~(see \cref{sec:labeling} for a brief introduction),
\item identifying the points that are inside the intersection via the inequality
  \begin{multline*}
  \frac{\left[(x-h)\cos{\phi}+(y-k)\sin{\phi} \right]^2}{a^2} + \\
    \frac{\left[(x-h) \sin{\phi}-(y-k)\cos{\phi}\right]^2}{b^2} < 1,
  \end{multline*}
  where $x$ and $y$ are the coordinates of the sampled points, and finally
\item approximating the area by multiplying the proportion of points inside the
  overlap with the area of the ellipse.
\end{enumerate}

With this in place, we are now able to compute the areas of all intersections
and their relative complements, $\omega$, up to numerical precision.

\subsection{Final optimization}
\label{sec:final-optimization}

We feed the initial
layout computed in~\cref{sec:initConfig} to the optimizer---once again
we employ \code{nlm()} from \pkg{stats} but now also provide the option to
use ellipses rather than circles, allowing the ``circles''
to rotate and the relation between the semiaxes to vary, altogether
rendering five parameters to optimize per set and ellipse (or three if we
restrict ourselves to circles). For each iteration of
the optimizer, the areas of all intersections are analyzed and a measure of loss
returned. The loss we use is the same as that in \pkg{venneuler}~\cite{Wilkinson_2012},
namely \emph{stress},
\begin{equation}
 \frac{\sum_{i=1}^n (A_i - \beta \omega_i)^2}{\sum_{i=1}^n A_i},
\label{eq:stress}
\end{equation}
where \[\beta = \frac{\sum_{i=1}^n A_i\omega_i}{\sum_{i=1}^n \omega_i^2}.\]
This is equivalent to linear regression through the origin, where $\beta$ is the
slope of the regression line~(\cref{fig:origin-regress}).

<<origin-regress, fig.env="marginfigure", fig.cap="Optimizing via stress is analogous to least-squares linear regression through the origin. $\\omega$ is the set of unique quantities in the input~(\\cref{def:omega}) and $A$ the respective areas in the diagram.",fig.width=2.4,fig.height=2.4>>=
fit <- euler(c(X = 3, Y = 1, Z = 2,
               "X&Y" = 2, "X&Z" = 1.2, "Y&Z" = 1.5,
               "X&Y&Z" = 1))
d <- data.frame(A = fitted(fit), omega = fit$original.values)
xyplot(A ~ omega, d, groups = as.factor(1:7), asp = "iso",
       xlim = c(-0.1, 3.9), ylim = c(-0.1, 3.9),
       xlab = expression(omega),
       panel = function(x, y, ...) {
         panel.abline(a = 0, b = sum(d$A*d$omega)/(sum(d$omega^2)),
                       rotate = TRUE, at = 0.8, pos = 1,
                       col = "grey60")
         panel.xyplot(x, y, col = trellis.par.get("superpose.symbol")$col, ...)
         panel.text(x, y, labels = names(resid(fit)),
                    col = trellis.par.get("superpose.symbol")$col,
                    pos = c(3, 3, 3, 2, 3, 4, 2))
       })
@

\subsection{Last-ditch optimization}
\label{sec:last-ditch}

If the fitted diagram is still inexact after the procedure in
\cref{sec:final-optimization}, we offer a final step in which we
pass the parameters on to a last-ditch optimizer. The weapon of
choice\footnote{We conducted thorough benchmarking, that we opt not to report here, to
decide upon an algorithm for this step.} is \emph{Generalized simulated annealing}
from the R package \pkg{GenSA}~\citep{Xiang_2013}. The solutions offered by
\pkg{GenSA} often avoid local minima.

By default, this last-ditch step is activated only when we have a three-set diagram
with ellipses and a diagError~\eqref{eq:diagError} above
0.001\footnote{The choice of if and when this last-ditch optimizer is
activated is left to the user via simple commands to the main function
of the package.} The
reason being that the method is considerably more computationally intensive.

\section{Goodness of fit}
\label{sec:gof}

Every Euler diagram must be investigated for its adequacy in representing the
input. Exact Euler diagrams are not always possible
When \pkg{eulerr} cannot find a perfect solution, it offers an approximate one
instead, the adequacy of which has to be measured in a standardized way. For
this purpose we adopt two measures: \emph{stress}~\citep{Wilkinson_2012}, which
is also the loss metric we use in our final optimization step and is
used in \pkg{venneuler}, as well as
\emph{diagError}~\citep{Micallef_2014a}, which is used by \pkg{eulerAPE}.

The stress metric is not easily grasped but
can be transformed into a rough analogue of the correlation coefficient via
$r = \sqrt{1-\text{stress}^2}$.

diagError, meanwhile, is given by
\begin{equation}
{\max_{i = 1, 2, \dots, n}}\left|
  \frac{\omega_i}{\sum_{i=1}^n \omega_i} - \frac{A_i}{\sum_{i=1}^n A_i} \right|,
\label{eq:diagError}
\end{equation}
which is the maximum \emph{absolute} difference of the proportion of any
$\omega$ to the respective unique area of the diagram.

\section{Plotting}
\label{sec:plotting}

Once we have ascertained that our Euler diagram fits well, we can turn to
visualizing the solution. For this purpose, \pkg{eulerr} relies on the
\pkg{grid} graphics system~\citep{RCT_2017} and offers intuitive and
granular control over the output.

Plotting the ellipses is straightforward using the parametrization of a rotated
ellipse,
%
\begin{equation*}
\begin{bmatrix}
  x \\ y
\end{bmatrix} =
\begin{bmatrix}
  h + a \cos{\theta} \\
  k + b \sin{\theta}
\end{bmatrix},\quad \text{where } \theta \in [0, 2\pi],\, a,b>0.
\end{equation*}
%
Most users will also prefer to label the ellipses and their intersections
with text and this, however, is considerably more involved.

\subsection{Labeling}
\label{sec:labeling}

Labeling the ellipses is complicated since the shapes of the intersections
often are irregular, lacking well-defined centers; we know of no analytical
solution to this problem. As usual, however, the next-best option turns out to
be a numerical one. First, we locate a point that is inside the required region
by spreading points across the discs involved in the set intersection. To
distribute the points, we use a modification of
\emph{Vogel's method}~\citep{Arthur_2015,Vogel_1979} adapted to ellipses.
Vogel's method spreads points across a disc using
\begin{equation}
p_k =
\begin{bmatrix}
  \rho_k \\
  \theta_k
\end{bmatrix} =
\begin{bmatrix}
  r \sqrt{\frac{k}{n}}\\
  \pi (3 - \sqrt{5})(k - 1)
\end{bmatrix}\quad\text{for } k = 1, 2,\dots, n.
\label{eq:vogel}
\end{equation}
In our modification, we scale, rotate, and translate the points formed
in~\eqref{eq:vogel} to match the candidate ellipse. We rely, as before, on
projective geometry to carry out the transformations in one go:
\[
p' =
\begin{bmatrix}
  x' \\
  y' \\
  1
\end{bmatrix} =
\begin{bmatrix}
  1 & 0 & h \\
  0 & 1 & k \\
  0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
  \cos{\phi}  & \sin{\phi} & 0 \\
  -\sin{\phi} & \cos{\phi} & 0\\
  0           & 0          & 1
\end{bmatrix}
\begin{bmatrix}
  a & 0 & 0 \\
  0 & b & 0 \\
  0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
  \hat{x} \\
  \hat{y} \\
  1
\end{bmatrix},
\]
where $h,k$ translates, $\phi$ rotates, and
$a,b$ stretches the ellipse.

After we spread our points throughout the ellipse and find a point, $p'_i$, that
is contained in our desired intersection, we proceed to optimize its position
numerically. The position we are looking for is that which maximizes the
distance to the closest ellipse in our diagram to provide as much margin as
possible for the label. This is a maximization problem with a loss function equal to
\begin{equation}
\mathcal{L}(x,y) = \min_{i=1,2,\dots,N} f(x,y,h_i,k_i,a_i,b_i,\phi_i)
\label{eq:lossDist}
\end{equation}
where $f$ is the function that determines the distance from a point ($x,y$) to
the ellipse defined by $h,k,a,b$ and $\phi$.

Similarly to fitting Euler diagrams in the general case, there appears to be no
analytical solution to computing the distance from a point to an ellipse. The
numerical solution we use has been described by~\citet{Eberly_2016a} and
involves solving the roots to a quartic polynomial via a robust
bisection optimizer.

To optimize the location of the label, we employ a version of the
\emph{Nelder--Mead method}~\citep{Nelder_1965}, which has been translated from
a \proglang{Matlab} code by \citet{Kelley_1999} and adapted for \pkg{eulerr} to
ensure that it converges quickly and that the simplex remains within the
intersection boundaries (since we want the local maximum). The method is
visualized in~\cref{fig:vogel}.
\begin{marginfigure}
<<vogel, fig.show = "hold">>=
h <- c(0.2, 0.6)
k <- c(0.4, 0.4)
a <- c(0.4, .5)
b <- c(0.3, 0.6)
phi <- c(2.5, -1)
n <- 250
seqn <- seq(1, n, 1)
theta <- seqn*pi*(3 - sqrt(5))
rad <- sqrt(seqn/n)
x <- rad*cos(theta)
y <- rad*sin(theta)
p <- rbind(x, y, 1)

# Scale, rotate and translate to match the ellipse
p0 <- translation(c(h[1], k[1])) %*%
   rotation(-phi[1]) %*%  scaling(c(a[1], b[1])) %*% p

inside <- ((p0[1, ] - h[2])*cos(phi[2]) + (p0[2, ] - k[2])*sin(phi[2]))^2/a[2]^2 +
  ((p0[1, ] - h[2])*sin(phi[2]) - (p0[2, ] - k[2])*cos(phi[2]))^2/b[2]^2 <= 1

# Set up limits
tx <- atan2(-b*tan(phi), a)
ty <- atan2(b*tan(pi/2L - phi), a)

xlim <- range(h + a*cos(tx)*cos(phi) - b*sin(tx)*sin(phi),
              h + a*cos(tx + pi)*cos(phi) - b*sin(tx + pi)*sin(phi))
ylim <- range(k + b*sin(ty)*cos(phi) + a*cos(ty)*sin(phi),
              k + b*sin(ty + pi)*cos(phi) + a*cos(ty + pi)*sin(phi))

xyplot(
  x ~ y, aspect = "iso", pch = 20, xlab = "", ylab = "",
  col = 1,
  xlim = extendrange(xlim, f = 0.01),
  ylim = extendrange(ylim, f = 0.01),
  par.settings = list(axis.line = list(col = "transparent")),
  scales = list(draw = FALSE),
  panel = function(x, y, ...) {
    panel.points(t(p0[1:2, !inside]), pch = 19, col = "grey80", cex = 0.5)
    panel.points(t(p0[1:2, inside]), pch = 19, col = "steelblue2", cex = 0.5)
    panel.points(t(p0[1:2, !inside][, 14]), pch = 19, col = 1, cex = 0.5)
    panel.polygon(ellipse(c(a[1], b[1]), c(h[1], k[1]), phi[1]))
    panel.polygon(ellipse(c(a[2], b[2]), c(h[2], k[2]), phi[2]))
  })

centers <- eulerr:::locate_centers(h, k, a, b, phi, 1:2)

xyplot(
  x ~ y, aspect = "iso", pch = 19, xlab = "", ylab = "",
  col = 1,
  xlim = extendrange(xlim, f = 0.01),
  ylim = extendrange(ylim, f = 0.01),
  par.settings = list(axis.line = list(col = "transparent")),
  scales = list(draw = FALSE),
  panel = function(x, y, ...) {
    panel.arrows(p0[1, !inside][14], p0[2, !inside][14],
                 centers[1, 1], centers[2, 1],
                 length = 0.05, col = "grey70")
    panel.points(t(p0[1:2, !inside][, 14]), pch = 19, col = 1, cex = 0.5)
    panel.polygon(ellipse(c(a[1], b[1]), c(h[1], k[1]), phi[1]))
    panel.polygon(ellipse(c(a[2], b[2]), c(h[2], k[2]), phi[2]))
  })

xyplot(
  x ~ y, aspect = "iso", pch = 19, xlab = "", ylab = "",
  col = 1,
  xlim = extendrange(xlim, f = 0.01),
  ylim = extendrange(ylim, f = 0.01),
  par.settings = list(axis.line = list(col = "transparent")),
  scales = list(draw = FALSE),
  panel = function(x, y, ...) {
    panel.text(t(centers[1:2, 1]), labels = "42")
    panel.polygon(ellipse(c(a[1], b[1]), c(h[1], k[1]), phi[1]))
    panel.polygon(ellipse(c(a[2], b[2]), c(h[2], k[2]), phi[2]))
  })

#grid.arrange(p1, p2, p3, ncol = 1)
@
\caption{The method eulerr uses to locate an optimal position for a label in
three steps from top to bottom: first, we spread sample points on one of the
ellipses and pick one inside the intersection of interest, then we begin moving
it numerically, and finally place our label.}
\label{fig:vogel}
\end{marginfigure}

\subsection{Aesthetics}
\label{sec:aesthetics}

Euler diagrams display both quantitative and qualitative data. The quantitative
aspect is the quantities or sizes of the sets depicted in the diagram and is
visualized by the relative sizes, and possibly the labels, of the areas of the
shapes---this is the main focus of this paper. The qualitative aspects,
meanwhile, consist of the mapping of each set to some quality or category, such
as having a certain gene or not. In the diagram, these qualities can be
separated through any of the following aesthetics:
%
\begin{itemize}
\item color,
\item border type,
\item text labelling,
\item transperancy,
\item patterns,
\end{itemize}
%
or a combination of these. The main purpose of these aethetics is to separate
out the different ellipses so that the audience may interpret the diagram with
ease and clarity.

Among these aesthetics, the best choice (from a viewer perspective) appears to
be color~\citep{Blake_2016}, which provides useful information without
extraneous chart junk~\citep{Tufte_2001}. The issue with color, however, is that
it cannot be perceived perfectly by all---8\% of men and 0.4\% of women in
European Caucasian countries, for instance, suffer the most common form,
red--green color deficiency~\cite{Birch_2012}. Moreover, color
is often printed at a premium in scientific publications and adds no information
to a diagram of two shapes.

For these reasons, \pkg{eulerr} defaults to distinguishing ellipses with color
using a color palette generated via the R package
\pkg{qualpalr}~\citep{Larsson_2016}, which automatically generates qualitative
color palettes based on a perceptual model of color vision that optionally
caters to color vision deficiency. This palette has been manually modified
to fullfil our other objectives of avoiding using colors for two
sets. The first eight colors of the pallete are visualized in
\cref{fig:colorexample}.

\begin{figure}[htbp]
\caption{The eight first colors of the default color palette.\label{fig:colorexample}}
<<colorexamle, fig.width = 4.2, fig.height = 3>>=
set.seed(2)
con <- c(A = 1, B = 1, C = 1, D = 1, E = 1, F = 1, G = 1, H = 1,
         "A&B" = 0.2, "B&C" = 0.2, "C&D" = 0.2, "D&E" = 0.2, "E&F" = 0.2,
         "F&G" = 0.2, "G&H" = 0.2)
plot(euler(con), labels = as.character(1:8))
@
\end{figure}

\section{Normalizing dispered layouts}
\label{sec:layout}

If there are disjoint clusters of ellipses, the optimizer will often
spread these out more than is necessary, wasting space in our diagram. To
tackle this, we
use a SKYLINE-BL rectangle packing algorithm~\citep{Jylaenki_2010}
designed specifically for \pkg{eulerr}. In it, we surround each ellipse cluster
with a bounding box, pack these boxes into a bin of appropriate size and
aspect ratio, and adjust the coordinates of the ellipses in the clusters to
compact our diagram. As a bonus, this increases the chance of having
similar layouts for different function calls.

\clearpage
\begin{fullwidth}
\bibliography{eulerr}
\bibliographystyle{unsrtnat}
\end{fullwidth}

\end{document}
