\documentclass[10pt]{elsart}

\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{color}

\newtheorem{theorem}{Theorem}
\newtheorem{remark}{Remark}
\newtheorem{lemma}{Lemma}
\newtheorem{corollary}{Corollary}



\begin{document}
\begin{frontmatter}
\title{Cubature rules with positive weights and internal  nodes on spline curvilinear domains}



\thanks[project]{Work partially supported by the DOR funds and the biennial project BIRD163015 of the University of Padova, and by the GNCS-INdAM.}

 \author[padua]{A. Sommariva},
 \ead{alvise@math.unipd.it}
 \author[padua]{M. Vianello}
 \ead{marcov@math.unipd.it}
 
 \address[padua]{Dept. of Mathematics, University of Padova\\
via Trieste 63, 35121 Padova (Italy)}




\begin{abstract}
In this paper we introduce an algorithm that computes an algebraic cubature rule of degree $n$ with at most $(n+1)(n+2)/2$ internal nodes and positive weights, over spline curvilinear regions $\cal{S}$. The key ingredients are a theorem by Davis on Tchakaloff sets in compact domains and a new {\it{in-domain}} algorithm in such $\cal{S}$. Numerical Matlab codes are provided.
\end{abstract}



\begin{keyword}
Algebraic cubature, splines, curvilinear polygons, moment computation, compression, monnegative-least squares, point containment, winding number.
 \end{keyword}
 
\end{frontmatter}

 {\scriptsize{
 {\bf{2010 AMS subject classification}}: Primary 41A10, 65D32; Secondary: 65D10.
}}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





\section{Introduction}
During the last years, there has been an increasing interest on algebraic cubature rules over polygonal and curvilinear domains. 

In particular, in Virtual Elements methods, there has been requests of formulae on regions that have non standard shapes, but whose boundaries can be typically well approximated by parametric splines (see e.g. \cite{BRV18} and references therein). In these instances, the algebraic degree of precision, term often shortened by the acronym ADE, is usually not too high, but rules of PI-type, i.e. with nodes in the domain and positive weights, are in many senses more appealing.

In the past we wrote several contributions on these topics. In \cite{SV07}, \cite{SV09}, \cite{SSV11} we discussed cubature rules respectively over simple polygons, spline curvilinear and simple domains whose boundary can be described by parametric functions, also providing Matlab routines.

In all these works, the rules were of PI-type in convex domains but this property was not always guaranteed in non convex regions.

To fulfil these requests, we have recently introduced in \cite{BSV19} an algorithm that determines formula of PI-type and prescribed degree of precision, even on polygons with complicated shapes, that can be non convex, not simple and not simply connected domains. In particular, if the cardinality of the rules is an issue, we were able to provide rules of PI-type, ADE equal to $n$ and a number of nodes at most $N_n=(n+1)(n+2)/2$ by means of a numerical implementation of Tchakaloff theorem (see \cite{T57}, \cite{T15}, \cite{SV15}). The relative Matlab codes make use of the new {\tt{polyshape}} environment and can be retrieved at \cite{SH}.

The spline curvilinear counterpart is more delicate and we propose here an algorithm that provides rules of PI-type not only over polygons but also on {\it{simple}} spline curvilinear domains. 

{\color{blue}
Here we propose an algorithm that first, by means of a new in-domain routine for spline curviliear polygons, determines a {\sl{sufficiently dense}} set of bivariate points $\{\xi_i\}$ belonging to the spline curvilinear domain $\cal{S}$ and then it evaluates the Vandermonde-like matrix $V=\{V_n(X)\}_{i,j}=\phi_j(\xi_i)$, w.r.t. a certain polynomial basis $\{\phi_i\}$ of total degree $\delta$ over $\cal{S}$. Next, it computes the moments $\gamma_j=\int_{\cal{S}} \phi_j d{\cal{S}}$, and finally uses the numerical implementation of Tchakaloff theorem, computing a sparse solution of the system $V \cdot w=\gamma$ where the vector of the weights $w$ has at most $(n+1)(n+2)/2$ non null components. The existence of such solution relies on a result on Tchakaloff theorem, proved in {\cite{W}}.

The paper is organized as follows. In Section 2, we introduce our {\it{in-domain}} routine for spline curvilinear domains ${\cal{S}}$. In Section 3 we briefly recall some main results about Tchakaloff theorem and its application In Section 4 we explain our new method for determining cubature rules of PI-type in ${\cal{S}}$, pointing out how to compute the required moments $\gamma_k$ by means of Gauss-Green theorem as well as reporting the fundamental results of {\cite{W}}, describing in detail the proposed algorithm and some approximation properties. Finally in Section 5 we show the numerical results on two non convex domains, regarding either the cubature either the {\it{in-domain}} procedure.
}

\section{In domain algorithm for spline curvilinear regions}
Let ${\cal{S}} \subset {\mathbb{R}}^2$ be a simply connected domain and suppose that its boundary $\partial {\cal{S}}$ 
\begin{enumerate}
\item can be described by parametric equations $x={\tilde{x}}(t)$, $y={\tilde{y}}(t)$, $t \in [a,b]$, ${\tilde{x}},{\tilde{y}} \in C([a,b])$ such that ${\tilde{x}}(a)={\tilde{x}}(b)$ and ${\tilde{y}}(a)={\tilde{y}}(b)$;
\item is not self-intersecting (i.e. ${\cal{S}}$ is {\it{simple}});
\item there is a partition $\{I^{(k)}\}_{k=1,\ldots,M}$ of $[a,b]$, and partitions $\{I_j^{(k)}\}_{j=1,\ldots,m_k}$ of each $I^{(k)}$, such that the restrictions of ${\tilde{x}}$, ${\tilde{y}}$ on each $I^{(k)}$ are splines of degree $\delta_k$, w.r.t. the subintervals $\{I_j^{(k)}\}_{j=1,\ldots,m_k}$.
\end{enumerate}

Typical instances are {\it{approximations}} of more general and piecewise regular, not self-intersecting, simply connected domains $\Omega$ by means of spline curvilinear regions ${\cal{S}}$. To this purpose, given the {\it{vertices}} $V_k \in \partial \Omega$, $k=1,\ldots,M+1$, then $\partial \Omega \approx \partial {\cal S}:=\cup_{k=1}^M V_k \frown V_{k+1}$ where each curvilinear side $V_k \frown V_{k+1}$ is tracked by a spline curve of degree $\delta_k$, interpolating an ordered subsequence of control points $P_{1,k}=V_k,P_{2,k},\ldots,P_{m_k-1,k},P_{m_k,k}=V_{k+1}$ with a suitable parametrization that determines each $I_j^{(k)}$ (and thus each $I^{(k)}$).

