(sampler_stochastic_method_constrained_lhs_heuristic)=
# The heuristic

The idea behing this empirical heuristic is to rely on permutations and to decide on the best 
permutation to be done thanks to the content of the constraint matrix and also the distributions of 
solutions along the row and columns. This will be discussed further in the rest of this section but 
first a focus is done on the definition of a constraint. If one considerers a simple constraint, its 
implementation can be decompose in the following steps:

1. define it simply with an equation, for instance one wants to reject all locations for which 
$x_0 < x_1$;

2. define a constraint function that can compute the margin of success, for instance in our simple 
case $c(x_0, x_1) = x_0 - x_1$;

3. define a characteristic constraint function that only states whether the constraint is fulfilled, 
for instance in our simple case $ {\rm 1\kern-0.28emI}_{c}(x_0, x_1) = \begin{cases} 0 & \text{if} 
x_0 < x_1 \\ 1 & \text{if} x_0 > x_1 \end{cases} $

If formally the constraint function can provide more information when reading it (in terms of margin), 
one needs to know the way to apply a selection on these results which is not the simplest aspect to 
provide which explains why in the rest of this section the focus will be put on the characteristic 
contraint function. In order to illustrate our method, one can start from a provided LHS {{doe}}, 
called hereafter $\mathcal{L} = \{ \mathbf{x}^{i}, i=1,\ldots,n_S\}$ where $\mathbf{x}^{i}$ is the 
i-Th location which can be written as $\mathbf{x}^{i}=(x^{i}_1, \ldots ,x^{i}_{n_X})$: it is a sample 
of size $n_S$ in the $n_X$ input space. From there, the constraint matrix is indeed defined as done 
below: 

```{math}
\mathbf{C}(n_{S},n_{S})= \left ( \begin{array}{cccc} {\rm 1\kern-0.28emI}_{c}(x^1_{row}, x^1_{col}) & 
{\rm 1\kern-0.28emI}_{c}(x_{row}^{1}, x_{col}^{2}) & \cdots & {\rm 1\kern-0.28emI}_{c}(x_{row}^{1}, 
x_{col}^{n_{S}})\\ {\rm 1\kern-0.28emI}_{c}(x_{row}^{2}, x_{col}^{1}) & {\rm 1\kern-0.28emI}_{c}
(x^2_{row}, x^2_{col}) & \cdots & {\rm 1\kern-0.28emI}_{c}(x_{row}^{2}, x_{col}^{n_{S}})\\ \vdots & 
\vdots & \ddots & \vdots \\ {\rm 1\kern-0.28emI}_{c}(x_{row}^{n_{S}}, x_{col}^{1}) & 
{\rm 1\kern-0.28emI}_{c}(x_{row}^{n_{S}}, x_{col}^{2}) & \cdots & 
{\rm 1\kern-0.28emI}_{c}(x_{row}^{n_{S}}, x_{col}^{n_{S}}) \\ \end{array} \right)
```

The constraint matrix is, as visible from the definition above, a $(n_S,n_S)$ matrix which only 
contains 0 and 1 depending on whether the current configuratio of the $(x_{row},x_{col})$ plane 
fulfill the constraint. This shows why restraining the constraint to only two-variables function is 
a reasonable approach: for a given configuration the number of permutations will blow up if one will 
considerer larger dimension constraint. Giving this object, one can introduces more useful objects:

- the constraint row solution vector, that sums up for every row the number of solutions, *i.e.* the 
number of columns that allow to fulfill the constraint $\lbrace c^i_{row} = \sum_{j=1}^{n_S} 
\mathbf{C}^{ij}, \forall i \in [1, n_S] \rbrace$;

- the constraint column solution vector, that sums up for every column the number of solutions, 
*i.e.* the number of row that allow to fulfill the constraint $\lbrace c^i_{col} = \sum_{i=1}^{n_S} 
\mathbf{C}^{ji}, \forall i \in [1, n_S] \rbrace$;

- the constraint diagonal vector *i.e.* the diagonal of the constraint matrix $\lbrace c^i_{diag}(n_S) 
= \mathbf{C}^{ii}, \forall i \in [1, n_S] \rbrace$.

Bearing in mind that the objective is to have a constraint fulfilled, the heuristic will use the 
following criterion as a stopping signal: $ c_{stop} = \sum_{i=1}^{n_S} c^i_{diag} = n_S$. This 
heuristic starts from a provided LHS {{doe}} for which one can consider three cases:
          <variablelist>

**A**
: it cannot fulfill the constraint, meaning that there is no sets of permutation to have $c_{stop} 
= n_S$;

**B**
: it can fulfill the constraint with only one set of permutation leading to the only solution 
$\mathcal{L}^{constr}$;

**C**
: it can fulfill the constraint with many different sets of permutation leading to a very large 
number of configurations and so to a very large number of constrained {{doe}}.

Practically, the heuristic is organised in a step-by-step approach in which the variable used as 
first argument in the constraint definition will be used as row indicator (it will be called $x_{row}$) 
and it is considered fixed. This implies that the permutation will be done by reorganising the second 
variable, called $x_{col}$. The heuristic is then described below:

1. The constraint matrix is estimated along with all the objects discussed above $(\mathbf{C}, 
c_{row}, c_{col}, c_{diag})$.

    - If $c_{row}$ contains one 0, it means that the {{doe}} cannot fulfill the constaint. If this 
    {{doe}}, has been provided, the method stops, on the other hand if it was generated on-the-fly, 
    then another attempt is done (up to 5 times).

    - A bipartite graph method is called to check that there can be at least one solution for every 
    row. If not, if the {{doe}} has been provided, the method stops, on the other hand if it was 
    generated on-the-fly, then another attempt is done (up to 5 times).

    This step allows to sort out the {{doe}} that will fall into the category A (defined above) from 
    those that might fall either in the B or C ones.

