---
title: "Fitting Model with constraints"
author: "D. Renard"
date: "26 avril 2022"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
rm(list=ls())
library(RGeostats)
```

# Introdution

This is an answer to a post from J. Yarus looking for information on the way to fit a Model **with constraints**.

# Creating data

We first create a data set using simulation on a regular grid (150 by 100). The model is a nested variogram with nugget effect (sill=1), an exponential structure (sill=2 and range=30) and a cubic variogram (sill=1 and range=60)

```{r}
grid = db.create(nx=c(150,100))
print(grid)

model = model.create(vartype="Nugget Effect",sill=1)
model = model.create(vartype="Exponential",range=30,sill=2,model=model)
model = model.create(vartype="Cubic",range=60,sill=1,model=model)
print(model)
```

We perform the (non-conditional) simulation over the grid

```{r}
grid = simtub(,grid,model,nbtuba=200)
grid = db.rename(grid,grid$natt,"Simu")
plot(grid)
```

# Caculate experimental variogram

We calculate the experimental variogram along the two main directions of the grid, benefiting from the grid organization of the data (i.e. using **vario.grid** instead of **vario.calc**).

```{r}
vario = vario.grid(grid,nlag=80)
plot(vario)
```

# Fitting under constraints

Now we wish to demonstrate the capacities of the fitting procedure, i.e. perform fitting with constraints on the various basic structures

## Natural fitting

We first let the system perform the fitting automatically

```{r}
model = model.auto(vario)
model
```

The resulting model is composed of a nugget effect (sill=0.998), a cubic structure (sill=1.3 and range of 51.7 along X and 57.2 along Y) and a spherical (sill=0.766 and range of 68 along X and 45 along Y).

## Constrained fitting

This paragraph gives indications on the way to perform the fit adding constraints of the different ingredients of the resulting model.

As a first remark, note that the constraints grammar covers a much larger field than the one of fitting a single model for a single variable. It can be used when fitting simultaneously several Gaussian Random Functions for several variables and using several basic structures, hence the complexity of the terminology.

The best way to introduce this grammar is to consider tha the user wishes to obtain a model composed of pre-defined basic structures: for example a nugget effect and an isotropic cubic.

```{r}
model.auto(vario,struct=c("Nugget Effect","Cubic"))
```

The result is a model composed of the nugget effect (sill=1.028) and a cubic (sill=1.912 with range 51.525 along X and 48.52 along Y).

Note that the number of structures in the resulting model is not predictable. As a matter of fact, this procedure is allowed to suppress structures whose sill would be too small with respect to the total sill. 
This is the case in the next example where a Gaussian structure has been added. As the Gaussian (with a sill of 1.908) and the Cubic variograms are quite similar, the system prefers the Gaussian and drops out the Cubic.

```{r}
model.auto(vario,struct=c("Nugget Effect","Cubic","Gaussian"))
```

If we absolutely want to keep all the provided basic structures, the option **flag.noreduce** can be added. Then the system will balance the continuous part of the structure between the Cubic (sill=0.092) and the Gaussian (sill=1.937)

```{r}
model.auto(vario,struct=c("Nugget Effect","Cubic","Gaussian"),flag.noreduce = TRUE)
```

This remark is important as we will have to identify each structure on which we want to define a constrain by its rank.

## Add ing a constraint as an interval

Let us return to the case of a combination of Nugget Effect and Cubic (we add  **flag.noreduce** as a safety).

```{r}
model.auto(vario,struct=c("Nugget Effect","Cubic"),flag.noreduce=TRUE)
```

Let us imagine that we want to constraint the sill of the Cubic structure to be larger than 2. We refer to the lower bound ("keyword *lower*) of the second structure (*M2*) and to the sill of the first variable (*V1*)

```{r}
model.auto(vario,struct=c("Nugget Effect","Cubic"),flag.noreduce=TRUE,
           lower = c("M2V1=2"))
```

As expeccted, the sill of the Cubic structure is equal to 2. The other parameters are fitted without any constraints.

We can use the same technique in order to constraint the sill of the nugget effect to lie within the interval [0.2; 0.9]

```{r}
model.auto(vario,struct=c("Nugget Effect","Cubic"),flag.noreduce=TRUE,
           lower = c("M1V1=0.2"), upper = c("M1V1=0.9"))
```

The simple c() simply indicates that the constraints can be stacked on the different structures. If they only concern a single constraint on a single structure, the previous line can be written for short

```{r}
model.auto(vario,struct=c("Nugget Effect","Cubic"),flag.noreduce=TRUE,
           lower = "M1V1=0.2", upper = "M1V1=0.9")
```

## Constraints defined as equality

We can also ut a constraint as an equality. For example, we wish to set the sill of the nugget effect equal to 1.2

```{r}
model.auto(vario,struct=c("Nugget Effect","Cubic"),flag.noreduce=TRUE,
           equal = c("M1V1=1.2"))
```

## Constraints on the total sill

Another way to put a constraint on the total sill of the model (say the total is set to 3.2).

```{r}
model.auto(vario,struct=c("Nugget Effect","Cubic"),constraints=3.2)
```