In this section we introduce an algorithm that establishes if a point $P \in {\mathbb{R}}^2$ is inside, outside or on the boundary of $\partial {\cal{S}}$. 

\begin{figure}[!htbp]
  \centering
   {\includegraphics[scale=0.25,clip]{monotone_boxes.png}}
 \caption{A domain $\cal S$ of (spline) curvilinear type and its {\it{monotone boxes}}.}
 \label{fig_1}
 \end{figure}

We start computing a domain that consists of union of rectangles and contains $\partial {\cal{S}}$. To this purpose, let $I_j^{(k)}=[t^{(k)}_j, t^{(k)}_{j+1}]$ the partitions introduce above of the parametrization interval $[a,b]$. Furthermore, if ${\tilde{x}}^{\prime}$ changes sign in $(t^{(k)}_j, t^{(k)}_{j+1})$, then define as ${\cal{N}}_j^{(k)}=\{t^{(j,k)}_{i}\}_{i=1,\ldots,l_{j,k}}$ the set of all the points $t^{(j,k)}_{i} \in (t^{(k)}_j, t^{(k)}_{j+1})$ such that ${\tilde{x}}^{\prime}(t^{(j,k)}_{i})=0$ (observe that the restriction to $I^{(j,k)}$ is a polynomial of degree $\delta_k$, hence its derivative exists), otherwise put ${\cal{N}}_j^{(k)}=\emptyset$.   Let $T^{(j,k)}=\{t^{(k)}_j, t^{(k)}_{j+1}\} \cup {\cal{N}}_j^{(k)}$, where we suppose that its elements $T^{(j,k)}_i$ are in increasing order.
Next introduce the rectangles ${\cal{B}}^{(j,k)}_i$, that we will call {\it{monotone boxes}}, 

{\color{blue}
$${\cal{B}}_i^{(j,k)}:=
[\min_{t \in \{T^{(j,k)}_i \cup T^{(j,k)}_{i+1}]\}} {\tilde{x}}(t), \max_{t \in \{T^{(j,k)}_i \cup T^{(j,k)}_{i+1}]\}} {\tilde{x}}(t)] 
\times [\min_{t \in [T^{(j,k)}_i,T^{(j,k)}_{i+1}]} {\tilde{y}}(t), \max_{t \in [T^{(j,k)}_i,T^{(j,k)}_{i+1}]} {\tilde{y}}(t)].$$
}

Observe that by definition if ${\cal{N}}_j^{(k)}=\emptyset$, then there is only the monotone box ${\cal{B}}^{(j,k)}_1$. 

Furthermore, since $\tilde{y}$ restricted to $[T^{(j,k)}_i,T^{(j,k)}_{i+1}]$ is a polynomial of degree $\delta_k$, the evaluation of 
$$\min_{t \in [T^{(j,k)}_i,T^{(j,k)}_{i+1}]} {\tilde{y}}(t), {\mbox{   }} \max_{t \in [T^{(j,k)}_i,T^{(j,k)}_{i+1}]} {\tilde{y}}(t)$$ can be easily determined once the derivative ${\tilde{y}}^{\prime}$ is at hand, by computing its zeros in $[T^{(j,k)}_i,T^{(j,k)}_{i+1}]$ and the evaluation of ${\tilde{y}}$ at $T^{(j,k)}_i$ and $T^{(j,k)}_{i+1}$.

By construction, the restriction of ${\tilde{x}}$ to each monotone box ${\cal{B}}^{(j,k)}_i$  is a monotone function.

Once the set ${\cal{B}}:=\{ {\cal{B}}^{(j,k)}_i \}$ is at hand, we apply the crossing theorem to determine if $P=(P_x,P_y)$ is in the domain $ {\cal{S}}$. We consider a point $Q=(Q_x,Q_y) \not \in {\cal{S}}$ with $Q_y < P_y$, and see how many times the segment $QP$ crosses $\partial \cal{S}$. 

To this purpose we the consider the monotone boxes 
$${\cal{B}(P)}= \{   B = [\alpha_1,\beta_1] \times  [\alpha_2,\beta_2] \in {\cal{B}}: P_x \in [\alpha_1,\beta_1] , P_y \geq \alpha_2 \}$$ 
Take the generic monotone box $B_l= [\alpha^{(l)}_1,\beta^{(l)}_1] \times  [\alpha^{(l)}_2,\beta^{(l)}_2] \in {\cal{B}(P)}$. If $P_y > \beta^{(l)}_2 $ then the point $P$ is {\it{over}} the monotone box $B_l$ and necessarily the segment $QP$ crosses the boundary ${\partial {\cal{S}}}$ once in $B_l$ and {\it{below}} $P$, due to the monotonicy of $\tilde{x}$ in $B_l$. 

Otherwise the point $P \in B_l$. Let $t^*$ be the unique solution of the polynomial equation ${\tilde{x}}(t)=P_x$. Next, 
\begin{itemize}
\item if ${\tilde{y}} (t^*) < P_y$ then the segment $QP$ crosses the boundary ${\partial {\cal{S}}}$ once in the monotone box, {\it{below}} $P$;
\item if ${\tilde{y}} (t^*) > P_y$ then the segment $QP$ does not cross the boundary ${\partial {\cal{S}}}$ once in $B$, {\it{below}} $P$;
\item if ${\tilde{y}} (t^*) = P_y$ then $P$ is on the boundary $\partial {\cal{S}}$.
\end{itemize}
 
Thus counting all these crossings, we are {\it{usually}} able to determine if a point is inside, outside or on the boundary of ${\cal{S}}$. 