2. If $c_{stop} < n_S$, then all the rows for which the diagonal elements is 0 are kept aside and 
sorted out by increasing number of solutions over all the columns (their value of $c_{row}$), which 
defines $\Omega_{row} = \lbrace i \in [1, n_S], c^i_{diag}=0 \rbrace$. The row considered, whose 
index will be written k for kept hereafter, is the one with the lowest number of solutions (the most 
urgent one in a way), so it can be written $k = \min_{c_{row}} \Omega_{row} $

3. For $x_{row}^{k}$ the chosen value, all the columns that provide a solution are kept aside and 
sorted out by increasing order[^heuristic] of solutions over all rows (their value of $c_{col}$), 
which defines $\Omega^k_{col} = \lbrace i \in [1, n_S], \mathbf{C}^{ki}=1 \rbrace$. This sorting 
provides information on the marging, as the highest values of solutions over the rows means that this 
value of $x_{col}$ is compatible with many other $x_{row}$ instances. A loop is performed over all 
these solutions, so that $\forall t \in \Omega^k_{col}$:

    - By definition we know that $\mathbf{C}^{kk}=0$ and $\mathbf{C}^{kt}=1$, so if $\mathbf{C}^{tk}
    =1$ then the permutation will only increase the stopping criterion ($c_{stop}$), as the actual 
    column $k$, will become the new column $t$ after permutation. In this case, one can move to the 
    permutation step. This can be written, by calling $s$ the actual index of the column (for 
    selected), $s = \min_{c_{col}} \lbrace t \in \Omega^k_{col}, \mathbf{C}^{tk}=1 \rbrace $. From 
    there, one moves to the permutation step.

    - If none of the solution $t$ under investigation can satisfy the requirement $\mathbf{C}^{tk}=1$, 
    then $s = \min_{c_{col}} \lbrace t \in \Omega^k_{col}, t \not\in \Omega^k_{col,perm} \rbrace $. 
    In this definition, the ensemble $\Omega^k_{col,perm} $ is the ensemble of column index which has 
    already been used previously (as this is an interative heuristic) in a permutation for the row 
    $k$ under investigation. This precaution has been introduced in order to prevent from having a 
    loop in the permutation process. From there, one moves to the permutation step.

    - If $ \lbrace t \in \Omega^k_{col}, t \not\in \Omega^k_{col,perm} \rbrace = \varnothing $, 
    meaning that all possible solutions have been tested, then the selected column index is the 
    result of a random drawing in the nsemble of solutions, meaning that 
    $s = {\rm rand}(\Omega^k_{col})$. From there, one moves to the permutation step.

4. Now that both $k$ and $s$ are known, the permutation will be done, it consists in:

    - changing content of the {{doe}}, meaning doing $x_{col}^k \leftrightarrow x_{col}^s$;

    - changing the columns in the constraint matrix, meaning doing $\mathbf{C}^{ik} \leftrightarrow 
    \mathbf{C}^{is}, \forall i \in [1,n_S]$;

    - changing the content of the constraint column solution vector, meaning doing $c_{col}^k 
    \leftrightarrow c_{col}^s$;

    - fill once more the constraint diagonal vector ($c_{diag}$) and compute once more the stopping 
    criterion ($c_{stop}$) to check whether the permutation process has to be continued. If so, then 
    one moves to the second step.

[^heuristic]: the decreasing order has also been tested, but it has show a lower discrepancy in the
resulting {{doe}}.

This presentation has been simplified as, here, there is only one constraint applied to the {{doe}}. 
Indeed, when there are more than one constraint, let's call $(\mathbf{C}, c_{row}, c_{col}, c_{diag})$ 
and $(\mathbf{D}, d_{row}, d_{col}, d_{diag})$ the objects associated to two constraints defined in 
both planes $(x_{row}^c, x_{col}^c)$ and $(x_{row}^d, x_{col}^d)$, few cautions and extra steps have 
to be followed as long as both stopping criteria, $c_{stop}$ or $d_{stop}$ are different from $n_S$ 
depending on the variables in the constraint plane definition:

- if planes $(x_{row}^c, x_{col}^c)$ and $(x_{row}^d, x_{col}^d)$ have no common variable, or if 
$x_{row}^c = x_{row}^d$ but $x_{col}^c \neq x_{col}^d$ then the constraints can be considered 
orthogonal;

- if $x_{col}^c$ is the same variable as $x_{row}^d$ then any permutations for the constraint $c$ 
have to be propagated to the objects of the constraint $d$, meaning that on top of the list in step 
4, the extra steps consist in:

    - changing the rows in the constraint matrix for the constraint $d$, meaning doing $\mathbf{D}^{ki} 
    \leftrightarrow \mathbf{D}^{si}, \forall i \in [1,n_S]$;

    - changing the content of the constraint row solution vector for constraint $d$, meaning doing 
    $d_{row}^k \leftrightarrow d_{row}^s$;

    - fill also the constraint diagonal vector ($d_{diag}$) and compute once more the stopping 
    criterion ($d_{stop}$) to check whether the permutation process has to be continued.

- if $x_{col}^c = x_{col}^d$, then any permutations from on the constraint can undo a previous one 
defined from the other one, meaning that one can create a loop that will (thanks to our heuristic 
definition) will end up into the random permutation area. To prevent this, so far, one can not create 
a second constraint using the same variable as second argument.