Unfortunately $\partial \cal{S}$ can have points $S_k=({\tilde{x}}(t_k),{\tilde{y}}(t_k))$ where $$\lim_{t \rightarrow t^-_k}{\tilde{y}}^{\prime}(t) \lim_{t \rightarrow t^+_k}{\tilde{y}}^{\prime}(t) < 0,$$
that geometrically means that there is locally a vertical turn of boundary from left to right (or conversely from right to left). Even worst, there can be points $P$ where the boundary is locally a vertical segment and it is not straightforward to determine how many times the line $l$ crosses $\partial {\cal{S}}$ {\it{below}} $P$. In these {\it{singular}} cases we have used an algorithm that is based on  the well-known {\it{winding}} theorem to ascertain whether $P$ belongs or not to $\cal{S}$. 
In particular, to determine if $P=(P_x,P_y)$ belongs to $\cal{S}$, we computed by a (shifted) Gauss-Legendre rule of sufficiently high degree of precision the so called {\it{winding number}} ${\mbox{wind}}(P,{\tilde{x}},{\tilde{y}}) \in {\mathbb{Z}}$,
$$
{\mbox{wind}}(P,{\tilde{x}},{\tilde{y}}) := {\frac{1}{b-a}}\int_a^b  \frac {{\tilde{y}}^{\prime} ( {\tilde{x}}(t)-P_x  ) - {\tilde{x}}^{\prime} ( {\tilde{y}}(t)-P_y  )} { ({\tilde{x}}(t)-P_x)^2  +  ({\tilde{y}}(t)-P_y)^2  } dt.
$$
If this quantity is odd then the point belongs to $\cal{S}$ otherwise is not inside such domain. We point out that the Hormann-Agathos algorithm is based on a clever computation of this quantity whenever $\cal{S}$ is a polygon.

\begin{remark} A common istance is to determine what points of a set $X \subset {\mathbb{R}}^2$ belong to $\cal S$. In this case, the construction of the rectangles $\cal B$ is made once for all. 
\end{remark}

\begin{remark} 
{\color{blue}
In our implementation, each monotone box ${\cal{B}}^{(j,k)}$ has abscissae ranging in the interval 
$$
\left[
\min (  {\tilde{x}}(T^{(j,k)}_i),  {\tilde{x}}(T^{(j,k)}_{i+1})  ),
\max (  {\tilde{x}}(T^{(j,k)}_i),  {\tilde{x}}(T^{(j,k)}_{i+1})  )
\right]
. $$ 

If we partition $[T^{(j,k)}_i , T^{(j,k)}_{i+1}]$  as $\cup_{s=1}^{n_{i,j,k}-1} [T^{(j,k)}_{i,s},T^{(j,k)}_{i,s+1}]$
where
$$
T^{(j,k)}_i=T^{(j,k)}_{i,1} < T^{(j,k)}_{i,2} < \ldots T^{(j,k)}_{i,n_{i,j,k}}=T^{(j,k)}_{i+1}
$$ 
and define the monotone boxes
$$
{\cal{B}}^{(j,k)}_{i,s}:=[T^{(j,k)}_{i,s},T^{(j,k)}_{i,s+1}] \times [\min_{t \in [T^{(j,k)}_{i,s},T^{(j,k)}_{i,s+1}]} {\tilde{y}}(t), \max_{t \in [T^{(j,k)}_{i,s},T^{(j,k)}_{i,s+1}]} {\tilde{y}}(t)]
$$
we have again that $\cup {\cal{B}}^{(j,k)}_{i,s}  \subseteq {\cal{B}}$ contains $\cal S$, but being usually a strict subset of ${\cal{B}}$, in general the solution of an inferior number of polynomial equations is needed to establish what points of $X$ belong to $\cal{S}$.
In passing, we observe that there is numerical evidence that a too high number of boxes slows down the {\it{in domain}} process.  
}
\end{remark}

\begin{figure}[!htbp]
  \centering
   {\includegraphics[scale=0.25,clip]{monotone_boxesII.png}}
 \caption{A domain $\cal S$ of (spline) curvilinear type and a {\it{thinner}} set of {{monotone boxes}} (w.r.t. the previous figure).}
 \label{fig_1}
 \end{figure}
 

\begin{remark} In the algorithm, to see if a point $P=(P_x,P_y)$ belongs to the boundary, we solved a {\it{polynomial}} equation ${\tilde{x}}(t^*)=P_x$ and then tested that 
${\tilde{y}}(t^*)=P_y$. In view of numerical errors, we can only establish that a point is very close to the boundary, that is ${| {\tilde{y}}(t^*)-P_y|}$ is below a certain threshold.  
\end{remark}
\begin{remark} In the case each degree $\delta_k$ is equal to $1$, i.e. $\cal{S}$ is a polygon, it is more profitable to use to the well known {\it{Hormann-Agathos}} algorithm, implemented in Matlab by the {\tt{inpolygon}} routine. 
\end{remark}

\section{Caratheodory-Tchakaloff subsampling}

The purpose of this section is to show how from a rule with high cardinality, positive weights, and $ADE=\delta$ we can extract another with the same degree of precision, positive weights, but with cardinality equal at most to the dimension $N_{\delta}$ of the vector space ${\mathbb P}_{\delta}$.

In some recent papers, the authors have applied a mathematical tool named CATCH (acronym of (Caratheodory-Tchakaloff) subsampling), for such compression of discrete measures, proposing its application to discrete polynomial Least squares by sparse moment matching.  In this framework, the method selects from a large discretization of a given region a much smaller number of (weighted) sampling points, even on a complex shape as can be the polygons investigated in this paper, keeping numerically invariant the Least Squares approximation estimates.

The key point is the following discrete version of Tchakaloff theorem \cite{T57} proved by Caratheodory theorem on finite-dimensional conic/convex combinations \cite{C07}, \cite{PSV17},

\vspace{0.3cm}

\begin{theorem}
Let $\mu$ be a multivariate measure whose support is a ${\mathbb{P}}_{\delta}$-determining finite set $X = \{\xi_i\} \subset {\mathbb{R}}^d$ (i.e., $n$-degree polynomials vanishing there vanish ev- erywhere), with correspondent positive weights (masses) ${\bf{\lambda}} = \{\lambda_i\}$, $i = 1,\ldots,M$, $M = card(X) > N_{\delta}$, $N_{\delta}$ being the dimension of ${\mathbb{P}}_{\delta}$.
Then, there exist a cubature formula for the discrete measure $\mu$, with nodes $T_n = \{t_j\} \subset X$ and positive weights ${\bf{w}} = \{w_i\}$, $1 \leq j \leq m$, with $m \leq N_{\delta}$, such that
$$
\int_X p(x) d\mu = \sum_{i=1}^M \lambda_i p(x_i) = \sum_{j=1}^m w_j p(t_j), {\mbox{  }} \forall p \in {\mathbb{P}}_{\delta}.
$$
\end{theorem}

From the numerical side, given any polynomial basis $\{\phi_k\} $ of $ {\mathbb{P}}_{\delta}$,  define the Vandermonde-like matrix $V=\{V_n(X)\}_{i,j}=\phi_j(\xi_i)$  and let $\gamma$ the vector of moments of the polynomial basis $\{\phi_j\}$ with respect to the original discrete measure and consider the underdetermined moment system $V^Tu=\gamma$.

Caratheodory/Tchakaloff theorem asserts that there exists a sparse nonnegative solution $u$ to the system above, whose nonvanishing components (i.e., the weights $\{w_j \}$) are at most $N_{\delta}$ and determines the corresponding reduced sampling points $T_\delta = \{t_j\} \subseteq X$, that we may term the Caratheodory-Tchakaloff (CATCH) points of X.

The computation of these nodes has been considered in several papers (see e.g. \cite{SV15},  \cite{T15}). To our knowledge, essentially two approaches have been developed to get these rules, i.e. via Linear Programming (LP) and Quadratic Programming (QP).  

About the LP approach, it consists in solving via simplex-method
\begin{equation}
\left\{
\begin{array}{l}
\min_{u\geq 0} c^T u \\
V^Tu=b, {\hspace{0.5cm}} u \geq 0
\end{array}
\right.
\end{equation}

\noindent
where the constraints identify a polytope (the feasible region) in ${\mathbb{R}}^M$ and the vector $c$ is chosen to be linearly independent from the rows of $V^T$, so that the objective functional is not constant on the polytope \cite{PSV17}, \cite{T15}.

The QP based algorithm requires instead the solution of the NonNegative Least Squares (NNLS) problem
\begin{equation}\label{4.2}
{\mbox {compute }} u^* : \|Au^* - \gamma\|_2 =\min_u \|Au - \gamma\|_2 , u \geq 0, 
\end{equation}
in which $u^*$ can be obtained by the well-known {\sl{Lawson-Hanson}} active set optimization method {\cite{LH95}}, which determines a sparse solution to ({\ref{4.2}}). Its application gives a residual $\epsilon = \|Au^*- \gamma\|_2$ that is typically very small, say $< 10^{-14}$ for $\delta \leq 30$. 

Our numerical experience with now available Matlab software has shown that NNLS usually performs better than LP
in computing the CATCH weights, at least for moderate degrees $\delta$ (namely, $N_{\delta}$
for mild degrees) \cite{PSV17}. Consequently all our codes are based on this the application of Lawson-Hanson method to compute the Caratheodory/Tchakaloff sets.

We point out that there are several versions of NNLS codes available in Matlab. One is the built-in function {\tt{lsqnonneg}}, based on the Lawson-Hanson algorithm while an open-source version present in the package NNLSlab in \cite{SH}. Other alternatives are often obtained by MEX files and for generality purposed will not be used here.

\begin{remark}  As alternative, one can compress the rule via the QR algorithm with column pivoting proposed in 1965 by Businger
and Golub {\cite{BG65}}, which is implemented for example by the Matlab {\sl{backslash}}
operator. This approach is equivalent to a greedy selection of the columns of $A$ in order to
maximize the successive volumes, and eventually the (absolute value of the)
determinant (the column selection problem being NP-hard). If the
matrix has full rank, the final result is a weight vector $w^*$ where only $N \leq N_{\delta}=(\delta+1)(\delta+2)/2$ components are nonzero, so that we can extract the cubature nodes $T_\delta = \{t_j\}$ from $X$ by
the column indexes corresponding to such components. A drawback is that the
resulting weights $w^*=\{w^*_k\}$ may not be all positive, but typically the
negative ones are few and of small size, so that the relevant stability parameter
$$
{\mbox{cond}}(\{w_k^*\})=\frac{\sum_{k=1}^N |w^*_k|}{|\sum_{k=1}^N w^*_k|} \geq 1
$$
is not far from $1$ (see also {\cite{SV15}}).
\end{remark}

\begin{remark} In the cubature framework an algorithm termed {\sl{Recursive
Halving Forest}}, based on a hierarchical SVD, has been proposed in \cite{T15}. Performances are reported for large scale problems (say that the order of $N_{\delta}$ is  $10^3$, $10^4$).  Unfortunately the software is not available and thus cannot be applied here as comparison.
\end{remark}

\section{Algebraic cubature on curvilinear domains}

Let $\cal{S}$ be a simple and simply connected domain and suppose that the boundary is defined parametric by splines $\tilde x$, $\tilde y$, as required in the previous chapter.

In this section we describe an algorithm that computes an algebraic cubature formula of degree $n$ over $\cal{S}$, i.e. such that
$$
\int_{\cal S} p_n (x,y) d {\cal S}=\sum_{i=1}^{N_n} w_i p_n(x_i,y_i)
$$
for any bivariate polynomial $p_n$ of total degree at most $n$.

In particular the nodes $(x_i,y_i)$ will belong to $\cal S$, are at most $(n+1)(n+2)/2$ and the respective weights $ w_i$ are all positive i.e., using a common acronym, the rule is of $PI$ type that stands for {\it{positive}} weights and points {\it{inside}} the domain.

We proceed as follows. 
\begin{enumerate}
\item First we compute the {\it{moments}} $\gamma_i$ of a certain polynomial basis $\{\phi_i\}={i=1,\ldots,(n+1)(n+2)/2}$ of degree $n$ over $\cal{S}$, i.e.
$$
\gamma_i=\int_{\cal S} \phi_i (x,y) d {\cal S}.
$$ 
These values can be efficiently obtained by means of Gauss-Green theorem without the need of numerical rules over $\cal{S}$.

In our Matlab software as polynomial basis $\{\phi_j\}$ we have used the total-degree product Chebyshev basis $\{T_p(\alpha_1(x)) T_q(\alpha_2(y))\}$, $(x,y) \in [a_1,b_1] \times [a_2,b_2]$ (the smallest cartesian rectangle containing the domain, easily available when the all the monotone boxes are determined), with $T_h(\cdot)=\cos(h \arccos(\cdot))$ is the $h$-degree Chebyshev polynomial and $\alpha_i(s)=(2s-b_i-a_i)/(b_i-a_i)$, $s \in [a_i,b_i]$, $i=1,2$.

By Gauss-Green theorem (see e.g. \cite{A69}), 
\begin{eqnarray}
\gamma_i&=&\int_{\cal S} \phi_i (x,y) dx dy=\oint_{\partial \cal S} \Phi_i (x,y) dy
\end{eqnarray}
where, for some $p,q$, 
\begin{eqnarray}
\Phi_i (x,y) &=& \int \phi_i(x,y) dx = T_q(\alpha_2(y)) \int T_p(\alpha_1(x)) dx \nonumber
\end{eqnarray}
In particular,  for $p=0$
$$
 \int T_p(\alpha_1(x)) dx = x,
$$
for $p=1$
$$
 \int T_p(\alpha_1(x)) dx = \frac{b_1-a_1}{4} \cdot \alpha^2_1(x),
$$
while for $p \geq 2$
$$
 \int T_p(\alpha_1(x)) dx= \frac{b_1-a_1}{2} \cdot \left( \frac{p}{p^2-1}  T_{p+1}(\alpha_1(x)) - \frac{x}{p-1} T_{p}(\alpha_1(x))\right) .
$$
Setting $I^{(k)}_j=[t^{(k)}_j,t^{(k)}_{j+1}]$, $P_{k,j}=(\tilde{x}(t^{(k)}_j), \tilde{y}(t^{(k)}_j)  )$ and denoting with $P_{k,j} \frown P_{k,j+1}$ the arc of $\partial {\cal S}$ joining $P_{k,j}$ with $P_{k,j+1}$, we have
\begin{eqnarray}
\gamma_i=\oint_{\partial \cal S} \Phi_i (x,y) dy &=& \sum_{k,j} \int_{P_{k,j} \frown P_{k,j+1} } \Phi_i (x,y) dy \nonumber \\
&=& \sum_{k,j} \int_{t^{(k)}_j}^{t^{(k)}_{j+1}} \Phi_i ({\tilde{x}}(t),{\tilde{y}}(t)) \,\, {\tilde{y}}^{\prime}(t) dt.
\end{eqnarray}
Observing that
\begin{itemize}
\item since $\phi_i$ is a polynomial of total degree $n$ then $\Phi$ is a polynomial of total degree $n+1$,
\item in the interval $I^{(k)}_j$ both ${\tilde{x}}$, ${\tilde{y}}$ are polynomials of degree $\delta_k$,
\end{itemize}
each integrand is a polynomial of degree $(n+1)\delta_k+\delta_k-1=(n+2)\delta_k-1$ that can be exactly integrated by a (shifted) Gauss-Legendre formula with $\lceil \frac{(n+2)\delta_k}{2} \rceil$ points.

\item At the $j$-th iteration of the algorithm we define a set ${\cal{P}}_j$, consisting in the union of ${\cal{P}}_{j-1}$ with some new $M_{j}$ points that belong to $\cal{S}$ to ${\cal{P}}_{j-1}$. To this purpose, it is fundamental to apply the {\it{in-domain}} algorithm previously defined.
Next, we apply the algorithm that computes the Tchakaloff formula with nodes $\{(x^{(j)}_i,y^{(j)}_i)\}_{i=1,\ldots,N_n}$ and weights $\{w^{(j)}_i\}_{i=1,\ldots,N_n}$, finally testing if the so obtained rule is such that
$$
\gamma_s^{(j)}=\sum_{i=1}^{N_n} w^{(j)}_i \phi_s(x^{(j)}_i,y^{(j)}_i), {\mbox{   }}s=1,\ldots,N_n
$$
well approximates $\gamma_s$, i.e.
\begin{equation}\label{eqtoll}
\|\gamma_s^{(k)}-\gamma_s\|_{2} \leq {\tt{toll}}
\end{equation}
where {\tt{toll}} is the tolerance fixed by the user. If (\ref{eqtoll}) does not hold we iterate the procedure until (\ref{eqtoll}) is satisfied or a maximum number of iterations is reached, providing in this case an error message.

In our implementation of the Caratheodory-Tchakaloff compression we used the Matlab built-in routine {\tt{lsqnonneg}} but one can alternatively use the algorithm proposed in \cite{SH}.
\end{enumerate}

{\color{blue}
The fact that this procedure terminates in a finite number of iterations comes from a theorem proved in {\cite{W}} and that we state below. 

{\vspace{0.1cm}}

{\em{Let $\Phi$ be the linear span of continuous, real-valued, linearly indepentent functions $\{\phi_k\}_{k=1,\ldots,n}$ defined on a compact set $D$. Assume $\Phi$ satisfies the Krein condition (i.e. there is at least one $f \in \Phi$ which does not vanish on $D$) and that $L$ is a positive linear functional, i.e. $Lf \geq 0$ whenever $f \geq 0$. If $\{x_i\}_{i=1}^{+\infty}$ is an everywhere dense subset of $D$, then for sufficiently large $m$, the set $\{x_i\}_{i=1}^{m}$ is a Tchakaloff set, i.e.
$$
Lf=\sum_{i=1}^m \lambda_i f(x_i), {\mbox{   }} f \in \Phi
$$
where $\lambda_i > 0$ and $\{x_i\}_{i=1,\ldots,m} \subseteq D$.}}

Setting $D={\cal{S}}$, $Lf=\int_{{\cal{S}}} f \,\,\, d {\cal{S}}$, and $\Phi={\mathbb{P}}_n$, i.e. the space of bivariate polynomials of total degree $n$, this theorem implies that our algorithm determines the desired cubature rule of PI-type after a finite number of iterations, extracting the formula by the Lawson-Hanson algorithm from a set of certain $m$ points of $D$.
}


\begin{remark}
 ({\it{Convergence rate}}). Concerning the approximation quality, if $f$ is a continuous function over $\cal S$, $\mu(S)$ is the measure of $\cal{S}$, and we denote by $I_{\cal{S}}(f)$, $I_{n}(f)$ respectively the integral to be computed and its numerical approximation via this new cubature rule having degree of precision $n$, denoting by $p^*_n$ the best uniform approximation of degree $n$, $\|I_{\cal{S}}\|_{\infty}=\mu({\cal{S}})$, $\|I_{n}\|_{\infty}=\sum_i |w_i|=\mu({\cal{S}})$ (here we exploited the fact that all the weights are positive),
\begin{eqnarray}
|I_{\cal{S}}(f)-I_{n}(f)| & \leq & |I_{\cal{S}}(f)-I_{\cal{S}}(p^*_n)|+|I_{\cal{S}}(p^*_n)-I_{n}(p^*_n)| \nonumber \\
&\leq& \|I_{\cal{S}}\|_{\infty} \|f-p^*_n\|_{\infty}+\|I_{n}\|_{\infty} \|f-p^*_n\|_{\infty}\nonumber \\ &=&(\mu({\cal{S}})+\sum_i |w_i|) E_n(f,{\cal{S}}) \nonumber \\ &=&2 \mu({\cal{S}}) E_n(f,{\cal{S}}) \nonumber
\end{eqnarray}
where $E_n(f,{\cal{S}})$ is the best approximation error at degree $n$ over ${\cal{S}}$. We observe that with respect to the error of the algorithm in \cite{SV09} these results offer in general a tighter estimate.

In particular by multivariate versions of Jackson theorem, we have
$$
E_n(f,{\cal{S}})={\cal{O}}(n^{-(p+\theta)}), \,\, f \in C^{p+\theta}({\cal{S}})
$$
for every function $f$ that has $\theta$-H\"older continuous $p$-th partial derivatives, with $\theta \in (0,1]$.
\end{remark}

\begin{remark} ({\it{Boundary Approximation}}). As observed in \cite{SV09}, $\cal{S}$ can be an approximation of a certain domain $\Omega$. 

As stated before, assume that the vertices $V_k \in \partial \Omega$, $k=1,\ldots,M+1$, are computed and that $\partial \Omega \approx \partial {\cal{S}}:= \cup_{k=1}^M V_k \frown V_{k+1}$ where each curvilinear side $V_k \frown V_{k+1}$ is tracked by a spline curve of degree $\delta_k \leq 3$, interpolating an ordered subsequence of control points $P_{1,k}=V_k,P_{2,k},\ldots,P_{m_k,k}=V_{k+1}$.

Furthermore, suppose that it is adopted a {\it{cumulative chordal}} parametrization, i.e. each $V_k=({\tilde{x}}(t^{(k)}),{\tilde{y}}(t^{(k)}) )$ so that $I^{(k)}=[t^{(k)},t^{(k+1)}]$, being 
$$
t^{(k+1)}=t^{(k)}+\sum_{j=1}^{m_k-1}|P_{j+1,k}-P_{j,k}|.
$$

In these assumptions, also supposing that the boundary has a global $C^4$ pa\-ra\-me\-trization, that $H=\max_{j,k}|P_{j+1,k}-P_{j,k}|$, and setting ${\Omega \Delta {\cal{S}}}=({\Omega \backslash {\cal{S}}}) \cup ({{\cal{S}} \backslash \Omega })$, with an approach similar to that adopted in {\cite{SV09}}, by a result in {\cite{24}}, we get
\begin{eqnarray}
|I_{\Omega}(f)-I_n(f)| &\leq& I_{\Omega \Delta {\cal{S}}} (|f|)  + |I_{{\cal{S}}}(f)-I_n(f)| \nonumber \\
&\leq& \|f\|_{\infty,{\Omega \Delta {\cal{S}}}} \,\, \mu({\Omega \Delta {\cal{S}}})+2 \mu({\cal{S}}) E_n(f,{\cal{S}}) \nonumber \\
&\approx& \|f\|_{\infty,{\Omega \Delta {\cal{S}}}} \,\, l(\partial \Omega) O(H^4) + 2 \mu(\Omega)E_n(f,\Omega) \nonumber.
\end{eqnarray}
where $E_n(f,{\cal{S}})$ is the best polynomial approximation error of degree $n$ of the function $f$, over $\cal S$ (notice that differently from {\cite{SV09}} all the nodes are inside the domain ${\cal S}$ and the weights are positive).

\end{remark}
\begin{remark} ({\it{Moments approximation}}). 
Let $\{\phi_k\}$ be a basis of the space of polynomials of total degree $n$, and $\gamma_k=\int_{\Omega} \phi_k d{ \Omega}$, being $\Omega$ a general domain. In our examples $\Omega$ corresponds to $\cal{S}$, but it can be a more general compact domain. Suppose that the formula has nodes $(x_i,y_i)$, weights $w_i$, $i=1,\ldots,N_n$ and is of PI-type. If ${\tilde{\gamma}}_k=\sum_{i=1}^{N_n} w_i \phi_k(x_i,y_i)$ and $\epsilon_{\tiny{mom}}=\|{\tilde{\gamma}}_k-{{\gamma}}_k\|_2$, then, for $C \leq 2(\mu({\Omega}) + \sqrt{\mu{({\Omega})}})$, we have
$$
|\int_{\Omega} f(x,y) d \Omega-\sum_{i=1}^{N_n} w_i f(x_i,y_i)| \leq C E_n(f, \Omega) + \|f\|_{{\cal L}_2(\Omega)} \epsilon_{\tiny{mom}}
$$
for any $f \in C(\Omega)$ (see {\cite{SV15}} for details, using the Lebesgue measure $d\lambda$ instead of the product Chebyshev measure on the smallest rectangle containing $\Omega$, having sides parallel to $x$ and $y$ axes).
\end{remark}

\section{Numerical experiments}

The purpose of this section is twicefold, on a side we test the new cubature algorithm and on the other the in-domain routines.

We observe the Matlab routine {\tt{splinegauss}} proposed in \cite{SV09} is not necessarily of $PI$-type on non convex domains $\cal S$, since in these instances some nodes could be outside the domain and some weights could be even negative.

However, since it was based on a sequential use of tensorial Gauss-Legendre formula, the determination of the cubature rule was rather fast, w.r.t. the geometry of the domain and the degree of precision required.

For mild degrees, say $n \leq 15$, the new routine {\tt{splcub}}, implementing the ideas developed in the previous sections, has similar cputimes but produces rules of PI-type with a number of nodes $N_n \leq (n+1)(n+2)/2$.

Being all the weights positive, these new formulas are optimally stable, i.e. the cubature condition number
$$
{\mbox{cond}}(\{w_i\})=\frac{\sum_{i=1}^{N_n}|w_i|}{\sum_{i=1}^{N_n} w_i}
$$
is equal to $1$, the best possible result for rules with ADE at least equal to $0$.
Furthermore, they are appropriate for problems were the sampling of the integrand is not possible outside the domain. 

We implemented an algorithm that can be sketched as follows:
\begin{enumerate} 
\item define a mesh on the smallest rectangle containing the domain $\cal S$, say $[a_1,b_1] \times [a_2,b_2]$; in particular if $ADE=n$, then setting $k=\lfloor n^{1.5}\rfloor$, we considered the points $P=((P_x)_i,(P_y)_j)$ where
$$
(P_x)_i=a_1+i \frac{a_2-a_1}{k-1}, (P_y)_j=b_1+j \frac{b_2-b_1}{k-1}
$$
i.e. a tensor grid mesh, based on $k$ equispaced data in each direction. 
\item establish those $P$ strictly inside the domain, say ${\cal{P}}^*=\{P_1^*,\ldots,P_m^*\}$;
\item compute the moments $\gamma_k$ over $\cal{S}$ of the total-degree product Chebyshev basis $$\{T_p(\alpha_1(x)) T_q(\alpha_2(y))\}, {\mbox{ }}(x,y) \in [a_1,b_1] \times [a_2,b_2],$$ by means of Gauss-Green theorem,
\item extract the nodes and the weights using the technique introduced in {\cite{SV15}} and described in section II.
\end{enumerate} 

In case of failure, i.e. the norm 2 of the moment error is bigger than a fixed tolerance, say $10^{-12}$, we proceed by defining a finer grid, determining those points in ${\cal S}$, providing the union with the previous $P_1^*,\ldots,P_m^*$, and retrying points 3 and 4 of the algorithm above.

Due to the fact that for higher $ADE$, ill-conditioning of Vandermonde matrix $V$ on ${\cal{P}}^*$ may occur,  whenever it was necessary, we have reothortogonalised the basis by the {\it{economy size}} QR algorithm.


All the Matlab routines for the numerical experiments are available at {\cite{SH}}.

We start our numerical test by analysing the region ${\cal{S}}_1$, whose boundary is defined by $\tilde x$, $\tilde y$, that in $I^{(1)}=\{I_j^{(1)}\}_{j=1,\ldots,4}$ are linear splines, i.e. $\delta_1=2$, while in $I^{(2)}=\{I_j^{(2)}\}_{j=1,\ldots,5}$ are cubic splines satisfying {\it{not-a-knot}} conditions. The region is depicted in Figure {\ref{fig_1N}} and it consists of a polygon with a side substituted by a curve.

\begin{figure}[!htbp]
  \centering
   {\includegraphics[height=5cm]{dom1N6.png}}
   {\includegraphics[height=5.15cm]{dom2N6.png}}
 \caption{\footnotesize{The spline curvilinear domains ${\cal{S}}_i$ with $i=1,2$, the grid points $P$ outside the domain (in red), those inside the domain (in green) and the nodes of a cubature formula of PI-type for $n=6$ (black dots).}}
 \label{fig_1N}
 \end{figure}


In particular if the routine in {\cite{SV09}} is applied to these tests, several nodes are external to the domain, and negative weights are present.

 
 
 
\begin{table}[ht]
\begin{center}
\label{tab01}
{\renewcommand{\arraystretch}{1.1}
\begin{tabular}{|c| c c | c  | c | }
\hline
{\mbox{ADE}} & $\#$ & $\#_{O}$ & {\mbox{cond}} & {\mbox{res}} \\
\hline
 	  $ 5 $ & $  21$ &  $0$  & $1$ &  $7.2e-17$ \\
 	  $10 $ & $ 66$ & $0$  & $1$ &   $1.3e-16$ \\
 	  $15 $ & $136$ & $0$ & $1$ &   $1.1e-16$ \\
 	  $20 $ & $231$ & $0$ & $1$ &   $1.2e-16$ \\
\hline
\end{tabular}
        }
\caption{\footnotesize{Cardinality $\#$ of the nodes, number of nodes $\#_{O}$ outside the domain, conditioning of the rule and residual of the rule w.r.t. the moments $\{\gamma_k\}$ of domain ${\cal{S}}_1$.}}
\end{center}
\end{table}




 
 
 
 
 
\begin{table}[ht]
\begin{center}
{\renewcommand{\arraystretch}{1.1}
\begin{tabular}{|c| c c | c  | c | }
\hline
{\mbox{ADE}} & $\#$ & $\#_{O}$ & {\mbox{cond}} & {\mbox{res}} \\
\hline
 	  $ 5 $ & $  21$ &  $0$  & $1$ &  $6.7e-18$ \\
 	  $10 $ & $ 66$ & $0$  & $1$ &   $3.8e-18$ \\
 	  $15 $ & $136$ & $0$ & $1$ &   $3.9e-18$ \\
 	  $20 $ & $231$ & $0$ & $1$ &   $4.8e-18$ \\
\hline
\end{tabular}
        }
\caption{\footnotesize{Cardinality $\#$ of the nodes, number of nodes $\#_{O}$ outside the domain, conditioning of the rule and residual of the rule w.r.t. the moments $\{\gamma_k\}$ of the hand-like domain ${\cal{S}}_2$.}}
\label{tab02}
\end{center}
\end{table}








It not surprising that, with the exception of the residuals, the results are identical, in view of the fact that we have computed rules with at most $N_n=(n+1)(n+2)/2$ nodes.
All the rules are of $PI$-type, and the residuals are particularly small.

As second domain ${\cal{S}}_2$, we considered a hand-like domain, whose boundary is defined by $\tilde x$, $\tilde y$, that in $I^{(1)}=\{I_k^{(1)}\}_{k=1,\ldots,37}$ are both cubic splines satisfying {\it{not-a-knot}} conditions.

\begin{remark} (Nested rules). We point out that in both the examples, as shown in \cite{SV18}, if it is available a formula of PI-type, with ADE equal to $n$, then one can extract those with ADE less than $n$, i.e. {\it{nested rules}}, by further application of Caratheodory-Tchakaloff compression, with the advantage of having estimates of the cubature error, without further function evaluations. 
 \end{remark}
 
 
As final tests, we show the cputimes necessary to process $m$ in-domain operations, on the regions ${\cal{S}}_i$, $i=1,2$. To this purpose we have defined the Matlab routine {\tt{incurvpolygon}}, that implements the ideas introduced in section II.






\begin{table}[ht]
\begin{center}
{\renewcommand{\arraystretch}{1.1}
\begin{tabular}{|c| c | c |}
\hline
  m   &  ${\cal{S}}^{(1)}$   & ${\cal{S}}^{(2)}$  \\
  \hline
  $10^2$ & $6e-03$ & $1e-02$ \\
    $10^3$ & $1e-02$ & $3e-02$ \\
      $10^4$ & $4e-02$ & $7e-02$ \\
        $10^5$ & $6e-01$ & $6e-01$ \\
        \hline
\end{tabular}
        }
\caption{\small{Cputime of the application of the {\it{in-domain}} algorithm for $m$ points, respectively in ${\cal{S}}^{(1)}$ and ${\cal{S}}^{(2)}$.}}
\label{tab03}
\end{center}
\end{table}







The results are particularly positive, and allow the application of the cubature algorithm at higher degrees. In this case the bottleneck is the application of the Lawson-Hanson algorithm. If negative weights are not an issue, the alternative provided by QR method provides better cputimes and formulae that usually have quite good cubature condition number.


\section{Acknowledgements}

This work is partially supported by the DOR funds of the University of Padova, and by the GNCS-INdAM. Their research has been accomplished within the RITA {\sl{Research ITalian network on Approximation}}.
 
\bibliographystyle{plain}
\begin{thebibliography}{99}

\bibitem{A69} % CITED
 {\sc T.M.~Apostol},, Calculus, 2nd ed., vol. II, Blaisdell, 1969.
 
\bibitem{BSV19} % CITED
 {\sc B.~Bauman, A.~Sommariva and M.~Vianello},  Compressed cubature over polygons with applications to optical design, submitted.

\bibitem{BRV18} % CITED
 {\sc L.~Beir${\tilde{{\mbox{\sc{a}}}}}$o da Veiga, A.~Russo and G.~Vacca},  The Virtual Element Method with curved edges, ESAIM Mathematical Modelling and Numerical Analysis 53 (2019), pp.375--404.
 
 
\bibitem{BG65} % CITED
 {\sc P.A.~Businger, G.H.~Golub}, Linear least-squares solutions by Householder transformations, Numer. Math. 7 (1965), pp.269--276.

\bibitem{C07} % CITED
 {\sc C.~Caratheodory}, {\"U}ber den Variabilittsbereich der Koeffizienten von Potenzreihen, die gegebene Werte nicht annehmen, Math. Ann. 64 (1907), pp.95--115.

\bibitem{D67} % CITED
 {\sc P.J.~Davis}, A construction of nonnegative approximate quadratures, Math. Comp. v.21 (1967), pp.578--582.


%\bibitem{MFLS17} % CITED
%{\sc J.M.~Myre, E.~Frahm, D.J.~Lilja, M.O.~Saar}, TNT-NN: A Fast Active Set Method for Solving %Large Non-Negative Least Squares Problems, Procedia Computer Science, 108 (2017), pp.%755--764. See {\tt{https://github.com/MyreResearchGroup/tnt-nn}} for the related software.





\bibitem{HA01}
 {\sc K.~Hormann and A.~Agathos}, The point in polygon problem for arbitrary polygons, Comput. Geom. Theory Appl. v.20, n.3 (2001), pp. 131--144.


\bibitem{LH95}
 {\sc C.L.~Lawson and R.J.~Hanson}, Solving least squares problems, Classics in Applied Mathematics 15, SIAM, Philadelphia, 1995.


\bibitem{PSV17} % CITED 
{\sc F.~Piazzon, A. ~Sommariva and M. ~Vianello}, Caratheodory-Tchakaloff Subsampling, Dolomites Res. Notes Approx. DRNA 10 (2017), pp.5--14.


\bibitem{24} {\sc M. ~Sakai, R.A.~ Usmani}, On orders of approximation of plane curves by parametric cubic splines, BIT Numerical Mathematics 30 (1990) 735--741.

\bibitem{SSV11} % CITED 
{\sc G.~Santin, A. ~Sommariva and M. ~Vianello}, An algebraic cubature formula on curvilinear polygons, Appl. Math. Comput. 217 (2011), pp.10003--10015.


\bibitem{SH} % CITED 
{\sc M.~Slawski}, Non-negative least squares: comparison of algorithms, {\tt{https://sites.google.com/site/slawskimartin/code}}.

\bibitem{SH} % CITED
 {\sc A.~Sommariva}, {\tt{http://www.math.unipd.it/$\sim$alvise/software.html}}.
 
\bibitem{SV07} % CITED
 {\sc A.~Sommariva and M.~Vianello}, Product Gauss cubature over polygons based on Green's integration formula, BIT Numerical Mathematics, 47 (2007), pp. 441--453.
 
\bibitem{SV09} % CITED
 {\sc A.~Sommariva and M.~Vianello}, Gauss-Green cubature and moment computation over arbitrary geometries, Journal of Computational and Applied Mathematics, 231 (2009), pp. 886--896.
 
\bibitem{SV15} % CITED
 {\sc A.~Sommariva and M.~Vianello}, Compression of multivariate discrete measures and applications, Numer. Funct. Anal. Optim. 36 (2015), pp.1198--1223.
 
 \bibitem{SV18} % CITED
 {\sc A.~Sommariva and M.~Vianello}, Nearly optimal nested sensors location for polynomial regression on complex geometries, to appear on Sampling Theory in Signal and Image Processing.
 

\bibitem{T57} % CITED
 {\sc V.~Tchakaloff}, Formules de cubatures m\'ecaniques \`a coefficients non n\'egatifs, (French) Bull. Sci. Math. 81 (1957), pp.123--134.
 
 \bibitem{T15} % CITED
 {\sc M.~Tchernychova}, {\sl{Caratheodory}} cubature measures, Ph.D. dissertation in Mathematics (supervisor: T. Lyons), University of Oxford, 2015.

\bibitem{W}
{\sc Don.R.~Wilhelmsen}, A Nearest Point Algorithm for Convex Polyhedral Cones and Applications to Positive Linear approximation, Math. Comp. vol. 30, n. 133, January 1976, pp.48--57-


\end{thebibliography}
\end{document}